Occurrences from a study of a changing Lutz spruce (Picea x Lutzii) hybrid zone on the Kenai Peninsula, Alaska

Study Extent

Our study area was the Kenai Peninsula.

Sampling

Initially, we were primarily interested in confirming the taxonomy of the spruce that occurred on the southern Kenai Peninsula prior to deforestation due to spruce bark beetle and wildfire. We acquired two Landsat TM images from USGS Earth Explorer website (http://earthexplorer.usgs.gov/), dated 12 September 1986 (Landsat 5) and 25 September 2014 (Landsat 8). We used ENVI to pre-process the imagery to normalize the data (Radiometric Calibration Tool to calculate reflectance) and remove atmospheric effects (Dark Subtraction). We imported the pre-processed multi-spectral image into ERDAS Imagine. We conducted an unsupervised classification for each time-step and then visually identified the classes that related to spruce forests. Pixels that were spruce forest in 1986 but not in 2014 were classified as deforested. We converted the deforested pixels into a shapefile and dissolved “donut holes” ≤ 0.4 ha. Polygons were buffered and joined to adjacent polygons, excluding the 2014 Funny River Fire. Our focal area for genetics sampling was a 37 790 ha union of major fire polygons south of Tustumena Lake on the southern Kenai Peninsula (1994 Windy Point Fire, 1996 Crooked Creek Fire, 2005 Fox Creek Fire, and 2007 Caribou Hills Fire). Within this area, our sample frame became the centroids of the 250 m pixels from the Alaska eMODIS product, selecting every 12th pixel in both north-to-south and east-to-west axes, making a grid of 58 points spaced at 3 km intervals. This grid ensured we had a representative sample of spruce growing on the southern Kenai Peninsula. We collected spruce needles for genetic analysis from a total of 446 spruce seedlings and adult trees from 56 sites within the sample frame. Two additional sites outside the sample frame were intended to serve as reference populations for parental genomes: white spruce at Silver Lake Trail on the western peninsula and Sitka spruce at Lost Lake Trail on the eastern peninsula. Needles were frozen and stored at −80 ◦C until DNA was extracted.

Method steps

1. DNA was extracted from two spruce needles per individual using the Gentra PureGene Tissue Kit. Prior to extraction, needles were dried for 24 h in a lyophilizer, and ground to a fine powder in a bead beater. We screened 24 microsatellite primer pairs, selecting those that reliably amplified and produced the most scorable bands. Thirteen microsatellite loci, using 12 primer pairs, were amplified from all individuals. One locus proved too difficult to score, so analyses are based on 12 loci from 11 primer pairs (Table 1). All forward primers had an addedM13 tail (CACGACGTTGTAAAAC), which allowed priming of a third, universal M13 primer that was fluorescently labeled, and used for visualization of PCR products. We performed microsatellite polymerase chain reaction (PCR) amplification in 10 μL volume reactions containing 4 ng genomic DNA, 1 U of Kapa3G Plant DNA polymerase, KAPA Plant PCR Buffer, 0.3 μMof each locus-specific primer, and 0.15 μMof the fluorescently labeled M13 primer. PCR amplification conditions were: 1 cycle of 95 ◦C for 3 min; 10 cycles of 98 ◦C for 20 s, a locus-specific annealing temperature (Table 1) for 30 s, and 72 ◦C for 30 s; 28 cycles of 98 ◦C for 20 s, 55 C for 15 s, and 72 ◦C for 30 s; a final extension of 72 ◦C for 15 min. PCR products were sized on an ABI 3730xl genetic analyzer at the Cornell University Institute of Biotechnology. In the genus Picea, mitochondria are maternally inherited, while plastid are paternally inherited. To determine which parental species mitochondrial and plastid genomes were inherited from, we sequenced DNA from one mitochondrial locus (nad5a) and one plastid locus (the trnT-L intergenic space) using the Illumina miSeq. PCR primers were developed by aligning GenBank sequences from all spruce species on the Kenai Peninsula (P. glauca, P. mariana, and P. sitchensis) to ensure that the region amplified encompassed regions where the species differed. We used Primer3plus to identify good priming sites that did not vary between species, and that were less than 500 bp apart, ensuring overlap between the 300 bp forward and reverse miSeq sequence reads. Primers contained iTru tails that allowed the addition of Illumina adapters containing P5 and P7 sequencing primer sites and Adapeterama II indices used to demultiplex individual samples after sequencing. The trnT-L PCR amplifications were carried out in one reaction per sample that included both locus-specific primers, and the indexed iTru primers. However, a single PCR reaction produced too much primer-dimer with nad5a primers, so we did two reactions per sample. The first reaction used the tailed locus-specific primers, and the second used the index primers. Samples were cleaned with Solid Phase Reversible Immobilization (SPRI) beads to remove primers and primer dimer prior to the second PCR reaction. After the second PCR, all samples were pooled, cleaned with SPRI beads, and sent to the University of Alaska Core lab for sequencing on the Illumina miSeq. Locus-specific primer sequences for trnT-L (including iTru tails) were iTru_trnT_3f: 5'-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT AGC TAA GCA GGC TCA ATG GA-3' and iTru_trnT_3r: 5'-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC TTA CTC CCC TTC TCT CGC CAT-3' and primers for nad5a were iTru_nad5_1f: 5'-ACA CTC TTT CCC TAC ACG ACG CTC TTC CGA TCT GAA GGA AGA AGG GGC CCA AG -3' and iTru_nad5_1r: 5'-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC TCG AGC TCT GTT ACC CTT GCA-3'. PCR for trnT-L was carried out in 10 μL volume with 4 ng genomic DNA, 5 μL KAPA HiFi Readymix, and 0.3 μM of each primer (4 primers total). PCR amplification conditions were: 1 cycle of 95 ◦C for 3 min; 25 cycles of 98 ◦C for 20 s, 57 ◦C for 15 s, and 72 ◦C for 15 s; 8 cycles of 98 ◦C for 20 s, 60 C for 15 s, and 72 ◦C for 15 s; a final extension of 72 ◦C for 1 min. The first PCR for nad5a was carried out in 10 μL volume with 4 ng genomic DNA, 5 μL KAPA HiFi Readymix, and 0.3 μM of each primer. PCR amplification conditions were: 1 cycle of 95 ◦C for 3 min; 30 cycles of 98 ◦C for 20 s, 57 ◦C for 15 s, 72 ◦C for 15 s, and a final extension of 72 ◦C for 1min. Cleaned PCR products were used as templates in the second PCR with the same conditions, but only 8 cycles, and an annealing temperature of 60 ◦C.
2. Microsatellite peaks were scored in STRand 2.2.30, and binned in MSatAllele 1.05. Peaks and bins were manually checked and corrected as necessary. The fraction of ancestry from each parental species for each individual was estimated from microsatellite data using the admixture model in the Bayesian clustering software STRUCTURE 2.3.4. Allele frequencies were assumed to be uncorrelated among clusters, and no prior population information was used in clustering. It was assumed that null alleles could be present, but PCR reactions that produced no peaks were assumed to be missing data. Dirichlet parameters (α) for degree of admixture were estimated and different clusters were allowed to have different values (POPALPHAS = 1), allowing for asymmetric admixture. The burnin for the Markov Chain was set to 100 000 steps, and we collected data from 100 000 subsequent steps. To identify the value of K (the number of ancestral populations or clusters) with the highest probability of explaining the data, we compared the mean likelihood of models with different values of K, ranging from K = 1 to 10. Each model was repeated 10 times. The best value of K was chosen as the smallest value of K where the mean log likelihood of the model plateaued, as suggested by the STRUCTUREmanual. We also evaluated K using delta K method. Results were visualized with DISTRUCT and PHYLOGEOVIZ. DNA sequence data from trnT-L and nad5a DNA was analyzed using the SeekDeep pipeline version 2.6.2, intended for analysis of targeted amplicon sequencing, where high-throughput sequencing is used with targeted genes amplified by PCRs. The reads of amplicons shorter than 150 bp were discarded at the filtering step of SeekDeep. Similarly, we discarded reads of amplicons longer than 478 and 550 bp for trnT-L and nad5a, respectively, on the assumption that these were off-target sequences. We used additional arguments “——converge——leaveOutSinglets” for the Qluster step, “——converge——illumina——pop-noErrors—— collapseLowFreqOneOffs” for the processClusters step. The option “——converge” causes the clustering to iterate until there is no more collapsing. These extra options remove the singlet clusters and unreliable and/or low-frequency one-off haplotypes. The extra clean-up options for the processClusters step were recommended for the cases where no PCR replicates are available to help correct for PCR noise. After haplotypes of all individuals were determined, we used fastq_masker of FASTX-Toolkit version 0.0.14 to convert bases with Phred Quality Score less than 20 to ambiguous nucleotides (N). Each sequence read was assigned to a haplotype in SeekDeep, and haplotypes were assigned to White, Black or Sitka spruce based on their match to Gen- Bank sequences. Spatially explicit maps of sampling sites were created using R (version 4.2.2) and packages maps (3.4.1), mapdata (2.3.1), and plotrix (3.8.2).

1. Picea
   rank: genus

The intended extend was the Kenai Peninsula.