The original source locations for the 223 P. ponderosa genotypes in the Chico orchard (Fig. 1) fall within just one of the several genetic subdivisions previously identified in ponderosa pine [66–68]. Fresh needles were collected from these individuals and placed in labeled tea bags over silica gel to quickly dry them and preserve the DNA for extraction.
Seeds were collected from a subset of 50 genotyped parent trees in the summer of 2018. Because pines are wind-pollinated and outcrossing , seeds from the same tree are mostly half-siblings, occasionally full-sibs. We placed 2–3 mature cones from each mother tree into paper bags and placed them in a warm dry place until seeds were released. The seeds were stored in a refrigerator until the greenhouse experiment was carried out (see greenhouse experiment section below). All voucher specimens were maintained at Moran’s Lab in University of California, Merced.
DNA was extracted from the dried needles using a modified Qiagen plant kits protocol and quantified using an Eppendorf BioSpectrometer (Eppendorf, AG, Germany). Samples were frozen and sent to the UC Davis Genome Center for processing. Four 48-plex GBS libraries consisting of 47 DNA samples and a negative control (no DNA) and one 36-plex GBS library consisting of 35 DNA samples and a negative control were prepared. The pool was quantified via qPCR using the KAPA Library Quantification Kit (Kapa Biosystems, Wilmington, MA, USA) for Illumina sequencing platforms, with 0.9X bead cleanup to remove small fragments (< 250 bp). Additional DNA purification using the Zymo DNA Clean & Concentrator kit (Zymo Research, Irvine, CA) was performed to increase the purity of the extracted DNA. The libraries were then sequenced (single-end read 90 bp or 100 bp) using an Illumina HiSeq 4000 (Illumina, San Diego, CA), one library per lane.
No reference genome is available for ponderosa pine (Pinus ponderosa), but one does exist for loblolly pine (Pinus taeda) [69, 70]. Of the conifers that have been sequenced to date, P. taeda is the most closely related to P. ponderosa [71, 72]. Furthermore, the P. taeda reference genome was used to successfully used to design probes for sequence capture in P. contorta [73, 74]. Based on preliminary analyses, we selected the Stack v.2.2 pipeline  with this reference genome (https://treegenesdb.org/FTP/Genomes/Pita/) for SNP calling (Shu & Moran, in review). Each step in the Stacks reference pipeline is performed internally in Stacks algorithms except alignment with BWA v.0.7.17  and the Samtools v.1.9  step used to get read position. Default settings were used in Stacks, BWA and Samtools.
After calling the SNPs, we ran SnpEff  to identify the location of gene that the SNP locates. We built the data base with the annotated genome and the reference genome of loblolly pine v.2.01 in TreeGenes (http://treegenesdb.org/FTP/Genomes/Pita/v2.01/). The location of each SNP is listed in the output file of SnpEff as one of six primary location categories, including intragenic variants, intergenic variants, upstream SNPs, downstream SNPs, synonymous, and missense variants in the gene coding sequence. In Snp Eff, "intragenic" refers to SNPs in introns, while "missense" refers to any non-synonymous mutation in the transcribed region.
Many SNPs identified by GBS fall between genes and regulatory regions (in the intergenic regions) and likely have no direct effect on gene expression or function. In addition, because of the low amount of linkage disequilibrium in conifers [79, 80], any associations identified between such intergenic SNPs and a phenotype or environment of interest are likely false positives, rather than reflecting linkage between the SNP and a causal variant. We therefore filtered out the intergenic SNPs first before running the association analysis using a Python script (https://github.com/shumengjun/LFMM).
We obtained 30-year (1921–1950) averages of climate data for each genotype source location from the 270 m resolution California Basin Characterization Model (BCM) . The five variables were mean climatic water deficit (CWD, a measure of evaporative demand exceeding soil moisture); mean minimum winter (December-February) temperature (TMIN); mean maximum summer (June - August) temperature of summer (TMAX); mean monthly winter precipitation (PPTW); and mean April 1st snow pack (PCK4).
We used LFMM2, which was developed for G2E association and has been shown to outperform other similar approaches with several orders-of-magnitude faster computing , and which also controls for the effects of demographic processes and population structure on the distribution of genetic variation . This approach is robust to high amounts of missing data, such as GBS sequencing tends to produce, when sample sizes are > 100 . LFMM2 regression models combine fixed and latent effects with the following equation
Y = XBT + W + E.
where Y is a matrix of genetic information measured from p genetic markers for n individuals, and X is a matrix of d environmental variables measured for n individuals. The fixed effect sizes are recorded in the B matrix, which has dimension p * d. The E matrix represents residual errors with the same dimensions as the response matrix. The matrix W is a matrix of rank K, defined by K latent factors where K can be determined by model choice procedures. The K factors represent unobserved confounders - usually geographical structure in the genotypes of the samples – represented as an n*K matrix, U. V is a p × K matrix of loadings. The matrix U is obtained from a singular value decomposition (SVD) of the matrix.
W = UVT
To determine K, we used the two approaches implemented in the LEA v.2.6.0 R package: principal component analysis (PCA) and admixture analysis [85, 86]. First, we ran the LEA function pca to select the number of significant PCA components by computing Tracy-Widom tests with the LEA function tracy.widom . Second, we ran the LEA function snmf for K values between 1 and 5 with 10 repetitions each. The most likely K value was identified by minimizing the cross-validation error evaluated in the 10-fold cross-validation procedure (Frichot & Francois, 2014). We then chose significant associations based on a false rate of 5% (q⩽0.05) using the R package QVALUE .
We used fifty half-sib families that span the climatic range of ponderosa pine in California for the greenhouse experiment and G2P analyses. We aimed to have 10 seedlings from each maternal family in both wet and dry treatments, 1000 seedlings in total. During winter 2018, the seeds were stratified to break dormancy by placing them in aerated water for 48 hours, then surface-drying them and placing them in plastic bags in the refrigerator (~ 1.7°C) for 6 weeks. Forty-eight of the 50 families had enough seeds in their cones to be included in the experiment (Additional file 1: Fig. S4).
Because maximum seedling root length in a pilot experiment conducted in 2017 was more than 110 cm, we used plastic tubes with an 8-cm width and 120-cm depth for planting. The bottom of each tube was capped with mesh to prevent the soil from falling out while allowing drainage, and the lightweight clear tubes were wrapped in black plastic to keep roots in the dark. The planting soil was a mixture of 70% sand, 20% vermiculite, and 10% organic-rich potting mix to mimic the coarse texture of the soil of many Sierra Nevada conifer forests . To keep tubes upright we used PVC pipes to build 10 frames that could each hold 100 tubes. Two seeds from each family were planted in each tube in February 2019, and two tubes from each family were randomly placed within each frame. In April 2019, we replanted more stratified seeds of the correct family in any tubes without seedlings. All the tubes were watered every other day during the germination and seedling establishment period (February through June).
At the end of June 2019, all but one seedling per tube was removed, and alternating frames were assigned to the wet treatment and the dry treatment (5 frames containing up to 500 seedlings per treatment) (Additional file 1: Fig. S4). The wet treatment group was watered twice every week and the drought treatment group was watered once every three weeks until mid-October (3.5 months). While wild ponderosa pine seedlings would receive little to no precipitation during the summer months, this occasional watering was necessary in the greenhouse environment to prevent complete mortality. Temperatures inside the greenhouse in the low-elevation environment of Merced, CA reached as high as 37°C on the hottest days and the soil volume of the tubes was limited, with no access to groundwater, both of which make evaporation and drought stress more intense than the no-precipitation condition in the wild.
Multiple phenotypic traits were measured during and after the greenhouse experiment. We calculated shoot growth as final height minus height at the initiation of the treatments. The length of fresh roots was measured from soil surface to taproot tip immediately after the harvesting, to avoid shrinkage. Following harvest, needles, fresh stem and fresh roots of all the seedlings were separately put into paper lunch bags and dried at 75°C for 48 hours. We measured root dry mass (RW) as well as shoot weight (SW, total of stem and needles). We then calculated root-shoot ratio (R2S) as RW/SW. Specific Root length (SRL) was calculated as root length/root weight.
Before harvest, we also collected 3–4 fresh needles from living seedlings to calculate stomatal density. In pines, stomata are arranged into longitudinal rows. We put each needle on a slide and photographed it at 100x magnification using a Leica DME compound microscope equipped with a Leica DFC290 digital camera. All counts were conducted near the middle of the needle to avoid variation that might occur at the base and at the tip. Approximately 1.96 mm lengths of needle were surveyed for number of stomata and stomatal rows on their adaxial (upper) and abaxial (lower) surfaces. Needle width was measured in magnified images using the line measure tool in the Leica software. Then we calculated the stomata density on each side as the number of stomata divided by 1.96*width of needle. Individual seedling means were calculated by averaging adaxial (AD) stomatal density and number of stomatal rows on both sides (AB & AD) across sampled needles. Only 42 out of 48 mother trees had enough germination to carry out these measurements across both treatments.
Genotype-phenotype association analysis
The 42 individual mother trees had already been genotyped, and we used these same SNPs for the G2P association analysis, focusing on the traits significantly associated with drought treatments. The breeding value (BV) of a tree reflects the tendency of an individual to produce offspring with high values of that trait and is estimated by measuring relatives [16, 90]. For the wet treatment traits, we use the average trait value across all members of each family in the wet treatment as the BV. For the drought response traits, we deduct the average trait value for a given family in the wet treatment from the value for each offspring of that family in the drought treatment and then use the mean difference as the BV. We used LFMM 2 (Caye et al. 2019) to run the genotype to phenotype association analysis, and then identified associations based on p (< 10− 5) value.
After identifying the significantly associated SNPs, we aligned the gene sequences for these regions against the nonredundant protein sequences database using UniProt to identify the gene and protein with the implemented Blastx (2.9.0+, E < 1e− 10). The Gene Ontology Annotation Database [91, 92] was used to further identify the potential functions of the genes. If a SNP is located in the intragenic region, we performed a search by querying the flanking sequence 400 bp from the beginning position of the gene. This had to be done separately because the “start” and “end” positions given for the genes containing the introns were too far apart; no hits could be obtained by Blastx.