Kalmia buxifolia (sand-myrtle, Ericaceae) is disjunctly distributed across  the high-elevation rock outcrops of the southern Appalachians, upper monadnocks and pine savannas of the Carolina Piedmont and Coastal Plain, and the New Jersey Pine Barrens.  Here, we sampled plants from each region and reconstructed the phylogeographic history of  K. buxifolia to test a rock-outcrop Pleistocene refugium hypothesis, estimate the potential direction(s) and timing of migration, and date divergence from its alpine sister species, K. procumbens. We also assess whether isolation in these different environments has led to variation in intrinsic water-use efficiency (IWUE). Dating analysis challenges the current hypothesis that rock-outcrop species are relics of Pleistocene refugia (<18,000 ybp), placing the divergence of K. buxifolia and K. procumbens much earlier, in the late-Miocene (9.40 Ma).  Chloroplast haplotype analysis indicates four potential refugial sites, with the most ancient on Mount LeConte in the Great Smoky Mountains, and point to an Appalachian corridor as the likely Pine Barrens colonization route. The sister species divergence time and population level divergences within K. buxifolia generally coincide with major climatic shifts from the late-Miocene to mid-Pleistocene. Results from CID indicate that plant water-use varies geographically within K. buxifolia, as does leaf morphology, although it is unclear whether this variation is due to genetic adaptation or phenotypic plasticity. These patterns of phylogenetic divergence and resulting ecophysiological diversity within K. buxifolia are significant for clarifying long-held questions about the biogeographic history and trait differentiation within this species. Further, our results suggest that high-elevation rock outcrop communities may have been inhabitated by northern-affinity species for much longer than previously assumed, and that subsequent population disjunction and isolation may have resulted in ecophysioloigcal differentiation in these communities.

Sample Collection – Fifteen study sites (Fig. 1; Table 2) were selected to best represent the distribution of K. buxifolia, and generally followed the sampling strategy of Strand and Wyatt (1991). Two sets of samples were collected at each site at the start and end of the growing season, the first in May/June and the second in late September/October. Both sets of samples were used for carbon isotope analysis, but DNA was extracted only from those collected in May/June. At each site and in each sampling set, ~5cm clippings were taken from 5-10 individuals, placed into 15mL centrifuge tubes, and dried using Silica gel desiccant, -3+8 Mesh Granules (Alfa Aesar, Ward Hill, Masachussets). Ten leaves were also collected from each sampling site for stomatal analysis. Leaves chosen were near the ends of branches, fully exposed, and mature. The underside of each leaf was painted with clear, quick-drying nail polish, and allowed to dry before being placed in a resealable plastic bag. Soil depth and moss mat depth (if present) were also measured in three spots, triangulated around the canopy of the individual plant and then averaged. Individuals were also given a “health” rating from 1-5 based on the percent of the canopy that was browned or appeared dead. Individuals with no visible browning were rated at a 5 while individuals with ≥90% of the canopy browning or dead were rated a 1.

DNA Isolation, Amplification, and Sequencing - Silica-dried leaf tissue from three individuals per population was ground to powder in liquid nitrogen using a mortar and pestle. DNA was then extracted and isolated via the DNeasy plant mini kit (Qiagen, Valencia, California) following the manufacturer’s protocol. Four DNA regions were sequenced, both nuclear and chloroplast. The nuclear ribosomal internal-transcribed spacer (nrITS) was PCR-amplified in its entirety using the external primer pair “ITS4-5” (White et al. 1990). Three chloroplast non-coding regions were PCR-amplified, including the trnL-F intergenic spacer using primers c and d of Taberlet et al. (1991), the trnS-G spacer using primers “trnS (GCU)” and “trnG(UUC)” of Shaw et al. (2007), and the rpS4-trnT intergenic spacer using primers “rpS4R2” and “trnT(UGU)R” of Shaw et al. (2005). These regions were chosen based on previous sequencing done in K. buxifolia and reported high levels of variation in Shaw et al. (2005, 2007). Preliminary samples were sequenced bi-directionally using the same primers as for PCR, but due to the lack of variation found in all markers, remaining samples were sequenced in the single, forward direction using “ITS 4”, “c”, “trnS (GCU)” and “rpS4R2”. Twenty-five microliter PCR reactions for nrITS were prepared using 12.5μL PCR Taq Mixture (Omega Bio-tek Inc., Norcross, Georgia), 0.5 μL forward primer (10 μM), 0.5 μL reverse primer (10 μM), 1μL DNA template, and ddH₂O to fill. PCR reactions for cpDNA regions were prepared with 12.5 μL 2x Platinum^TM SuperFi^TM PCR Master Mix, 1 μL SuperFi^TM GC Enhancer (Invitrogen^TM, Thermo Fischer Scientific Corp., Carlsbad, California), 1.25 μL forward primer (10 μM), 1.25 μL reverse primer (10 μM), 5 μL DNA template, and ddH₂O to fill to 25 μL.

The ITS region was amplified as follows: initial denaturation cyle at 95^oC for 3 minutes; 30 cycles of denaturation at 95 ^oC for 45 seconds, annealing at 52 ^oC for 45 seconds, and extension at 72 ^oC for 2 minutes; followed by a final extension at 72 ^oC for 5minutes. Amplifications for the trnL-trnF intergenic spacer region were performed as follows: initial denaturation at 95 ^oC for 2 minutes; 35 cycles of denaturation at 95 ^oC for 50 seconds, annealing at 50 ^oC for 50 seconds, and extension at 72 ^oC for 1 minute 50 seconds; followed by a final extension at 72 ^oC for 7 minutes. Amplifications for the trnS-trnG intergenic spacer region were performed as follows: initial denaturation at 80 ^oC for 5 minutes; 30 cycles of denaturation at 95 ^oC for 1 minute, annealing at 65 ^oC for 1 minute, and extension at 66 ^oC for 3 minutes; followed by a final extension at 72 ^oC for 2 minutes. The rpS4 gene was also amplified as follows: initial denaturation at 98 ^oC for 30 seconds; 30 cycles of denaturation at 98 ^oC for 10 seconds, annealing at 50 ^oC for 30 seconds, and extension at 72 ^oC for 30 seconds; followed by a final extension at 72 ^oC for 1 minute.

The PCR products were run on a 1% w/v agarose gel with a DNA ladder for sizing, then cleaned and prepared for sequencing using the ExoSAP-IT enzymatic cleaner (Thermo Fisher Scientific, Santa Clara, California). PCR primers were diluted to 5 μM and 15 μL sequencing reactions were prepared using one of the following formulas, depending upon band-strength. Samples with strong bands were prepared with 5 μL, primer, 2 μL PCR product (trnS-G) or 1 μL PCR product (ITS, trnL-F, rpS4-trnT), and ddH₂O to fill. Samples with weak bands were prepared with 5μL primer, 4 μL PCR product (trnS-G) or 2 μL PCR product (ITS, trnL-F, rpS4-trnT), and ddH₂O to fill. Prepared samples were then sequenced by the external vendor GeneWiz (South Plainfield, New Jersey).

Sequence Alignments and Phylogenetic Analyses - Resulting primer sequences for each sample were downloaded into and viewed in Sequencher (GeneCodes Corp, Ann Arbor, Michigan) to compare and confirm sequences, and samples with both forward and reverse primer sequences were combined into a consensus. Outgroup sequences for Phyllodoce nipponica Makino, Phyllodoce aleutica (Spreng.) A. Heller, Kalmiopsis leachiana (Henderson) Rehd, Kalmia latifolia L. and Kalmia procumbens (L.) Gift & Kron were obtained for nrITS, trnS-G, and trnL-F (Phyllodoce sp. unavailable) from GenBank (NIH, Bethesda, MD) (Table 3). Outgroup sequences were then combined with our sequences, aligned in Sequencher, and a Nexus file generated. FindModel (Posada and Crandall 1998) was used to identify the appropriate nucleotide substitution model for each gene region, as determined by the AIC criterion (Akaike 1974). Models were identified for each region as follows: nrITS – GTR + Γ; trnS-G – GTR; trnL-F – F81; rpS4 – F81.  As no variation was found among any of our K. buxifolia nrITS sequences, this region was excluded from further analyses.

Divergence dating for K. buxifolia and K. procumbens was performed using the Bayesian analysis program BEAST v2.4.8 with sequences for the trnL-F and trnS-G gene regions (Bouckaert et al. 2014).  Site model parameters were set and Nexus file converted to a BEAST XML file in BEAUti with partitions for each gene region (Bouckaert et al. 2014). The clock model was set to “Relaxed Clock Log Normal” and the tree prior set to the Birth-Death model with uniform birth and death rates, as this prior is one of those shown to produce the most accurate dating results for datasets containing a mixture of inter- and intraspecific sampling (Ritchie et al. 2016). The Kalmia crown node was calibrated using the fossil Kalmia saxonica, dated to 15.97-125myr (Mai 2001). The calibration prior used a log-normal distribution and the offset was set as a hard-minimum constraint at 15.97 (Ho 2007). The M and S variables were then set as 54.5 and 2.0 respectively to achieve a mean of 70.5 and 95% CI. Markov-chain monte-carlo (MCMC; Yang and Rannala 1997) analyses consisted of 100,000,000 generations, sampling every 10,000 steps. The BEAST output was then assessed in Tracer v1.6 (Rambaut et al. 2013) for chain stationarity and high ESS (over 1,000) and trees were drawn in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/). Clade support is reported as posterior probabilities.

Phylogeographic relationships were analyzed using both TCS (which applies statistical parsimony for network estimation; Clement et al. 2000) and MrBayes (Bayesian inference software) (Huelsenbeck, J. P. and F. Ronquist. 2001). Sequences from genetic markers trnL-F, trnS-G, and rpS4 were used in both analyses with K. procumbens as the outgroup. TCS was run twice, once treating gaps as a fifth character state and once treating gaps as missing data; results were the same for both analyses (see below). For the Bayesian phylogenetic analysis, gaps were treated as missing data, since most gaps were phylogenetically uninformative, or in two cases, located within poly-A strings and ambiguous. Bayesian inference of phylogeny among populations was performed in MrBayes, run remotely on the CIPRES Science Gateway online server (v. 3.1) (Miller et al. 2010). As different site models were identified for these regions (trnS-G – GTR; trnL-F and rpS4 – F81) two runs of MrBayes were performed with a different model applied in each for comparison. MCMC was set to 6,000,000 generations with two simultaneous runs and four chains per run. The percentage of trees discarded as burn-in was 0.25. To assess stationarity of likelihood at the end of each run, the value of potential scale reduction factor (PSRF) was confirmed to be approaching 1.0 with a standard deviation of split frequencies approaching 0.0. The ESS statistic was confirmed to be above 100, indicating that the tree sampling from the likelihood distribution was adequate, and the overlay plots of generation vs. log likelihood values showed stationarity for both runs.

Carbon Isotope Analysis - Five individuals were selected from each of the 12 sampling localities in the Carolinas (excluding the two NJ populations). Seven of these populations represented mountain sites, while five were from the Piedmont/Coastal Plain. The pair of samples collected at the start and end of the growing season from each individual were prepped for carbon-isotope analysis. Samples were oven-dried at 55^oC for at least 48hrs and up to 96hrs, then ground to a homogeneous powder using a Mini-Beadbeater-8 (MBB-8) (Biospec Products Inc., Bartlesville, Oklahoma). Leaf-tissue was loaded into 2ml screw-cap microtubes with o-ring seals, along with three glass beads (2.5mm). The MBB-8 was then set to ‘Homogenize’ and run for 2min.

After samples were ground to a powder they were encapsulated and weighed. Tin capsules (4.75x 11mm) were weighed on a micro-balance accurate to 0.001 mg, tared, and then loaded with sample to 0.7 – 1.5mg. Capsules were then sealed and crimped using forceps down to a spherical package (<2 mm), and dropped to confirm there was no leakage. All samples were loaded into a 96-well microtiter plate and sent to the Stable Isotope Analysis Lab in the College of Marine Sciences at the University of South Florida (St. Petersburg, Florida) and analyzed for dC13 (δ¹³C) in a continuous flow isotope ratio mass spectrometer. The carbon isotope composition (δ) was calculated relative to the international Pee Dee Belemnite (PDB) standard (Farquhar et al. 1989) as δ_plant = (R_sa – R_sd)/R_sd X 1000, where R_sa is the ¹³C/¹²C ratio of the sample and R_sd is that of the standard. Carbon isotope discrimination (∆) was then estimated as ∆ = (δ_air – δ_plant)/(1+ δ_plant/1000), where -8.0‰ = δ_air (¹³C composition of atmospheric CO₂) (Farquhar et al. 1989; Adriegjo et al. 2014). IWUE was then interpreted from ∆, as generally higher discrimination (resulting from higher internal CO₂) represents higher IWUE.

A principal components analysis (PCA) was performed to examine the association of ∆ with environmental variables (soil depth, moss mat, and “health” rating). Samples were identified by geographical region as either “Mountains” (seven sampled populations) or “Piedmont/Coastal Plain” (five sampled populations) (Table 2). Analysis was performed in PC-ORD^TM 7 (McCune and Mefford 2016). Additionally, a parametric two-way ANOVA with Bonferroni correction was used to test the effects of region, season, and region x season on ∆ (α=0.05).

Stomatal Analysis - Stomatal density was assessed as an ancillary method for determining variation in transpiration across regions and compared to results from ∆. Leaf peels were generated from the samples collected from each of the 12 Carolina localities. In the lab, the nail polish layer on the underside of the leaf was over-laid with clear double-sided tape, then carefully peeled off and mounted on a slide with an imprint of the lower epidermis. Using a compound microscope at 400x, the total number of stomata in the view field (0.45mm area) were counted in three different areas and then averaged for the approximate stomatal density of the leaf. The radial length and width of each leaf imprint was then measured from the center to compare density with leaf area.  Area of the leaf was then calculated as an ellipse using the equation A=πab,  where area (A) is the product of π , the longest radial length (a), and the longest radial width (b). Two-sample t-tests were used to test for significance between the mean stomatal density and mean leaf area of the two regions (α = 0.05). Finally, a linear regression analysis was performed to test for the influence of stomatal density on CID (∆).

The WUEData.xlsx file contains all data used for the PCAs, seasons separated by fill color (white = spring, grey = fall sampling). CID was converted to delta as (outlined in the methods) for all other analyses and these values can be found in the LeafCharacters+CID. xlsx, along with measures of stomatal density and leaf area. The nexus file contains all aligned sequences for all three gene regions used in phylogeographic analyses.