Climatic oscillations throughout the last ice age played a significant role in shaping many of the plant and animal distributions we see today. The expansion-contraction (EC) model provides a template for understanding the impact of climate cycling on post-Pleistocene distributions of flora and fauna (Taberlet et al. 1998; Hewitt 1999; Provan and Bennett 2008). During glaciation, ranges contract as organisms track suitable habitat towards refugia in ice-free latitudes. Species that do not undergo range shifts fast enough can go extinct, and populations that do establish at ice-free latitudes are susceptible to genetic bottlenecks and founder effects (Taberlet et al. 1998; Hewitt 1999; Provan and Bennett 2008). As the climate warms and ice sheets retreat, species then expand back toward the poles. This cycle has been repeated many times throughout the Quaternary and can become complex when ecological factors such as niche specialization are also considered in predicting or explaining species' establishment patterns (Hewitt 1996, 2000; Taberlet et al. 1998; Ikeda and Setoguchi 2006; Provan and Bennett 2008). The EC model has long been used to explain the isolation of northern species in high-elevation communities of the southern Appalachians (latitude 35°–37°) (Ramseur 1960; Wiser 1994; Godt et al. 1996; Wiser et al. 1996). During post-glacial expansion periods, populations are thought to have become isolated (or ‘left behind') on high mountain slopes as lower elevation populations disappeared because they could not acclimate or were pushed out by expanding hardwood forests. This has resulted in a distinctive flora for high-elevation communities, comprised of rare and disjunct species such as Sibbaldiopsis tridentata (Aiton) Rydb. and Minuartia groenlandica (Retz.) Ostenf. (Wiser 1994).

One characteristic species of such high-elevation habitat is Kalmia buxifolia (Bergius) Gift & Kron (sand-myrtle, Ericaceae), a long-lived, woody, flowering shrub that is disjunct across high-elevation rock outcrops of western North Carolina, the upper monadnocks of the North Carolina Piedmont, pine savannas of the Carolina Coastal Plain, and the New Jersey Pine Barrens. Due to its distribution and variability in morphological traits such as leaf size and shape, circumscription of sand-myrtle has long been debated (Table 1; Small 1903; Camp 1938; Wood 1961; Wilbur and Racine 1971). In their investigation, Strand and Wyatt (1991) found significant among-region variation in seven morphological characters of K. buxifolia (= Leiophyllum buxifolium (Bergius) Elliott), but allozyme variation was lower than that found in other long-lived woody species, and no clear pattern of genetic divergence could be determined. In 1996, Kron and King moved sand-myrtle into the genus Kalmia based on phylogenetic analyses of the chloroplast rbcL gene and the nuclear ribosomal internal transcribed spacers (ITS) region, which showed that both Leiophyllum and the monotypic Loiseleuria were included in a Kalmia clade. Within Kalmia, K. buxifolia and its sister species K. procumbens (= Loisleuria procumbens (L.) Desvaux; alpine azalea) are derived and united in having petals not fused the length of the corolla and lacking the anther pockets and tension-loaded pollen dispersal mechanism that characterizes the rest of the genus (Gillespie and Kron 2010, 2013). Kalmia procumbens, an alpine-circumarctic species found in the high peaks of New England, USA, and northern boreal regions around the world (Ikeda et al. 2017), is the only member of Kalmia not restricted to North America and Cuba. The most recent common ancestor (MRCA) of K. buxifolia and K. procumbens is thought to have occurred at southern latitudes in North America during the late Pliocene/early Pleistocene, with divergence in the expansion period that followed (Ikeda et al. 2017).

Table 1.

Twentieth century taxonomic treatments of Kalmia buxifolia by distribution (adapted from Strand and Wyatt 1991).

Given the current disjunct distribution, several methods for colonization of the Pine Barrens by K. buxifolia following divergence could be hypothesized. First, although Ikeda et al. (2017) suggest a more recent, southern divergence, the Pine Barrens populations may be relicts of the ancient arcto-tertiary flora that were never actually connected with southern populations and persisted through the last glacial maximum (LGM) in their current location. Although unlikely, palynological, paleoclimate, and geological evidence all suggest that the Pine Barrens may have remained ice-free during the Pleistocene and served as an important refugium for northern and southern affinity flora (Potzger 1945; French et al. 2003; Stanford 2015). Alternatively, two possible corridors (coastal and mountain) may have connected southern populations to their northern habitat. Climatic oscillations throughout the Pleistocene are thought to have resulted in the rise and fall of sea level at least every 100,000 yr, which would have opened a pine-forest corridor along the eastern continental shelf, connecting the Carolina coast to New Jersey for several thousand years at time, facilitating species exchange (Pisias and Moore 1981). Another possibility is that the Appalachian Mountains instead served as a cooler, rocky corridor from the Carolinas to the Pine Barrens. This corridor is supported by historical documentation of K. buxifolia's presence in eastern Kentucky (Whitely Co.; observations and herbarium specimens) and eastern Pennsylvania (Monroe Co. and Northampton Co.; observations only, now extirpated), suggesting a historically more robust and interconnected Appalachian population that has since been extirpated (Pennsylvania Natural Diversity Inventory; Medley 1993; Rhoads and Klein 1993; Thompson 1997; Campbell and Medley 2012; Kartesz 2014). As the historical Pennsylvania and modern New Jersey populations are separated by fewer than 100 miles, species exchange may have been supported among them had there once been suitable habitat. Jump dispersal, however, is unlikely as the fruits are small, dehisce before abscission, and lack elaiosomes or any other unusual features, suggesting gravity as the most likely dispersal mechanism (Strand and Wyatt 1991). One key to deciphering which colonization route was utilized is to locate refugial populations. Applying the EC model, refugial populations would be expected in the southernmost part of a distribution and they should harbor the greatest genetic diversity, while recolonized populations should retain only part of the genetic structure from their refugial ancestors (Ikeda and Setoguchi 2006).

Chloroplast DNA (cpDNA) can be particularly useful in locating refugia and tracing post-glacial colonization in plants, as cpDNA is typically maternally inherited, which greatly reduces the effective population size, and thus retains finer geographical structure (Vendramin et al. 1999; Petit et al. 2002; Ikeda and Setoguchi 2006). Additionally, cpDNA's non-recombining characteristics can preserve haplotypes over generations with relatively few changes among sites (King and Ferris 1998; Petit et al. 2002; Ikeda and Setoguchi 2006). Non-coding and spacer regions in both nuclear and cpDNA are particularly helpful in tracing migrations and assessing relationships within a genus or among populations due to their high mutation rates (White et al. 1990). Gillespie and Kron (2013) found a combination of nuclear and cpDNA gene regions, both coding and non-coding (rbcL, matK, ndhF, trnS-G-G, nrITS, and WAXY), sufficient for resolving relationships within the tribe Phyllodoceae of Ericaceae. Phylogeographic studies of K. procumbens and other members of Phyllodoceae have also found the cpDNA trnL spacer to be highly variable (Taberlet et al. 1991; Ikeda and Setoguchi 2006; Ikeda et al. 2017). The application of such molecular tools to resolve the phylogeographic relationships within K. buxifolia and its divergence from K. procumbens may be key in predicting how the species will respond to variation in current and future climatic pressures.

In the southern Appalachians, the high-elevation rock outcrop habitat of K. buxifolia has unique environmental characteristics and community make-up. Outcrops are fairly uncommon on the landscape and are isolated from each other, which contributes to the high occurrence of rare and endemic species in their community make-up (Wiser et al. 1996). Greater than 17% of southern Appalachian outcrop flora are northern-affinity endemics or disjunct taxa, hypothesized to be remnants of glacial tundra or alpine vegetation, whose ranges have retracted in post-Pleistocene warming (Ramseur 1960; Wiser 1994, 1998; Godt et al. 1996; Wiser et al. 1996). Typical rock outcrops have shallow, acidic soils, and limited water retention and availability, which results in water stress and limited plant growth (Culatta and Horton 2014). Frequent cloud immersion (fog) may be essential to ameliorating this stress by decreasing the leaf-to-air vapor pressure deficit (VPD); a lower VPD reduces water loss through transpiration and increases plant water-use efficiency (WUE) (Culatta and Horton 2014). High-elevation forests (> 1500 m) of the southern Appalachians have been found to experience significant cloud immersion on 60–75% of days during the growing season, with that fog contributing up to 31% of the overall water budget in high-elevation forests (Berry and Smith 2012; Berry et al. 2013; Culatta and Horton 2014). Direct foliar uptake of fog has been documented in other rock-outcrop species that experience daily cloud immersion and may be an important source of moisture for K. buxifolia (Limm and Dawson 2010; Berry and Smith 2014; Culatta and Horton 2014). Rock outcrops are also high-light environments, and cloud immersion may be important for reducing insolation, which can limit plant productivity when soil moisture is low (Horton and Culatta 2016).

Although geographically separated, the coastal sites on which K. buxifolia is found share many commonalities with their mountain counterparts. These commonalities include poor soil water retention and nutrient poor, acidic soils. As Coastal Plain soils drain water quickly, sites frequently alternate between saturation and prolonged drought, which intensifies in summer months and is likely to become more frequent with climate change. Although all populations of K. buxifolia are somewhat drought-tolerant as indicated by the soils on which they live, it is reasonable to predict that Coastal populations may use water differently from their mountain counterparts, which are adapted to daily cloud immersion and possible foliar uptake. Intrinsic water-use efficiency (IWUE) relates net photosynthesis to stomatal conductance (A/g) and can be useful in assessing longterm water use and drought stress in individuals. As water becomes limited and stress rises, individuals often close their stomata to reduce water loss, which reduces CO₂ uptake and eventually internal CO₂ and net photosynthesis. Intrinsic WUE is frequently inferred from carbon-isotope discrimination (CID, Δ) in the leaves, or the ratio of δ¹³C to δ¹²C, as Δ is interpreted as inversely correlated with IWUE (Marshall et al. 2007). Over months or seasons, greater Δ, or discrimination against δ¹³C, can indicate higher internal CO₂ concentrations associated with greater stomatal conductivity, and thus lower IWUE. Variation in IWUE may also be apparent in functional traits such as leaf stomata, where plants with ample water will often have high stomatal density that allows them to capitalize on photosynthesis. In comparison, water-stressed plants tend to have a lower stomatal density to reduce water loss (Milburn and Weatherley 1971), although agreement between the observed spatial patterns of morphological characters (i.e. stomatal density) and physiological characters (i.e. IWUE) is not necessarily expected as they may be under different selection pressures (Ferguson 1980; Strand and Wyatt 1991).

Here, we combine the tools of phylogeography with those of ecophysiology to resolve long-held questions of the biogeographic history and diversity within K. buxifolia. The main objectives of this study were to: 1) date divergence from the sister species K. procumbens and reconstruct the phylogeographic history of K. buxifolia to determine the direction(s) of post-refugial migration; and 2) assess whether isolation in environments with differing drought-stress regimes (high-elevation rock outcrops and Coastal Plain) has led to variation in carbon isotope discrimination (Δ), and by inference, intrinsic water-use efficiency (IWUE) in disjunct K. buxifolia populations. We tested the EC hypothesis that range contraction and expansion during and after the Pleistocene has left scattered remnant populations on high-elevation rock outcrops, isolated monadnocks across the Piedmont, along the Suffolk escarpment of the Carolina Coastal Plain, and at Coastal sites. Additionally, looking to the future, we asked if current differences between mountain and Coastal Plain plants in traits related to IWUE (e.g. stomatal density) indicate the potential for these plants to adapt to climate warming effects such as a longer growing season or harsher summer drought. Such adaptative potential has already been documented in other characteristic rock outcrop species such as Hydatica petiolaris (Raf.) Small (=Micranthes petiolaris (Raf.) Brouillet & Gronall, cliff saxifrage) and Solidago simulans Fern. (granitic dome goldenrod), which are able to acclimate to decreasing cloud immersion through decreased transpiration and increased WUE (Horton and Culatta 2016). Understanding this potential in K. buxifolia is particularly important if the mountain populations are in fact a glacial refugium, and thus may retain the genetic diversity to select for coast-like drought tolerances as the climate changes. Further, understanding the divergence patterns of K. buxifolia and resulting ecophysiological variation may help to characterize other rock-outcrop species with similar mountain-coastal disjunct distributions (i.e. Zigadenus leimanthoides A. Gray (pine barren deathcamas)) and inform conservation and management planning for these unique systems.

Materials and Methods

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, ∼5 cm clippings were taken from 5–10 individuals, placed into 15 mL 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 a 5 while individuals with ≥ 90% of the canopy browning or dead were rated a 1. All data related to this investigation have been deposited in the Dryad Digital Repository (Quinlan et al. 2020) and GenBank accession numbers are listed in Appendix 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 intron using primers c and d of Taberlet et al. (1991), the trnS-G intergenic 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^™ SuperFi^™ PCR Master Mix, 1 µL SuperFi^™ GC Enhancer (Invitrogen^™, 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°C for 3 minutes; 30 cycles of denaturation at 95°C for 45 sec, annealing at 52°C for 45 sec, and extension at 72°C for 2 minutes; followed by a final extension at 72°C for 5 minutes. Amplifications for the trnL intron region were performed as follows: initial denaturation at 95°C for 2 minutes; 35 cycles of denaturation at 95°C for 50 sec, annealing at 50°C for 50 sec, and extension at 72°C for 1 minute 50 sec; followed by a final extension at 72°C for 7 minutes. Amplifications for the trnS-trnG intergenic spacer region were performed as follows: initial denaturation at 80°C for 5 minutes; 30 cycles of denaturation at 95°C for 1 minute, annealing at 65°C for 1 minute, and extension at 66°C for 3 minutes; followed by a final extension at 72°C for 2 minutes. The rpS4 gene was also amplified as follows: initial denaturation at 98°C for 30 sec; 30 cycles of denaturation at 98°C for 10 sec, annealing at 50°C for 30 sec, and extension at 72°C for 30 sec; followed by a final extension at 72°C for 1 minute.

Fig. 1.

Kalmia buxifolia sampling localities. Names given as 4-digit ID (Table 2), colors and symbols identify the region.

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, 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, 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 (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. The FindModel web implementation ( http://hiv.lanl.gov/content/sequence/findmodel/findmodel.html) of Posada and Crandall's (1998) Modeltest script was usedt o 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 – 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 v. 2.4.8 with sequences for the trnL 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–125 myr (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 1 × 10⁸ generations, sampling every 1 × 10⁴ steps. The BEAST output was then assessed in Tracer v. 1.6 (Rambaut et al. 2013) for chain stationarity and high ESS (over 1000) and trees were drawn in FigTree v. 1.4.4 ( http://tree.bio.ed.ac.uk/). Clade support is reported as posterior probabilities.

Table 2.

Kalmia buxifolia sampling localities and data collected from each (samples per population = n).

Phylogeographic relationships were analyzed using both TCS (which applies statistical parsimony for network estimation; Clement et al. 2000) and MrBayes (Bayesian inference software; Huelsenbeck and Ronquist 2001). Sequences from genetic markers trnL, 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 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 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.

Table 3.

GenBank accessions for outgroup taxa.

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°C for at least 48 hr and up to 96 hr, then ground to a homogeneous powder using a Mini-Beadbeater-8 (MBB-8) (Bio-spec Products Inc., Bartlesville, Oklahoma). Leaf tissue was loaded into 2 ml screw-cap microtubes with o-ring seals, along with three glass beads (2.5 mm). The MBB-8 was then set to ‘Homogenize’ and run for 2 min.

After samples were ground to a powder they were encapsulated and weighed. Tin capsules (4.75 × 11 mm) were weighed on a micro-balance accurate to 0.001 mg, tared, and then loaded with sample to 0.7–1.5 mg. 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 × 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). 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 “Mountain” (seven sampled populations) or “Piedmont/Coastal Plain” (five sampled populations) (Table 2). Analysis was performed in PC-ORD^™ 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 × 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.45 mm 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 (Δ).

Results

Divergence Dating Analysis—Divergence analysis in BEAST with the crown node of Kalmia calibrated with Kalmia saxonica (Mai 2001) dated the divergence of K. buxifolia and K. procumbens to 9.40 Ma (95% HPD ± 2.35–17.35; Fig. 2; Table 4). The wide HPD interval may have been the result of the uncertainty around the calibrating fossil date (15.97–125 Ma), which was treated conservatively during analysis by placing a hard-minimum constraint on the crown calibration (15.97 Ma) but a relatively soft maximum. Another possible cause is topological uncertainty in the tree due to the low number of informative sites and sequences. However, Zhu et al. (2015) showed that while divergence time estimation can be improved by adding more data, this does not eliminate uncertainty in the time estimates. Regardless of the source of the uncertainty, our overall conclusions are not based on precise divergence times, but rather on exclusion of species divergence during the Quartenary period. This analysis places the divergence of the two species in the late Miocene (∼5–24 Ma), with the tails of the estimate falling very early Pleistocene to early Miocene. Further population divergence in K. buxifolia occurred throughout the Pliocene-Pleistocene, with major splits dated to 3.89 Ma (± 0.45–9.50) and 2.22 Ma (± 0.20–5.93) and more recent divergence events occurring from 1.63–0.24 Ma (± 0.001–4.36).

Fig. 2.

BEAST analysis tree, dating the divergence of Kalmia buxifolia and K. procumbens. Node labels show mean estimated crown ages (bold), branch values show posterior probabilities, colors/symbols at tips indicate geographic region. Divergence of these sister species is dated to 9.40 Ma (95% HPD = 5.08–26.84 Ma).

Phylogenetic Analysis—Analysis of haplotype diversity and distribution in TCS resulted in a network of 10 haplotypes, with haplotype G containing the majority of samples (Fig. 3). Haplotype J (followed by H and I) is the most ancestral type relative to the outgroup K. procumbens (not shown in Fig. 3), which is connected to the H haplotype by 20 mutations. The Mount LeConte population (in GSMNP) contained the greatest haplotype diversity (three haplotypes) as well as the most ancestral type (J), suggesting it is a refugial population. The three other Mountain populations that contain more than one haplotype (Satulah, Laurel Knob, and Little Green Mountain) are all related to at least one of the types found on Mount LeConte, with two of the three (H, found on Laurel Knob; I, found on Satulah) likely derivatives of the most ancestral type (J). All but one of the coastal populations (Peachtree, Gordon Butler, and Orton Creek) also had two haplotypes, with Gordon Butler identical to those found on Laurel Knob and Orton Creek's type E one step away from type F found at Peachtree. Both New Jersey populations (Oswego River and Wharton Creek) exhibited the common haplotype (G), which was found within all regions and all but one population (Mount LeConte). The Oswego River population also contained a second haplotype (A), which is the most derived haplotype in the network and likely a derivative of the mountain haplotype D.

Table 4.

Summary of divergence time estimations for crown nodes in BEAST. Node (species name if outgroup, locality if population-level divergence within K. buxifolia), mean node age, and 95% HPDs from BEAST analyses.

Phylogenetic analysis in MrBayes resulted in nearly identical consensus trees regardless of model implemented (F81 or GTR) and thus trees were combined in Fig. 4. Both trees were largely inconclusive with nearly all samples clustering in a polytomy. However, an individual from the Mount LeConte (MTLC1) was unresolved relative to the outgroup, K. procumbens, rather than grouping with the rest of K. buxifolia in both runs (Fig. 4). This may support the findings of the haplotype network (Fig. 3), where the Mount LeConte haplotype (J) was the most ancestral compared to K. procumbens. Also consistent in both trees, the remaining Mount LeConte samples (MTLC5, MTLC10) group as sister (pp = 0.96, 0.97; Fig. 4), which may reflect the uniqueness of the haplotypes in the population found in TCS (Fig. 3). The grouping of GBNP1 and LAKN1 as sister to each other (pp = 0.77, 0.79; Fig. 4) and the larger polytomy also supports the findings of the TCS analysis in which these populations contained the same two haplotypes (G, H), where H is likely a derivative of the most ancestral type (J) found on Mount LeConte.

Fig. 3.

Kalmia buxifolia haplotype map. Haplotypes generated from sequences of three cpDNA markers (trnL intron, trnS-G intergenic spacer, and rpS4 intergenic spacer). A. Haplotype network, with each network line segment representing one mutation. B. Pie charts over sampling localities show the mix of haplotypes found in that population.

Stable Isotope and Stomatal Analysis—A PCA was conducted to examine the relationships between variation among populations in Δ and environmental variables (Fig. 5). PC axes 1 and 2 accounted for more than 65% of sampling variance (Spring: PC1 = 41.497%, PC2 = 23.592, PC3 = 23.031%, PC4 = 11.88%; Fall: PC1 = 41.767%, PC2 = 23.842%, PC3 = 21.200%, PC4 = 13.192%), and held eigenvalues close to 1 in both sampling sets (Fig. 2; Spring: PC1 = 1.660, PC2 = 0.944, PC3 = 0.921, PC4 = 0.475; Fall: PC1 = 1.671, PC2 = 0.954, PC3 = 0.848, PC4 = 0.528). Δ and average soil depth were strongly and positively correlated with PC1, suggesting sites to the right of the axis (i.e. Piedmont and Coastal populations) are characterized by lower IWUE (as it is negatively correlated with Δ) and deeper soils (Tables 5, 6). Similarly, average moss mat depth and average health score were strongly correlated with PC2, although the relationship was positively correlated with the Spring sampling and negatively correlated with the Fall due to the differently rotated axes. Regardless, these relationships followed a similar pattern in both samplings. Samples clustered by region along both axes, with Mountain samples exhibiting greater variation (particularly along PC2) than their Piedmont/Coastal counterparts (Fig. 5). In both PCAs, some samples from the Mountain set and Piedmont/Coastal set overlapped around 0 along PC1. However, this relationship was not indicative of any population but rather individuals from several populations, demonstrating the variation within both sites and regions rather than populations as a whole.

Further statistical analysis of carbon isotope discrimination showed there was a significant effect of region (Mountain vs. Piedmont/Coastal Plain) on Δ (F = 87.55, p < 0.001), with Piedmont/Coastal sites averaging higher Δ in both sub groups (Mountain: Spring mean Δ = 30.59 ± 1.28, Fall mean Δ = 31.04 ± 1.21; Piedmont/Coast: Spring mean Δ = 32.71 ± 1.26, Fall mean Δ = 33.10 ± 0.99). However, there was no significant effect of season (Spring vs. Fall) (F = 3.56, p = 0.062), or interaction between region and season (F = 0.01, p = 0.906). As IWUE is inversely proportional to Δ, we therefore infer that IWUE did not vary in K. buxifolia within a region over the growing season, but was consistently lower in the Piedmont/Coastal plants (Fig. 6).

Fig. 4.

MrBayes consensus tree for Kalmia buxifolia. Tree is based on three cpDNA markers (trnL intron, trnS-G intergenic spacer, and rpS4 intergenic spacer) with K. procumbens as the outgroup. Posterior probabilities displayed on nodes show results from F81 and GTR models, respectively. Symbols indicate geographical region.

Mean stomatal density was significantly greater (p < 0.00001) and mean leaf area significantly smaller (p < 0.006) in Mountain versus Piedmont/Coastal Plain populations (Fig. 7). Regression analysis between Δ and stomatal density indicated a significant (p < 0.001), but weak relationship with r² = 0.03 (Fig. 8).

Discussion

Dating analysis of the divergence of K. buxifolia and K. procumbens calls into question the hypothesis that southern Appalachian rock-outcrop species are of Pleistocene origin (maximum 2.58 Ma). Here, the divergence of these sister species is placed much earlier, in the late-Miocene (9.4 Ma ± 2.35–17.35), although we cannot definitively rule out an early Pleistocene divergence because of the broad confidence interval. This study suggests that alpine relatives inhabited these rock outcrops since at least early Pliocene (3.89 Ma ± 0.45–9.5), the inferred age of the oldest outcrop population on Mount LeConte. Such a Miocene origin of K. buxifolia is supported by dating of another common rock-outcrop species, Micranthes petiolaris (Raf.) Brouillet & Gornall, whose divergence was placed at ca. 14 Ma (K. Mathews unpubl. data). Following paleoclimatic evidence, it seems likely that a late-Miocene divergence for these species was largely brought on by climatic shifts. The period following the Terminal Eocene Event (35 Ma) until the late-Miocene was characterized by frequent oscillations between warm/humid and cool/dry climates that would have facilitated extensive range expansion and contractions (Zachos et al. 2008; Meseguer et al. 2015). The Middle-Miocene Climate Transition (∼14 Ma), which was a sudden decrease in average global temperature, is thought to have triggered diversification events across the northern hemisphere (Herbert et al. 2016). Both Wen et al. (2010) and Manos and Meireles (2015) found that the majority of divergences dated for North American – East Asian disjunct clades fall within the Miocene and attributed the splitting largely to climatic oscillations. Manos and Meireles (2015) in particular investigated the origins of over 250 southern Appalachian woody species representing 158 clades and found divergence times for over 60% of the clades placed in the last 20 Myr.

Fig. 5.

Principal components analysis (PCA) of Spring (top) and Fall (bottom) Kalmia buxifolia samplings. PCA plots PC1 vs. PC2 with eigenvectors.

Table 5.

Eigenvectors scaled to standard deviation, Kalmia buxifolia Spring sampling.

Table 6.

Eigenvectors scaled to standard deviation, Kalmia buxifolia Fall sampling.

Population-level divergence for K. buxifolia was dated to mid- to late-Pliocene (3.89 Ma ± 0.45–9.50) and more recent divergence events through the Pleistocene (2.22–0.24 Ma ± 0.001–5.94), a time span also strongly characterized by fluctuations in climate and geographic changes. A major warming period (the Mid-Pliocene Warming Event) occurred between 3–4 Ma, coinciding with the major splitting event between the Mount LeConte populations and all other K. buxifolia populations observed in our analysis. This warming not only resulted in habitat shifts but also a significant rise in sea level, which is thought to have peaked 30 m higher than it is today and likely contributed to extinction and isolation of some of the more coastal populations (Dwyer and Chandler 2009). Late Pliocene – early Pleistocene diversification has been documented in Sabatia, another plant group with a similar Appalachian – Coastal Plain disjunction (Mathews et al. 2015). The end of the Pliocene and beginning of the Pleistocene marked the transition from a warm, tropical climate to a much colder global climate than seen today in the northern hemisphere. Another major population divergence in our analysis, between a southern Blue Ridge escarpment population (Laurel Knob) and all others, coincides with this transition, occurring at the start of the Pleistocene (2.22 Ma ± 0.20–5.93), and was yet again likely in response to changes in the global climate and subsequent habitat range shifts.

Haplotype network analysis further explored the phylogeographic patterns at the population level. The Mount LeConte population in GSMNP not only resolved as the oldest population (3.89 Ma ± 0.45–9.50) of those assessed in this study, but also contained the greatest haplotype diversity (including the most ancestral and most derived types). Thus, in this analysis, Mount LeConte appears to be the oldest remaining refugium and the most historically isolated population. Interestingly, Mount LeConte also is the only population sampled from the western side of the Blue Ridge Divide, which may have contributed to its isolation and haplotype divergence. Meanwhile, although the populations on the eastern side of the mountains are dispersed over a wide geographic area, there are few barriers to their dispersal and it is plausible that they may have formed more interconnected populations regularly over the last several hundred millennia.

Fig. 6.

Boxplots of Δ (CID) by region and sampling season (Spring or Fall). Plots show median, 10th, 25th, 75th, and 90th percentiles as vertical boxes with error bars. Statistical analysis by parametric two-way ANOVA showed a significant effect of region (F = 87.55, p < 0.001) but not season (F = 3.56, p = 0.062) or interaction between region and season (F = 0.01, p = 0.906).

Fig. 7.

Mean differences in stomatal density (a) and leaf area (b) in Kalmia buxifolia by region, with stomatal density being higher in Mountain plants but leaf area higher in Piedmont/Coastal plants. Differences between regions were significant by both metrics (p = 1.84e^–5, t = 4.34, df = 379; p = 0.0058, t = 2.78, df = 295).

Fig. 8.

Regression of carbon isotope discrimination (Δ) and stomatal density in Kalmia buxifolia. Relationship is significant (p = 0.0008) but weak (r² = 0.03). Intrinsic water-use efficiency (IWUE) is inferred as the inverse relationship of Δ.

Evidence for interconnected populations can also be seen in the haplotype analysis, as all populations (excluding Mount LeConte) contained the same “G” haplotype, which is neither very ancestral nor recently derived. Additionally, a Mountain population (Laurel Knob) and a Coastal Plain population (Gordon Butler Nature Preserve) share a unique, ancestral haplotype (“H”), suggesting that at least those two populations were once interconnected. Although Mount LeConte appears to be the most ancestral population sampled, there is also support for at least two more recent Pleistocene refugia, in the mountains and along the Suffolk escarpment. Satulah Mountain, Laurel Knob, and Little Green Mountain all exhibited some haplotype diversity, with the common “G” haplotype and a second unique type. Of these, the unique types found on Satulah (“I”) and Laurel Knob (“H”) are both derivatives of the most ancestral type found on Mount LeConte (“J”), while Little Green (“D”) is a precursor to the derived types on Mount LeConte (“B” and “C”). Thus, these mountain sites (particularly Satulah and Laurel Knob) may be representative of a second refugium retaining these older haplotypes or possibly evidence for greater gene flow between the western and eastern sides of the mountains than expected.

The third potential refugium lies along the Suffolk escarpment, with Peachtree Rock Heritage Preserve (“F”) and Gordon Butler Nature Preserve (“H”) also exhibiting a second haplotype. Coastal flooding over the Pliocene and Pleistocene would have extirpated any populations that established east of this line, but the Suffolk escarpment itself was never underwater. Therefore, these populations may be remnants of a once more robust coastal population. The most coastal population, Orton Creek, also exhibits a unique haplotype (“E”), however this type is a derivative of the Peachtree Rock type and likely emerged after sea level subsided and the coast was recolonized from the remaining escarpment populations.

The two New Jersey populations sampled (Wharton State Forest and Oswego River Preserve) contained the common “G” haplotype as well as second type (“A”) found only in the Oswego River population. This “A” haplotype is one of the most derived in the network and a derivative of the “D” type found in the Little Green Mountain population. This supports the hypothesis that the Pine Barrens populations were colonized via an Appalachian mountain corridor, potentially by dispersal from the historical populations in Pennsylvania. Additionally, the Pine Barrens are known to never have flooded (being structurally high like the Suffolk escarpment; Stanford 2015) and although periglacial, were never under ice. It is possible that these sites were colonized via a mountain corridor during the early-Pleistocene, became isolated during the height of Pleistocene glaciation, and diverged.

Characterizing the directionality of population divergence and isolation in K. buxifolia has important implications for understanding the ecophysiological divergence between plants in the two geographical regions, demonstrated by carbon isotope discrimination analyses (Δ). Counter to our hypothesis that Mountain populations would have lower IWUE (higher Δ) due to frequent cloud immersion and lower drought stress, Piedmont and Coastal populations had lower IWUE over the 2017 season. This finding is similar to results from experiments in Arabadopsis thaliana, which found selection for lower WUE was stronger under drought conditions compared to well-watered conditions (Kenney et al. 2014). Our analyses also revealed regional differences in leaf morphology, as the Mountain plants had, on average, higher stomatal density but lower leaf area in comparison to their Coastal relatives. This follows the trend predicted by Milburn and Weatherley (1971): plants with ample water would have higher stomatal density, which allows them to capitalize on photosysnthesis, resulting in high A/g, or IWUE; in contrast, water-stressed plants would have lower stomatal density to reduce water loss, and thus lower IWUE if net photosynthesis is limited by transpiration. Differences in mean leaf area between Mountain and Coastal Plain/Piedmont plants may reflect additional habitat selection factors. Specifically, the higher mean leaf area observed in the Coastal plants, which is common in shade-plants, may reflect their more closed Pine-Savanna habitats in contrast to the exposed, sparsely-treed, high -elevation rock outcrops. Thus, although there is clear geographical variation in IWUE and related morphological traits, it remains unclear whether that variation is due to local adapation or phenotypic plasticity.

In addition to stomatal number and density, stomatal activity could influence IWUE. Lower IWUE in Piedmont/Coastal populations may be the result of plants closing their stomata early each day to cope with midday heat and water stress (e.g. Allium cepa, Quercus suber, etc.) (Meidner and Heath 1959; Tenhunen et al. 1984), thus limiting transpiration, but also reducing CO₂ uptake and daily net photosynthesis. If this were to continue into the future, Coastal plants may open their stomata for shorter and shorter periods as the climate warms, and eventually might reach a point where stomates are no longer open long enough to maintain carbon gain. The higher IWUE of Mountain plants may result from plants leaving their stomata open longer, especially in summer when lingering morning fog, cloud immersion, and the generally high humidity in the Mountain region can reduce water loss while maintaining CO₂ uptake and photosynthesis (Culatta and Horton 2014), thus producing greater daily IWUE (A/g). Ongoing climate change may threaten this trend. The Blue Ridge province consistently is predicted to experience a rise in average annual temperatures of at least 2.3°C to 3.7°C by 2100, a lengthened growing season, and although precipitation forecasts are more variable, they generally predict either historic or lower levels (2 mm to 72 mm decrease) with greater seasonal variability over the same period (McNulty et al. 2012). Cloud-base height is expected to increase and cloud immersion on outcrops to decrease (Richardson et al. 2003). Cloud-base height has already increased over three decades in the southern Appalachians (1970–2003), following the overall global trend (Chernykh et al. 2001; Foster 2001; Richardson et al. 2003; Sun et al. 2007). If periods of rain, fog, and/or cloud moisture continue to decrease or disappear from the mountains altogether, plants will become exposed to greater insolation and more frequent drought-like conditions and may begin regulating water loss like the Coastal plants with mid-day stomatal closures. Eventually, loss of daily cloud immersion may make these outcrops uninhabitable for K. buxifolia and other unique and endemic species. To test this hypothesis, observations of stomatal conductance and photosynthetic variation in these plants, including developing daily photosynthetic and transpiration curves, are needed over multi-seasonal variation in temperature and soil moisture.

In conclusion, we estimated the divergence time between K. buxifolia and its alpine sister species, K. procumbens, as late-Miocene, much later than predicted by the southern Appalachian rock-outcrop hypothesis. This divergence, as well as later population level divergences likely caused by extinctions in intervening populations, appear to coincide with major climatic shifts from the late-Miocene to mid-Pleistocene. Four potential refugia were identified, with the most ancient on the western side of the mountains at Mount LeConte, and more recent refugia on the eastern side of the mountains surrounding the Highlands Plateau, along the Suffolk escarpment, and in the New Jersey Pine Barrens. These sites should be studied as potential southeastern refugia for other mountain-coastal species with disjunct populations. Directions of postglacial colonization are difficult to infer due to relatively little within and among population variation in K. buxifolia. However, the haplotype and dating analyses suggest that the Pine Barrens populations are of Mountain (rather than Coastal) origin, becoming isolated during the Pleistocene. Population-level divergences and subsequent isolation has resulted in geographical variation in IWUE and leaf morphology, with Mountain populations exhibiting greater IWUE and stomatal density but lower leaf area compared to their Piedmont and Coastal Plain relatives. Further research is necessary to understand how and why this variation occurs, whether it is due to local adaptation or phenotypic plasticity, and the biological consequences, particularly in response to climate change. Finally, this study helps to clarify some long-held questions about the biogeographic history and trait diversity within K. buxifiolia. These results suggest that high-elevation rock outcrop communities in the Southern Appalachians may have been inhabitated by northern-affinity species like K. buxifolia for much longer than previously assumed, and that population disjunction and isolation as a result of past and ongoing climatic change may result in ecophysiological differentiation.