The role of glacial‐interglacial climate change in shaping the genetic structure of eastern subterranean termites in the southern Appalachian Mountains, USA

Abstract The eastern subterranean termite, Reticulitermes flavipes, currently inhabits previously glaciated regions of the northeastern U.S., as well as the unglaciated southern Appalachian Mountains and surrounding areas. We hypothesized that Pleistocene climatic fluctuations have influenced the distribution of R. flavipes, and thus the evolutionary history of the species. We estimated contemporary and historical geographic distributions of R. flavipes by constructing Species Distribution Models (SDM). We also inferred the evolutionary and demographic history of the species using mitochondrial (cytochrome oxidase I and II) and nuclear (endo‐beta‐1,4‐glucanase) DNA sequence data. To do this, genetic populations were delineated using Bayesian spatial‐genetic clustering, competing hypotheses about population divergence were assessed using approximate Bayesian computation (ABC), and changes in population size were estimated using Bayesian skyline plots. SDMs identified areas in the north with suitable habitat during the transition from the Last Interglacial to the Last Glacial Maximum, as well as an expanding distribution from the mid‐Holocene to the present. Genetic analyses identified three geographically cohesive populations, corresponding with northern, central, and southern portions of the study region. Based on ABC analyses, divergence between the Northern and Southern populations was the oldest, estimated to have occurred 64.80 thousand years ago (kya), which corresponds with the timing of available habitat in the north. The Central and Northern populations diverged in the mid‐Holocene, 8.63 kya, after which the Central population continued to expand. Accordingly, phylogeographic patterns of R. flavipes in the southern Appalachians appear to have been strongly influenced by glacial‐interglacial climate change. OPEN RESEARCH BADGES This article has been awarded Open Materials, Open Data Badges. All materials and data are publicly accessible via the Open Science Framework at https://doi.org/10.5061/dryad.5hr7f31.


| INTRODUC TI ON
Geographic barriers to dispersal, such as mountains and rivers, are considered major drivers of genetic divergence within and among species. The influence of climate change (e.g., glacial-interglacial oscillations during the Pleistocene) in generating phylogeographic structure is also widely recognized (Hewitt, 1996;Avise, 2000 and references therein). For example, in Europe, when ice sheets reached their maximum extent during glacials, this repeatedly resulted in range contraction into southern refugia, which subsequently served as key reservoirs for recolonization via northward expansion during interglacials (Hewitt, 1996(Hewitt, ,2004. In these regions at high latitudes, successive glacial-interglacial cycles were likely to reinforce the same genetic signatures of contraction and expansion (but see Gomez & Lunt, 2007;Shafer, Cullingham, Cote, & Coltman, 2010).
In contrast to landscapes that were repeatedly covered by ice sheet advances throughout the Pleistocene, those in temperate or tropical regions that remained unglaciated potentially contained numerous refugia (Byrne, 2008). Indeed, in montane areas with deeply dissected topography, latitude alone may be a poor proxy for the locations of refugial areas, as the steep environmental gradients that occur locally can exert a strong influence on persistence of habitat patches that can support viable populations. In such regions-in contrast to the traditional view of refuges being continuously occupied long-term stable areas-successive glacialinterglacial cycles are less likely to have repeatedly played out in the same way. Owing to stochastic processes, they may have instead been somewhat ephemeral. For instance, a refugium may have been only periodically occupied, with the process of shifting between alternative refugia from one glacial cycle to the next involving extinction at the trailing edge and colonization at the leading edge. Herein, we refer to this particular case of contraction-expansion dynamics as "distributional shift" and consider it a plausible model for the focal landscape setting. Indeed, consideration of how major shifts in geographic distributions contributed to population differentiation during the Pleistocene is important for understanding speciation processes  and references therein).
This topographically complex temperate region is characterized by steep environmental gradients, which have promoted population divergence in many species, particularly those with poor dispersal abilities (Hedin & Wood, 2002). Paleoclimatic (Loehle, 2007), biogeographic (Swenson & Howard, 2005) and comparative phylogeographic (Soltis, Morris, McLachlan, Manos, & Soltis, 2006) data indicate that the southern Appalachians remained free from Pleistocene ice sheet advances, and consequently, retained numerous refugial areas for forest-dependent biota during cool and dry glacial periods. Indeed, short-range endemism and high diversity have been well documented in plethodontid salamanders (Petranka, 1998) and other amphibians (Rissler & Smith, 2010). Similar patterns have also been reported for invertebrate groups such as crayfish (Crandall & Buhay, 2008), arachnids (Hedin & Wood, 2002;Thomas & Hedin, 2008), and millipedes (Marek, 2010). While the role of the southern Appalachian Mountains as a major barrier driving an east-west divide among lowland taxa is widely recognized (Soltis et al., 2006 and references therein), there have been surprisingly few biogeographic and phylogeographic studies of upland species that occupy the mid-and highelevation ridgelines, and research on invertebrates, in particular, is underrepresented.
The eastern subterranean termite, Reticulitermes flavipes, currently inhabits previously glaciated regions of the northeastern U.S., as well as the unglaciated southern Appalachian Mountains and surrounding areas. This species is a key ecosystem engineer that makes major contributions to dead-wood decomposition and nutrient cycling in forests (Myer & Forschler, 2019;Ulyshen, Wagner, & Mulrooney, 2014), and its distribution is influenced by humidity and temperature (Wiltz, 2015). This diploid eusocial species live in colonies that typically have a simple family structure, arising from an outbred primary reproductive pair that remains fertile for 6-11 years (Vargo & Husseneder, 2009). When the king or queen dies, some full-sib workers differentiate into male and female secondary reproductives, at which point the colony becomes inbred (Vargo & Carlson, 2006). However, in addition to temporal transitions from simple to extended families, there may also be spatial partitioning, whereby the initial reproductive center, with the primary reproductives, expands into satellite nests housing secondary reproductives (Thorne, Traniello, Adams, & Bulmer, 1999). Winged alates disperse away from the original colony and establish new colonies and then shed their wings. However, dispersal abilities are only moderate, with distances varying from a few meters to >1 km (Vargo & Husseneder, 2009). Such limited dispersal is conducive to strong historical inference (Cruzan & Templeton, 2000).
Reconstructing long-term population history is often achieved via analyses of geo-referenced DNA sequence data, using spatially explicit phylogenetic and/or coalescent-based analytical approaches (see Knowles, 2009;Hickerson et al., 2010 and references therein).
Increasingly, complementary nongenetic data are being employed to augment inferences or to generate hypotheses about past events and population processes. In particular, Species Distribution Models (SDM) are now widely used to locate glacial refugia (Richards, Carstens, & Knowles, 2007), or determine the influence of past climate change on current genetic structure (Alexandrino, Teixeira, Arntzen, & Ferrand, 2007). In some cases, similar conclusions about phylogeographic history have been drawn from SDMs and genetic data (Waltari et al., 2007). Briefly, SDMs relate occurrence records for a given species with the environmental conditions in those same locations in order to estimate geographic areas in which the species is likely to be found (Guisan & Thuiller, 2005). Given that historical climatic fluctuations can trigger range contractions and expansions-including wholesale distributional shifts (Pielou, 1991)-SDMs can form a framework for understanding the genetic consequences of glacial-interglacial climate change (Knowles & Alvarado-Serrano, 2010).
In this study, we investigated the genetic consequences of glacialinterglacial climate change on R. flavipes from the unglaciated southern Appalachian Mountains and surrounding areas, and considered distributional shifts as a plausible hypothesis (among others) to be assessed using SDMs and genetic data. Given the reliance of this species on dead-wood microhabitats, our expectation was that during the Pleistocene and earlier, R. flavipes closely tracked the changing distributions of forest habitats, and was strongly impacted by climatic fluctuations. Indeed, ecologically specialized low-mobility forest insects may be particularly well-suited for reconstructing past climatic impacts on montane forest landscapes, in part owing to their short generation times and ability to persist in habitat patches too small to support more mobile vertebrates (Garrick et al., 2004;Hugall, Moritz, Moussalli, & Stanisic, 2002;Sunnucks et al., 2006). Furthermore, owing to the limited dispersal ability of R. flavipes, we expected that relatively fine-scale genetic structuring would be detectable. To test these expectations, we modeled present and past distributions and used contrasts between these SDMs to make inferences about distributional shifts and to identify areas of stability (i.e., potential refugia). Based on this, we generated competing hypotheses about drivers of genetic divergence, and then tested these via analyses of DNA sequence data using coalescent simulations. In addition to the effects of historical climatic conditions, we also considered the influence, if any, of contemporary climatic conditions and dispersal-based spatial structure on genetic variation in R. flavipes.

| Phylogeographic framework
To address the aims of this study, we used the following workflow: Step 1-Model present-day and historical climate-based distributions of R. flavipes in order to identify potential refugia and generate expectations about directionality of range contractions or expansions, including distributional shifts; Step 2-Infer the number of distinct populations using spatial-genetic clustering, and cross-validate via principal component analysis, and phylogenetic reconstruction; characterize genetic variation within and differentiation among populations, and; estimate the amount of genetic variation explained by dispersal (spatial structure) and environment (contemporary climatic conditions); Step 3-Test alternative phylogeographic hypotheses to determine whether expansion out of refugia, distributional shifts, or vicariance was the underlying historical process generating the observed patterns of genetic variation within and among populations; estimate values of parameters included in the best-fit phylogeographic hypothesis; and assess evidence for changes in effective population size over time.  (Wang et al., 2009), termites were identified using a molecular assay (Garrick, Collins, Yi, Dyer, & Hyseni, 2015). Ultimately, R. flavipes were sampled from 50 rotting logs across 46 locations ( Figure 1; also see Supporting Information Table S1.1). From each log, 1-3 individuals were used for phylogeographic analyses. For out-group taxa, we included specimens representing three close relatives (Supporting Information Table S1.2): R. virginicus (n = 3 individuals), R. malletei (n = 1) and R. nelsonae (n = 1).

| Genetic data collection
Extraction of genomic DNA was performed using a DNeasy tissue kit (Qiagen, Valencia, CA) following the manufacturer's recommendations. Portions of the mitochondrial cytochrome c oxidase subunit I (COI) and II (COII) genes, and an intronic portion of the nuclear endo-beta-1,4-glucanase (EB14G) gene, were amplified via polymerase chain reaction using primers (Supporting Information   Table S2.3) and conditions reported in Supporting Information Appendix S2, and then sequenced at Yale University. Sequence alignments were performed using Geneious v.6.1.8 (Kearse et al., 2012), and manually edited as necessary.
We concatenated COI and COII and refer to this sequence (COI + COII) as the mitochondrial DNA (mtDNA) locus; we refer to EB14G as the nuclear DNA (nDNA) locus. For the latter, heterozygous sites were scored using the "Find Heterozygotes" plugin in Geneious. For a site to be considered heterozygous, we required that height of the secondary peak was at least 50% of the primary peak (sites with quality scores < 20, were coded as "N"). Allele haplotypes were inferred using PHASE v.2.1.1 (Stephens, Smith, & Donnelly, 2001), with the following settings: 90% phase certainty, 10,000 iterations, thinning interval = 10, burn-in = 1,000, and the default recombination model. PHASE was run three times to evaluate consistency of results.

| Step 1: Present and past geographic distributions
There are few published occurrence records of forest populations of R. flavipes with confirmed species-level identifications and adequate geospatial precision for SDM. Accordingly, in addition to the 46 sites that contributed to genetic analyses (above), the presence of R. flavipes at an additional 45 locations (surveyed from 2015 to 2016) was confirmed using Garrick et al.'s (2015) molecular assay, resulting in a total of 91 occurrence points (Supporting Information Table S1.1). To construct SDMs, we used the "biomod2" package (Thuiller, Georges, Engler, & Breiner, 2016) in R (R Core Team, 2018). Full details about SDM construction are given in Supporting Information Appendix S3.
Briefly, we used four machine learning algorithms to model distributions based on climatological data, presence records, and 20 independent sets of 100 pseudo-absence points ( Figure 2). The latter choice was based on work by Barbet-Massin, Jiguet, Albert, and Thuiller (2012), who showed that for machine learning methods it is better to use multiple replicates of pseudo-absence points, with the number of pseudo-absences in each replicate close to the number of occurrence points. We used environmental variables at a 1-km resolution for SDM construction. Present-day SDMs were based on mean climatological data spanning , and historical distributions were modeled for the Mid-Holocene (MH; ~6 kya), the F I G U R E 2 Species Distribution Modeling. Diagram showing the conceptual framework used to generate SDMs that enabled contrasts between successive time periods: "present" , Mid-Holocene (MH; ~6 kya), Last Glacial Maximum (LGM; ~22 kya), and Last Interglacial (LIG; ~120-140 kya)

| Distributional shifts and areas of stability
We used a threshold value to convert continuous occurrence probabilities to a binary classification of suitable (>0.2) versus unsuitable (≤0.2). The occurrence probability threshold was chosen based on the True Skill Statistic (TSS; Allouche, Tsoar, & Kadmon, 2006).
Specifically, we chose a threshold value that maximized the TSS, as this approach has consistently performed better than other thresholding methods (Liu, Berry, Dawson, & Pearson, 2005;Liu, Newell, & White, 2016;Liu, White, & Newell, 2013). However, since we used multiple pseudo-absence replicates, we had the opportunity to maximize TSS without risking under-prediction of presences, which results from choosing a high threshold value. Indeed, using distributions of TSS and threshold values, we were able to select the lowest threshold (0.2; Supporting Information Figure S3.1), below which TSS had a steep slope. To calculate the distributional shift between two successive time periods (e.g., LIG to LGM, or LGM to MH), we took the difference of the two binary maps, after multiplying the more recent time period by two in order to ensure that we obtain four categories in the distributional shift calculation: colonization (difference = 2), stability (1), absence (0), and extinction (−1; see Supporting Information Figure S3.2). Similarly, to estimate areas of stability (i.e., persistence in a location between successive time periods), we multiplied the binary occurrence maps (Supporting Information Figure S3.2) of the corresponding periods: locations where the product is 1 were considered to harbor stable habitats across time periods (stability = 1).
We assessed values of K (i.e., the number of clusters) ranging from 2-20, with 10 replicate runs each. We also examined evidence for geographically cohesive genetic groups by representing the variance in mtDNA sequences using Principal Components Analysis (PCA), performed with the "prcomp" function in R.

| Phylogenetic reconstruction and molecular dating
We reconstructed a mtDNA-based dated phylogeny to verify the existence of any genetic groups determined by BAPS, as well as to estimate divergence times. First, we used PartitionFinder 1.1.0 (Lanfear, Calcott, Ho, & Guindon, 2012) to determine the best partitioning scheme, and the Bayesian Information Criterion in jModelTest v.2.1.10 (Darriba, Taboada, Doallo, & Posada, 2012) to identify the optimal model of sequence evolution. The best-fit model for all three codon positions was HKY + I (Hasegawa, Kishino, & Yano, 1985). Then, to estimate a dated phylogeny, we used BEAST v.2.4.5 (Bouckaert et al., 2014), with a relaxed log-normal molecular clock (Drummond, Ho, Phillips, & Rambaut, 2006), and a coalescent tree prior. We used broad mutation rate priors. For the mtDNA locus, the range included Brower's (1994) commonly used insect rate of 1.15% sequence divergence per lineage per million years, and Luchetti, Marini, and Mantovani (2005) faster rate of up to 140% per million years, which was estimated from COII in European Reticulitermes taxa. Based on point estimates obtained using approximate Bayesian computation (ABC; Beaumont, Zhang, & Balding, 2002) assessments of competing phylogeographic hypotheses (described in Section 2-Step 3), we set the mean mutation rate at 12% per million years (see Section 3-Step 3) for the mtDNA locus. Since there was no mutation rate information available for the nDNA locus in Reticulitermes, we estimated the mean mutation rate in BEAST by conditioning on the mtDNA locus and setting the initial mean value at 0.6% with a range of 0.2%-2% (obtained using ABC; see Section 3-Step 3). BEAST was run for 50 million Markov chain Monte Carlo generations, with samples saved every 2,500 generations, after discarding the first 5 million generations as burn-in. We used Tracer 1.6 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018) to examine the stationarity of parameter estimates and to determine that effective sample sizes were greater than 500.
BEAST was run with and without the out-group Reticulitermes taxa using the same settings. Results were summarized via a Maximum Clade Credibility tree in TreeAnnotator v.2.4.4 (Bouckaert et al., 2014), with the first 25% of trees discarded as burn-in.

| Genetic variation influenced by environment and dispersal
To estimate the amount of genetic variation explained by spatial structure versus the environment, we used distance-based redundancy analysis (dbRDA; Legendre & Anderson, 1999). We computed the genetic distance matrix using the "dist.dna" function of the "ape" (Paradis & Schliep, 2019) package, and performed dbRDA using the "capscale" function of the "vegan" (Oksanen et al., 2018) package in R. To compute the response variable, genetic distances (i.e., matrix of pairwise mutational differences between DNA sequences) were estimated using the TN93 (Tamura & Nei, 1993) model of sequence evolution, allowing for different rates for transitions and transversions. For environmental predictors, we used the contemporary environmental factors obtained via factor analysis (see Section 2-Step 1). To obtain spatial structure predictors, we transformed Euclidean geographic distances to a continuous rectangular vector by Principal Coordinates analysis of Neighbor Matrices (PCNM) using the "pcnm" function in "vegan." Significance of the predictors was assessed using multivariate Fstatistics with 9,999 permutations. We first analyzed the relationship between the genetic distance matrix and each environmental factor separately, and then performed a partial dbRDA for each variable while controlling for the influence of spatial structure, using only significant PCNM eigenvectors. Similarly, we analyzed the relationship between genetic distances and PCNM eigenvectors, retained the significant eigenvectors, and then removed interactions with the environment to obtain the contribution of spatial structure alone.

| Competing scenarios
We used ABC, as implemented in the software DIYABC v.2.1.0 (Cornuet et al., 2014), to assess alternative hypotheses designed to determine whether expansion out of long-term stable refugia, distributional shifts, or vicariance was the major underlying process generating the present-day spatial distribution of genetic variation. mtDNA plus (phased) nDNA sequence data were used, and we conditioned these analyses on a posteriori knowledge of the existence of three distinct genetic clusters of R. flavipes (see Section 3-Step 2).
Because ABC analyses can suffer when a large number of candidate models are simultaneously considered (Pelletier & Carstens, 2014), we employed a two-tiered approach, where best-fit scenarios from separate analyses in the first tier are subsequently compared against each other in the second tier. This hierarchical or tournament-style approach has also been applied in other study systems (Espindola et al., 2016;Stone et al., 2017).
All scenarios in both tiers incorporated bottleneck events, because they all involved divergence of new populations from an existing population, and thus founder effects. Indeed, our inclusion of bottleneck events enabled specification of progenitor-descendant relationships between pairs of diverging populations (as in Garrick et al., 2014). Furthermore, the nonnegligible role of bottlenecks during climatically driven population divergence has been established. In one set of analyses in the first tier of ABC comparisons, we assessed scenarios in which R. flavipes persisted in a single major refugium (Supporting Information Figure S4.3), such that the other areas were colonized via successive expansions out of that refugium. We considered three different refugial locations (i.e., the north, south, or central portion of the study region; see Section 3-Step 2). In a second set of analyses within the first tier, we assessed scenarios that involved distributional shifts (Supporting Information Figure S4.4), whereby populations diverged in a stepping-stone fashion (i.e., one population gave rise to a descendant population, which later became the progenitor of the third population). Here, we considered all possible stepping-stone configurations (i.e., there was no assumption that only nearest neighbors can exhibit a progenitor-descendant relationship). In the second tier of ABC comparisons, the best-fit hypotheses from the refugial and distributional shift scenarios were directly compared, along with an additional hypothesis that incorporated vicariance (Supporting Information Figure S4.5). The reason for including this third hypothesis was to test the possibility that the original ancestral population no longer exists, having split into two new populations, one of them giving rise to a third population.
While there are other vicariance hypotheses that could have been compared in the first tier, we chose not to do this based on the sequence of divergence events best-fit refugial and distributional shift hypotheses had in common. This reduced the number of plausible vicariance hypotheses to one.

| ABC model specification, and model choice
Within the ABC framework, two classes of model parameters were used to characterize the phylogeographic hypotheses described above: effective population sizes (Ne), and divergence times (T).
We performed two rounds of modeling: (a) a preliminary round with broad priors, and (b) the final round with narrower priors (Supporting Information Table S4.5). Briefly, all competing scenarios had two divergence events: any two of T N , T C or T S, (where the subscript is the first letter abbreviation of the new cluster, i.e., Northern, Central, or Southern), the prior range for the more recent event encompassed the MH and the LGM whereas priors for the older event ranged from the LGM to the LIG assuming a 1-year generation time for R. flavipes. Full details of ABC priors on Ne and T parameters are given in Supporting Information Appendix S4. We set the mtDNA mutation rate priors from 5.0 × 10 −9 to 5.0 × 10 −7 , a broad range encompassing the Brower (1994) and Luchetti et al. (2005) rates (see Section 2-Step 2). Similarly, since no rates were available for the nDNA locus in Reticulitermes, we used broad priors for this locus, from 5.0 × 10 −10 to 2.5 × 10 −8 . Thus, the mean nDNA rate was an order of magnitude slower than the mean mtDNA rate, despite some overlap at the upper end of nDNA and lower end of mtDNA prior ranges.

HYSENI aNd GaRRICK
To characterize the empirical two-locus DNA sequence dataset, we used the following summary statistics: number of segregating sites (one-and two-sample) and private segregating sites (one-sample), mean (one-and two-sample) and variance of pairwise differences (one-sample), mean and variance of numbers of the rarest nucleotide at segregating sites (one-sample), Tajima's (1989) D (onesample), and F ST (Hudson et al., 1992) between two samples. ABC runs consisted of 1 × 10 6 simulated genetic datasets per competing phylogeographic hypothesis. We then compared the values of summary statistics calculated from simulated datasets to those from the empirical dataset. Following Cornuet et al. (2014), model checking was performed via principal components analysis, and then posterior probabilities were calculated via logistic regression (Fagundes et al., 2007) on 1% of simulated data most similar to the empirical data, to identify the best-fit model (Cornuet et al., 2008). We evaluated model performance (i.e., the ability to discriminate between the bestsupported and alternative scenarios), by estimating type I and type II error rates. To do this, we simulated 500 data sets and estimated the most likely model using a polychotomous logistic regression (Cornuet, Ravigne, & Estoup, 2010;Cornuet et al., 2008). The type I error rate was the proportion of data sets that were simulated under an alternative scenario but were incorrectly categorized under the best-supported scenario. The type II error rate was the proportion of instances in which the best-supported scenario was incorrectly selected as the most likely scenario. To calculate point estimates and confidence intervals for the values of parameters included in the best-fit model, we selected 1% of the simulated data closest to the observed data. Additionally, for the best-fit scenario, we estimated precision in parameter estimation (Cornuet et al., 2010) by computing the relative median of the absolute error for 500 simulated data sets with values drawn from posterior distributions.

| Population size changes over time
For each of the three R. flavipes genetic groups, we assessed evidence for population size changes versus stability by calculating Tajima's D, and Fu and Li's D* and F* (Fu & Li, 1993) from the mtDNA data, in DnaSP. To identify cases of departure from the null hypothesis of constant size, p-values for these statistics were obtained by computing 10,000 coalescent simulations based on θ from the observed data and assuming no recombination. We also calculated Ramos-Onsins and Rozas' (2002) R 2 statistic for which significantly small values indicate population growth, whereas significantly large R 2 values indicate size reduction. Statistical significance of deviation from the null hypothesis of constant population size was assessed by performing 10,000 coalescent simulations in DnaSP. To complement the above analyses, we also estimated mismatch distributions, where a unimodal distribution indicates growth, whereas a multimodal distribution is indicative of size constancy (Rogers & Harpending, 1992). Given that signatures of selection can mimic those of population size changes and therefore complicate interpretation of the above summary statistics, we examined evidence for nonneutrality using compound tests (Zeng, Shi, & Wu, 2007). We performed the compound tests using the program DH (http://zeng-lab.group.shef. ac.uk/wordpress). The significance (α = 0.05) of each test was determined using 100,000 simulations.

| Step 1: Present and past geographic distributions
When constructing SDMs, a strong correlation was observed among some of the 19 bioclimatic variables (Supporting Information Figure   S5.

| Step 2: Genetic variation and the role of environment and space in genetic structuring
The BAPS analysis identified three genetic clusters, each with largely separate geographic distributions (Figure 4a). Herein, we refer to them as the Northern, Central, and Southern clusters. We used the first three principal components (PCs) to represent these clusters in three dimensions (Figure 4b). The three PCs accounted for 53% of the variance at the mtDNA locus; they showed that the Northern cluster is most similar to the Central cluster. Phylogenetic reconstruction using BEAST produced a Bayesian tree (Figure 4c)  Although the Southern cluster comprised only four mtDNA haplotypes, this group had the most genetic variation (nucleotide diversity, π = 0.016; mean number of nucleotide differences, K = 17.50; (p = 0.012) of the variation. The interaction between the two explained an additional 12.4% of the genetic variation. After removing the effect of spatial structure, the factors with significant contribution to genetic variation were "temperature range" and "wet-season precipitation" (Supporting Information Figure S6.11).

| Step 3: Phylogeographic hypothesis testing and population size changes
In the two sets of first-tier ABC comparisons: (a) the refuge-based scenario with the highest posterior probability was the hypothesis that postulated the Northern region was the source from which the Southern cluster diverged first, followed by the Central cluster (scenario R3; Table 2; Supporting Information Figure S4.3); and (b) the distributional shift scenario that provided the best fit to the empirical data was the hypothesis that represented a case of Southernto-Northern-to-Central stepping-stone colonization (scenario DS1; Table 2; Supporting Information Figure S4.4). In the second tier of ABC comparisons, the best-fit scenario was DS1 (Table 2; Supporting Information Figure S4.5    Best-fit scenarios are highlighted in bold font. Note. ABC hypothesis testing was performed in two tiers. In the first tier, refugial and distributional shift scenarios were evaluated separately. In the second tier, these two scenarios, as well as a vicariance scenario (V; Supporting Information Figure S5.9), were compared.
Indeed, the predominant focus on vertebrates and vascular plants in conservation research and planning is likely to result in management strategies that fail to cater to a large proportion of biodiversity (Garrick, Newton, & Worthington, 2018 and references therein).
To understand drivers of phylogeographic patterns in R. flavipes, we examined evidence for distributional shifts using SDMs, and reconstructed the evolutionary and demographic history of R. flavipes using ABC analyses. Overall, we determined that the location of key refugia has changed over time (e.g., from one glacial period to the next), rather than a single refugium repeatedly serving as a reservoir of genetic diversity, whereby successive glacial-interglacial cycles reinforce the same genetic signatures of contraction and expansion.

| Climate change as a driver of distributional shifts and genetic divergence
Determining whether distributional shifts have occurred in the history of a species can lead to a better understanding of processes that have shaped present-day genetic variation. Our SDMs sug-  (Merz et al., 2013). Similarly, the red salamander, Pseudotriton ruber, persisted in the Coastal Plain in the early Pliocene, and then expanded its range toward Appalachian upland habitat as cooling trends started in the early Pleistocene (Folt, Garrison, Guyer, Rodriguez, & Bond, 2016 Despite the broad geographic range of the R. flavipes Central cluster (Figure 4a), this group contained the lowest genetic diversity (Table 1). We suggest that this is likely the result of founder effects associated with the relatively recent colonization of the central portion of the southern Appalachians from the north. Although subsequent population expansion seems to have occurred in the central region, more time may be needed to replace lost genetic variation.
Assessment of changes in Ne over time showed that the Central cluster had increased in size over the last 10,000 years (Figure 5b), which is consistent with inferences based on nongenetic data that indicated the amount of suitable habitat in the central region increased since the LGM (Figure 3).

| The potential role of environmental variables in promoting range expansions
Given the desiccation susceptibility of soft-bodied arthropods, range expansions and population growth in R. flavipes may be have been influenced by local-scale site-specific environmental variables such as precipitation. The southeastern U.S. was much warmer during the mid-Holocene (cf. LGM; Bartlein et al., 1998), when tupelo and oak forest types dominated over pine, indicating wetter conditions (LaMoreaux, Brook, & Knox, 2009). The R. flavipes Central cluster likely diverged from the Northern cluster following a cooling trend in the Younger Dryas (~12.9-11.7 kya). While this was a global cooling period, locally in the southeastern U.S., this period was characterized by a warmer and wetter climate, reflecting the trapping of heat in the western subtropical gyre due to reduced Atlantic meridional overturning circulation (Grimm et al., 2006). Accordingly, if high precipitation was important for facilitating range expansion, these conditions seem to have been in place at a time that coincides with colonization of the central region. Furthermore, seasonal differences in precipitation between the southern and northern portions of the study region (Hyseni & Garrick, 2019)  This contrast supports the view that fine-scale genetic structuring may be particularly prevalent in topographically complex montane areas (Garrick, 2011;Hedin & Wood, 2002;Thomas & Hedin, 2008).
Along a ~1,000 km transect traversing the southern Appalachians, a wood-feeding cockroach (Cryptocercus punctulatus) that is syntopic with R. flavipes consists of five distinct genetic groups (Everaerts et al., 2008;Garrick, Sabree, Jahnes, & Oliver, 2017). Interestingly, both of these saproxylic taxa have a zone of parapatry between genetic groups in the central region. Comparative phylogeographic analyses would be informative about the extent to which spatial-genetic patterns seen in dead-wood-associated insects correspond with shared microevolutionary processes that underpin them.

| Caveats and future directions
An early understanding of genetic consequences of Pleistocene range expansions came from study systems that either repeatedly experienced severe glaciation (Hewitt, 1996), or were relatively simplified linear systems (Nason, Hamrick, & Fleming, 2002). In these cases, unidirectional expansion out of a single major refuge was commonly inferred, often based on signatures of repeated founder effects and serial reduction in genetic diversity at the leading edge. However, an expanded view of the geography of range expansion may be needed when considering unglaciated, topographically complex, montane landscape settings. In this study, we considered distributional shift (see Section 1) to be a plausible phylogeographic scenario for the southern Appalachian Mountains.
However, further work is needed to understand the circumstances under which distributional shift scenarios are distinguishable from single-refuge contraction-expansion scenarios. Indeed, inferring Pleistocene distributional shifts using genetic data can be challenging, as multiple historical factors can contribute to current genetic variation.
Although our ABC analyses identified distributional shift as the bestfit scenario, it did not receive unambiguously superior support relative to the next-best scenario, and the estimated error in scenario choice was large (Supporting Information Table S7.8). Accordingly, we must consider our ABC-based inference to be a preliminary working hypothesis, to be re-evaluated and re-tested with new data. Notwithstanding some limitations of our ABC inferences, it is notable that a common feature of the best-fit and second-best hypotheses is the expansion of the Central cluster. Specifically, both scenarios include the Northern cluster giving rise to the Central cluster. Additionally, both scenarios include a direct long-distance dispersal event. Buckley (2009) advocated for an iterative approach to phylogeography, highlighting the value of working hypotheses for focusing subsequent analytical efforts on scenarios that have some empirical support. This study contributes to a growing body of literature that highlights an important role for multiple refugia-including those located further north than previously expected-in phylogeographic structuring of plants (McLachlan et al., 2005), vertebrates (Fontanella, Feldman, Siddall, & Burbrink, 2008), and invertebrates (Merz et al., 2013). Having characterized contemporary fine-scale spatial structure and historical climate-based distributions for R. flavipes, the present study has also revealed specific geographic locations that warrant dedicated sampling (e.g., the Southern genetic cluster has a relatively small range that requires better representation, and based on SDMs, sampling in the Gulf Coast and Coastal Plain areas would be particularly valuable).

ACK N OWLED G M ENTS
We thank V. Noguerales and O. Razgour, whose valuable comments helped improve the manuscript. We thank R.J. Dyer, R.E. Symula, and

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
Both authors conceived the study, and collected and curated samples from across the southern Appalachians; C.H. collected DNA sequence data, performed the analyses, and wrote the first draft; both authors edited subsequent drafts.