Fine- scale genome- wide signature of Pleistocene glaciation in Thitarodes moths (Lepidoptera: Hepialidae), host of Ophiocordyceps fungus in the Hengduan Mountains

The Hengduan Mountains region is a biodiversity hotspot known for its topologi - cally complex, deep valleys and high mountains. While landscape and glacial refugia have been evoked to explain patterns of interspecies divergence, the accumulation of intra- species (i.e


| INTRODUC TI ON
Climatic oscillations from the late Pleistocene (130-11 kya) glaciation and subsequent Holocene warming left their signatures on the genetic structuring of many species around the world (Bagley et al., 2017;Devitt et al., 2013;Dudaniec et al., 2012;Marske et al., 2020;Wachter et al., 2016;Yannic et al., 2014). The signature of population responses to glacial expansion and subsequent retreat since the last glacial maximum (LGM, 26 kya) are especially well-documented. During periods of extensive glaciation, large ice sheets sometimes acted as direct barriers of gene flow that caused population divergence (Marske et al., 2020). While some species persisted in periglacial refugia (Escobar García et al., 2012), others remained at high-elevation ice-free "nunataks", resulting in deep divergence between populations (Stehlik et al., 2002). Post LGM, some populations colonized newly available ice-free habitats and showed a signal of immediate population expansion (Wachter et al., 2016); or some populations declined due to founder effects (Wang et al., 2020). These demographic responses are not mutually exclusive; the genetic signatures of glaciation are often taxa and region-specific. To date, the majority of studies on genetic signatures of Pleistocene glaciation have been conducted in North America and Northern Europe (for review see Wallis et al., 2016).
By contrast, the eastern Himalayas and Hengduan Mountains house an estimated 40,800 km 2 of existing glaciers (which account for up to 70% of the seasonal water flow of the region, Xu et al., 2009;Bolch et al., 2012), but the genetic signatures of historical glaciation on divergent populations of this region have rarely been investigated.
The Hengduan Mountains region (also known as the "eastern Himalayas" or the "mountains of southwest China") is a global biodiversity hotspot (Marchese, 2015). The region is known for topologically complex deep valleys and high mountains. In explaining the extraordinary Rhododendron diversity of the region (more than 50% of the world's existing species), Irving and Hebda (1993) postulated that species inhabiting a deep valley diverged in response to recurrent glacial-interglacial climate changes. While insect diversity in the region is largely understudied (Ji et al., 2013;Liu et al., 2020), recent studies on the distribution of plants and vertebrates confirmed the presence of a Pleistocene glacial barrier and refugia across the entire range of Hengduan-Himalaya (Meng et al., 2015;Qiu et al., 2011;Wan et al., 2021). Most studies of biodiversity in the region have addressed the process of interspecies divergence, while the "missing link" lies in explaining the accumulation of intraspecific genetic diversity across the mountain-valley landscape.
We here focus on the process of population divergence of the ghost moth, Thitarodes shambalaensis (Lepidoptera: Hepialidae, , which is endemic to pockets of high elevation alpine meadows inside the seven valleys of Mt. Gongga (29°35′45″N, 101°52′45″E), the highest peak (7556 m) in the Hengduan Mountains ( Figure 1). In the following introduction we first review how genetic divergence is spatially and temporally-explicitly modelled in complex landscapes, then we introduce how these methods are relevant to our specific study system of T. shambalaensis at Mt. Gongga using historical moraine records.

| Measuring population divergence across complex landscapes
Landscape genetics examine how geographical and environmental features shape population structure along spatial and temporal axes (Manel et al., 2003). On the spatial axis, geographical and environmental variation may limit an individual's ability to disperse (Kidd & Ritchie, 2006). The resulting genetic divergence among distant populations due to reduced gene flow and random drift is called "isolation by dispersal limitation" (IBDL). Populations under IBDL exhibit patterns of isolation-by-distance (IBD), whereby genetic differences between populations are positively correlated with their geographic Euclidean distances (Wright, 1931) or log-transformed equivalent in 2D models (Rousset, 1997). Since IBD is not the result of natural selection, it can be readily measured with genetically neutral markers such as microsatellites (Vieira et al., 2016) and restriction site-associated DNA sequencing (RADseq, Andrews et al., more similar to each other than those in different environments. Empirically, populations that are separated by geographical distances (IBD/IBR) often inhabit different environments (IBE), making it difficult to determine which process drove the diversification within a given species. Qualitative methods have been developed to assess the relative contribution of environment and geography to genetic structure (Bradburd et al., 2013;. A number of studies have found that a combination of IBD/IBR and IBE best explains existing genetic variation among populations (Bagley et al., 2017;Dudaniec et al., 2012;Fant et al., 2014;Rutherford et al., 2021;Smith et al., 2019;Van Buskirk & Jansen van Rensburg, 2020).
The second, temporal axis of landscape genetics, emphasizes that extant population differentiation is the historical average of gene flow over a long period of time (Bradburd & Ralph, 2019).
Prominent topological features of the landscape such as drainage systems and mountains set hard boundaries of population gene flow (Mayr, 1942), while historical fluctuations in climate alter the habitable range, and thus the resistance landscape, of gene flow (Vasconcellos et al., 2019). Populations often exhibit patterns of isolation-by-history (IBH), in which contemporary patterns of genetic structure (or measurable gene flow between populations) reflect historical habitat ranges or orogenesis boundaries that are no longer present (Marske et al., 2020;Moreira et al., 2020;Silva et al., 2021;Vasconcellos et al., 2019).
The dual spatial and temporal considerations in landscape genetics call for models of gene flow that are explicit in the distribution of populations, their pedigrees, and their environmental landscapes past and present (Bradburd & Ralph, 2019). Some recent studies specify historical population distributions and dispersal landscapes in agent-based, spatially-explicit simulations (Becheler et Currat et al., 2019) and leverage approximate Bayesian methods (Wegmann et al., 2010) to select for historical scenarios that best explain empirical genetic structure (Bemmels et al., 2016;Brown & Knowles, 2012;He et al., 2013;Lacey Knowles & Alvarado-Serrano, 2010;Massatti & Knowles, 2016). Studies also use coalescent-based methods (Excoffier et al., 2013) that start with current genetic compositions and simulate pedigree backwards in time, with prior parameters that reflect gene flow restrictions due to historical events (Bagley et al., 2017;Lanier et al., 2015;Moreira et al., 2020;Myers et al., 2020). Both the spatially-explicit simulation approach and the coalescent approach often combine genomic data with environmental niche modelling (ENM) to infer potential population distributions and barriers to gene flow (Carstens & Richards, 2007), and are generally referred to as integrated distributional, demographic and coalescent (iDDC; He et al., 2013) modelling approaches.
These spatially and temporally explicit genomic simulations are usually conducted on populations with wide-ranging, welldocumented distributions in environmentally variable habitats. At more refined spatial scales (e.g., populations within 10 km of each other), the resolution of paleoclimate data is not sufficient for ENM.
Paleoclimate data for ENM are limited to 2.5 arc min (4.6 km at 0 latitude) for most time periods (Fordham et al., 2017), and at best 30 arc sec (0.93 km at 0 latitude) for the last glacial maximum (Karger et al., 2017). If genetic structure exists at a finer scale than ENM allows, researchers often rely on known changes in immediate environment or historical events to formulate their coalescent hypotheses (Geng et al., 2009;Helsen et al., 2015;Hevroy et al., 2018;Kitamura et al., 2018;Kobayashi et al., 2018;Pearson et al., 2020;Ziege et al., 2020).

| Study system
We focus on the process of population divergence across seven glacial valleys of Mt. Gongga ( Figure 1). The valleys drain into the Moxi river, which runs along a fault line (Kangding-Moxi fault) east of Mt.
Gongga. The valleys extend 30 km from north to south, with parallel neighbouring valleys less than 10 km from each other. ENM at this scale is difficult, but the extent of ice sheet retreat from the last local glacial maximum to the present day has been inferred by moraine dating (Li et al., 2010;Pan et al., 2012;Wang, Pan, et al., 2013). The local glaciation history relevant to this study includes: (1) the local last glacial maximum (LGM L ) at Mt. Gongga occurred about 58 kya, considerably earlier than the 26 kya timing of the "global" LGM in Europe and North America (Clark et al., 2009). At LGM L , an ice sheet extended from the peak of Mt. Gongga to the Kangdi-Moxi fault. (2) The path of retreat of the post-LGM L ice sheet overlaps with creation of present day river valleys, suggesting that these valleys were gradually carved out during the retreat process. These regional histories allow us to test specific hypotheses of population divergence in the region.
Populations of the ghost moth, Thitarodes shambalaensis, (Lepidoptera: Hepialidae, Wang, Zhuang, et al., 2019) are endemic to pockets of high elevation alpine meadows inside the seven valleys of Mt. Gongga (Table S1). The full life history of T. shambalaensis has never been documented, but is likely to be similar to other species of Thitarodes that live in high elevation meadows (3600-5000 m) of the Himalaya and Hengduan Mountains (Chu et al., 2004;Wang & Yao, 2011). These moths have strict habitat requirements (Hopping et al., 2018;Wang et al., 2020) and low vagility due to their short period of adult activity (Maczey et al., 2010;Zou et al., 2011): larvae develop under the soil and can take 1-1.5 years to become adults (Tao et al., 2016), which emerge, mate and lay eggs over a brief 1-2 week period in the summer months. Thitarodes larvae are sometimes parasitized by the fungus, Ophiocordyceps sinensis, (order Hypocreales; family Ophiocordycipitaceae), resulting in a fungal mycelium-insect exoskeleton complex commonly called "caterpillar fungus" (Sung et al., 2007;Shrestha et al., 2014). Caterpillar fungi have been collected across the Himalaya-Hengduan region as ethnobotanical medicine since the 15th century (Winkler, 2008(Winkler, , 2010. Recent molecular evidence suggests that Thitarodes across the Himalaya-Hengduan region diverged into roughly 50 species along major landscape barriers and coevolved with their parasitic fungal strains (Quan et al., 2014;Zhang et al., 2014). However, fine-scale population-level divergence within a single species of Thitarodes has yet to be investigated. The distribution of T. shambalaensis populations in glacial valleys of Mt.
Gongga is interesting because contemporary populations, each located in a small "pocket" or cul de sac at the end of a glacier (i.e., at the head of a river valley), are separated by less than 10 km of Euclidean distance (Figure 1b), but the mountains flanking each valley and dividing each population are likely to be major dispersal barriers ) that confine isolated populations within a relatively small valley habitat. We hypothesized that T. shambalaensis populations arrived at their current isolated habitats from a more connected landscape following the retreating ice sheet post- LGM L , and that their contemporary genetic structure should reflect a combination of landscape topology (IBR) and historical glaciation events (IBH).
We used genome-wide SNPs to investigate population struc- To do this, we incorporated moraine-based historical distance sets as migration parameters in coalescent simulations. We then evaluated whether these refinements provide better estimates of contemporary genetic structure. To further explore the demographic dynamics derived from moraine-based models, we also modelled the divergence time and migration time of each population separately.
We tested the hypothesis that populations of T. shambalaensis underwent historical migration around LGM L , followed by population expansion in relative isolation during the subsequent glacial retreat. Demultiplexed sequences were cleaned by removing reads that had Phred quality scores lower than 10 for over 50% of their base pairs. Reads were assembled de novo following the STACKS pipeline (Rochette & Catchen, 2017). We started by grouping reads with a minimum coverage depth from a single individual into unique stacks.

| MATERIAL S AND ME THODS
We then removed stacks that were repetitive (using the program ustacks). We called SNPs on each locus following a multinomial-based likelihood model (Catchen et al., 2013), and compiled a catalogue of consensus loci among all individuals (program cstacks). We targeted a minimum coverage of three reads for a stack and allowed a two base pair mismatch between each stack (a conservative approach to delineating loci, see Rochette and Catchen (2017). Stacks with no more than one base pair mismatch across individuals were merged into consensus loci. Alleles from each individual were then called against the consensus loci catalogue in the program sstacks. To identify missing data, we used an approach similar to that of de Medeiros and Farrell (2018) to visualize the presence/ absence of each SNP for all individuals after hierarchical clustering (Ward, 1963). The resulting matrix ( Figure S1) identified individuals with highly degraded DNA as well as those with enzymatic cut sites not shared by most T. shambalaensis. We also cross-referenced these data with an independently derived CO1-based phylogeny to remove any misidentified individuals that were not monophyletic with T. shambalaensis (see Wang, Zhuang, et al., 2019 for details). After removing these individuals, we used the "-r" option of the populations program in stacks to select loci with less than 10% missing data across all individuals.
We included 1 SNP per locus to minimize linkage disequilibrium (LD) between SNPs.

| Population structure and phylogeny
We calculated the nucleotide diversity (θ, estimated as π), observed and expected heterozygosity (H obs and H exp ) and inbreeding coefficient (F IS ) for each population using the populations program in STACKS to examine possible correlations between environmental variables and population genetic diversity (Table 1), whereas these population diversity-related variables might correlate with historical habitat stability (Yannic et al., 2014;Ortego et al., 2015). Nei's pairwise F ST between all populations was calculated with the pairwise.fst function in r package adegenet (Jombart, 2008). To test for IBD, we applied Mantel tests (Mantel, 1967) to look for any correlation between population pairwise F ST and geographic Euclidean distance. While fine-scale climatic variables are not available at the resolution of individual habitats, we used elevation as proxy for environment (highly correlated with most bioclimatic variables, see supplementary information in Hopping et al. (2018) and tested for IBE by measuring the Mantel correlation between population pairwise F ST and elevational differences.
We conducted three different analyses to investigate the genetic structure of T. shambalaensis ( Figure 1c, Figure 2). First, we con- We iterated from 1 to 20 clusters and selected the number of genetic clusters with the lowest Bayesian Information Criterion (BIC) score. We visualized the selected clustering result in a principle component analysis (PCA) for SNP data (Patterson et al., 2006).
Multivariate approaches to clustering such as DAPC make no assumptions about ancestral populations. Secondly, we estimated the ancestry coefficients of each individual using sparse non-negative matrix factorization algorithms (sNMF, Frichot et al., 2014, imple-mented in the LEA package in r: Frichot & François, 2015), which is a likelihood-based approach to investigating ancestry proportions.
We iterated over K-values from 1 to 20, each with 10 replicate runs, and selected the K-value with the lowest cross-entropy score, which indicated lowest error in assigning number of supported ancestry (Frichot et al., 2014.). Finally, we visualized all individuals on a phylogenetic network generated by the neighbour-net algorithm (Bryant & Moulton, 2004, implemented in splitstree version 4.16.1 Huson & Bryant, 2006).
We estimated historical relationships among the seven valley populations using treemix (v1.12, Pickrell & Pritchard, 2012). We generated 1000 bootstrap replicates, each with a 1000-SNP resampling block. We first estimated the maximum likelihood phylogeny without allowing for migration or admixture, then sequentially added migration events to test whether any of these significantly improved the likelihood fit. The p-value for the added weight of migration edges is derived from jackknife estimates (but see discussion in Pickrell & Pritchard, 2012).

| Modelling historical gene flow
To understand the historical process that generated contemporary genetic structure of T. shambalaensis populations, we conducted  Table 2a). Current population distributions are represented as dots on a 840 × 600 raster. We extrapolated the historical distribution of each population from its current location: post-LGM L glacial retreat carved out the existing drainage terrain at Mt. Gongga, so we assumed that each surviving population has persisted in the periglacial refugia cut by the receding glacier valley river. We further assumed that these populations did not move across the flanking mountains. Figure 3a shows how this was simulated backward in time, where the locations of historical populations were represented as a few pixels on map, being deposited along the river valleys as the glacier "advanced". We first calculated six sets of Euclidean distances between populations (one on the modern landscape, as well as five historical landscapes, see

| Isolation by resistance
A potentially more accurate way to model the effect of environment on gene flow is to consider the environmental "resistance landscape" to gene flow (McRae et al., 2008). Since all climate variables that predict the habitat of caterpillar fungi are highly correlated with elevation (Hopping et al., 2018), we used change in elevation along the "dispersal paths" as a proxy for environmental resistance: (1) In a first set of distance calculations, "dispersal paths" were the 2D "least cost" path derived in the previous section ("Glaciation history"), but we integrated the total distances traversed as if these paths were projected onto a 3D plane with elevational data. Compared with 2D Euclidean distances between habitats, this new measurement imposes much higher costs for individuals directly crossing the high mountains separating neighbouring cul de sac habitats.
(2) In a second set of distances, we applied Dijkstra's algorithm (Dijkstra, 1959) in circuit theory to find "least resistance paths" between populations on a 3D plane.
We modelled the "resistance" to gene flow as the inverse of the absolute difference in elevation between two adjacent cells in our 840 × 600 raster (i.e., symmetrically higher resistance for making bigger "jumps" to higher or lower elevations), and calculated the path connecting populations with the least amount of total "fric- between 2D and 3D planes, these resistance distances (generated in 3D) were also projected down to 2D (called "least resistance path 2D" in Figure 3d). All path calculations were repeated on six (contemporary and historical) landscapes, with glaciers modelled as impermeable barriers to gene flow.

| Environmental contribution to divergence
The above-mentioned models predict paths (of "least cost" or "least resistance") by which the populations would disperse across a topologically complex landscape but do not quantify the effect of environmental variance on divergence. High mountains and low valleys are stronger barriers to gene flow than a flat landscape not only because of extended topological dispersal distance, but also because organisms are not well-adapted to the environmental range traversed (Ghalambor, 2006;Janzen, 1967). We thus used Bayesian estimation of differentiation in alleles by spatial structure and local ecology (bedassle, Bradburd et al., 2013) to quantify the contribution of geographic and environmental isolation on genome-wide differentiation between our populations.
Following the model formulation of Bradburd et al. (2013) that specified the effect size of both the geographic matrix (pair-wise Euclidean distance) and the environmental matrix (in our case, pair-wise elevational differences), labelled as α D , and α E , we con-  Note: Only the 15 most likely models are shown here. See Table S4 for full table. Model details can be found according to model numbers in the Supplementary material. Each model is specified by the "time periods" for which the population location is inferred (from modern to LGM L ), which habitat range is applied to allow for gene flow (climatic or elevational range), and whether or not glaciers are set as impermeable barriers to gene flow. Paths in each model were found by Euclidean distance, shortest path ("least cost") or calculating "least resistance" path in circuit theory. The established paths were traced and tallied either on 2D plane or 3D (accounting for distance moving up and down the elevation), and the BEDASSLE coefficient can be added to account for contribution of environment to genetic divergence. We calculated the Mantel R of each set of distances with pairwise population F ST . Each model was run independently 100 times in fastsimcoal2, and we show the mean and minimum AIC value of each model. All models have the same number of parameters (50).  Wang, Pan, et al. (2013). The glacier rivers (blue lines) are traced from current water flows. (b) Example of sets of relative distances between populations, between population location: Euclidean distances, "least cost path" that avoids crossing over glaciers, "least cost path" within suitable habitat range (mean temperature of the coldest quarter -10 to -4°C) labelled in grey and "path of least resistance" that minimizes total change in elevation between paths connecting two populations. (c) Mean fastsimcoal2 simulation maximum likelihood estimates for models based on different glaciation periods (all models have the same number of parameters). Bar indicates the mid-50 percentile of the 100 independent runs. Models higher on the y-axis are more likely given the observed data. (d) Mean fastsimcoal2 simulation maximum likelihood estimates for LGM models by distance type, representing a subset of "MIS3 max" shown in (c). Models are shown as three separate categories of paths: those that consider landscape as flat (light turquoise), as 3D (light blue) and as 3D with environmental contribution to genetic divergence (dark blue). Models higher on the y-axis are more likely given the observed data In these set of models, we followed the "least cost paths" and "least resistance paths" previously calculated on the 3D landscape, but we applied an extra distance of environmental change to account for accumulated change in elevation in the path (e.g., All sets of distance matrices generated are presented in Table 2b (also see Supplementary material). Prior to coalescent simulations, we also applied Mantel tests (Mantel, 1967) to check whether any of these new distance matrices differed significantly from the matrix of Euclidean distances between contemporary populations (that is, the one we used to assess IBD). We also applied the same tests for correlation between population pairwise F ST and these new distance matrices.
We used the above-mentioned sets of distance matrices as con-

| Modelling divergence
To estimate the timing of migration and population expansions, we is modelled as population coalescence with a resizing event. Coalescence, backward from present to historical time, is shown from bottom to top. We also posited a recent resizing event at T Exp for each population. In the pre-expansion (historical) migration model (horizontal arrows indicated in "Hyp 1"), populations were allowed to migrate between T Exp and T HX . In the post-expansion (contemporary) migration model (horizontal arrows indicated in "Hyp 2"), populations were allowed to migrate from present to T Exp population's one dimensional SFS. For each population we tested the likelihood of 4 models: (1) constant population size; (2) constant population size followed by constant growth or retraction; (3) constant population size followed by exponential growth or retraction; (4) constant population size followed by instantaneous decline followed by growth or retraction (bottleneck). Models (2) and (3)

| Single valley analysis
To test whether subpopulation genetic structure is distinguishable within a single valley population, we repeated our DAPC (Jombart et al., 2010), sNMF (Frichot et al., 2014) and neighbour-net (Bryant & Moulton, 2004) analysis on 30 individuals collected within the YZG valley, which we selected for more detailed substructure analysis because it includes an accessible main glacier and a river that melts from a northern glacier. Our samples belonged to 5 localities within the YZG valley on different sides of the glacier and rivers (Figure 5a, Table S1).
We could not construct an ML phylogeny at the subpopulation level due to a lack of sufficient genetic differentiation. Instead, we used fastsimcoal2 and modelled six different scenarios of subpopulation history involving glaciers changes (see Table S3). The observed joint SFS of all subpopulations were calculated using the number of projections that maximizes the number of segregating sites in each subpopulation. Since our aim was to select the best model for lineage divergence history, we did not allow cross-subpopulation migration in any of the models, and all parameters were specified the same as in 2.3 and 2.4. These settings would probably increase the estimate for coalescence time, but were unlikely to affect the lineage topology of the most likely model. For each model, we conducted 150 independent runs to determine the simulation with maximum composite likelihood (see Supplementary material for settings).

| Sequencing and SNP-calling
After quality-filtering, we obtained on average 6.5 million reads per sample (n = 122, SD = 2.1 million; sample information see Table S2).
A total of 23 samples (19%), mostly from GBG valley, shared only a negligible number of enzymatic cut sites with T. shambalaensis ( Figure S1) and did not form a monophyletic group with T. shambalaensis in the CO1 phylogeny (see Wang, Zhuang, et al., 2019). These belong to an undescribed Thitarodes species only distantly related to T. shambalaensis), and were not included in our analysis. After filtering the SNP catalogue generated by cstacks, selecting loci with less than 10% missing data across all individuals, we obtained a total of 11,038 SNPs across 2373 loci. One SNP per locus was randomly selected. Our final data set contained 2373 SNPs from 99 individuals, with 7.09% missing data (see Table 1 for post-filtering sample counts, see Figure S2 for STACKS parameter sensitivity, see Supplementary material for final SNP matrix).

| Population structure and phylogeny
With distance measured as Euclidean distances, population pairwise F ST showed a significant signal of IBD (Mantel p = .03, r = .52, Figure S3A). With environmental distances measured as population elevational differences, we did not detect a signal of IBE (Mantel p = .56, r = -.02, Figure S3B). Population statistics were summarized in Table 1. Although F IS per population was high (mean = .17, SD = 0.05) and H obs was lower than H exp (mean difference = 0.03), we did not detect any loci that significantly deviated from HWE (p = 1.00 after Bonferroni correction). The genetic diversity (H obs and π) and F IS of each population had no significant correlation with geographical coordinates, elevation or sample size (multiple regressions, p π = 0.39, p Hobs = 0.37, p F = 0.13, Figure S4).

| Modelling the landscape of historical gene flow
We generated 47 sets of distance matrices between populations to account for historical landscapes and environmental resistance to gene flow ( Table 2b). Most of these distance sets are minor modifications of, and thus highly correlated with, the contemporary Euclidean distances ( Figures S6 and S7). In reconstructing the historical position of populations from moraine records, the relative distance between YZG and MZG populations, as well as that between HG and GBG populations increased since MIS3, while the relative distance between XG and HLG have decreased (Figure 3a). The absolute distances across all populations increased from postulated LGM L localities to contemporary localities by 36% (SD = 0. 27, Figure S8). Our BEDASSLE simulations reached a "conversion coefficient" (α E /α D ) of 1.4 (i.e., each 1 m change in elevation of T. shambalaensis habitat amounts to 1.4 m of geographical isolation in terms of its contribution to genetic variation, Figure S9). The strengths of IBD correlations (measured as Mantel r statistics, McRae & Beier, 2007) in historical models were significantly greater than those in contemporary distance models (t = 5.17, df = 26, p < .01); the strengths of IBD correlations in models based on least resistance paths were significantly greater than those in models based on Euclidean distances (t = -2.52, df = 35, p = .02).
Four models failed to coalesce in our fastsimcoal2 simulations (Table S4, Figure S10); these are models in which gene flow was limited to an elevational range or climate range, while setting glaciers as impermeable barriers of dispersal completely blocked paths between populations. All models have the same number of parameters (50). Among the 44 models that coalesced, those that based their distances on MIS2, MIS3 and LGM L population distributions had significantly higher likelihood than their modern and Holocene counterparts (Figure 3c; t = −6.29, df = 41.41, p < .01). The most likely models were based on the population locations of LGM L (Table 2b).
Among models that used the population distribution of LGM L , those that were based on Euclidean distances had the lowest likelihoods, significantly outcompeted by ones with the more realistic consideration of environmental resistance, such as avoiding glaciers, or

| Modelling divergence
Between the pre and post-expansion migration models of monophyletic populations of YZG, HLG, XG and MZG, we found greater support that migration among the populations happened before population resizing (Table S5). The most likely model placed the start of population resizing at 41 k generations ago (95% CI: 34-42 k, Table 3a). Assuming a 1-1.5 year generation time for Thitarodes (Tao et al., 2016), this places the timing of likely among-population migrations during LGM L , which is consistent with our result that most likely models of historical and environmental resistance of gene flow are based on the population location of LGM L .
Our fastsimcoal2 models suggest that with the exception of XG (whose population remained almost constant), the post-LGM L population experienced a 30-to 100-fold expansion (Table 3a). Similarly, 1D SFS simulations of each population from δaδi suggested that, with the exception of HLG, which showed population decline, scenarios of "constant growth" and "exponential growth" (4-20 fold expansion) are most likely for all populations (Table 3b, model selection see Table S6).

| Single valley
Among 30 individuals in five different localities within YZG, both DAPC and sNMF analysis found greatest support for two genetic clusters (Figure 5b,c, Figure S11). Individuals adjoining the northern glacial river (from locality 4 and YLP) form one genetic cluster, while individuals on both sides of the main glacier (HZD south of the main glacier, and habitat 2, 3 north of the main glacier) form another genetic cluster. This clustering is consistent with the phylogenetic network of neighbour-net, which separated the two clusters by their glacier association (Figure 5d).
Among six fastsimcoal2 models of different subpopulation histories, we obtained greatest support for a lineage history in which subpopulations in YLP and location 4 are sister to the monophyletic cluster of HZD, location 3 and location 2 ( Figure S12, Table S3). The hypothesized tree topology was based on the scenario in which an ancestral population near the periglacial refugia diverged into a subpopulation at locality 4 and YLP as the norther glacier cut into the valley. Later, in smaller glacier cycles in which the YZG main glacier advanced, HZD was physically separated from locality 2 and 3 ( Figure 6).

| DISCUSS ION
Our LGM L . Since many of our historical inferences are model-based, we applied different methods of inference to verify that our results are robust even when subjected to specific model parameter choices.
By observing which results hold up across different models, we can be more confident about the validity of our historical inference. We discuss each of our findings and their implications below.

| Fine-scale population structure and IBD
Thitarodes shambalaensis populations in this study come from parallel glacier valleys no more than 10 km apart, but we were still able to select for 6-7 well-supported genetic clusters that correspond to each population's geographical distribution as our most likely hypothesis of population structure. In other words, variation exhibited by 2300 SNPs was sufficient to trace each individual to its habitat. While genetic structure of populations within 10 km has been reported in plants (Geng et al., 2009;Helsen et al., 2015;Hevroy et al., 2018;Kitamura et al., 2018), amphibians (Kobayashi et al., 2018), reptiles (Pearson et al., 2020), fish (Ciannelli et al., 2010) and mammals (Ziege et al., 2020), fine-scale genetic structure for flying insects is rare. For example, a similar number of SNPs used to study ants (Boyle et al., 2019;Smith et al., 2019), beetles (Weng et al., 2021) and butterflies (MacDonald et al., 2020) identified meaningful genetic clusters between long-distance populations that reflect their geographical isolation, but in these studies, populations within 10 km of each other were genetically hard to distinguish. Factors intrinsic to Thitarodes biology, as well as external geographical factors, contributed to the maintenance of genetic structure in geographically adjacent populations that are often characteristics of low-vagility, short-range endemic terrestrial arthropods (e.g., Peres et al., 2018;Schwentner & Giribet, 2018). For example, fungi-associated Thitarodes are known to have unusually strict habitat requirements (Hopping et al., 2018;Wang et al., 2020) and low adult dispersal due to their short period of activity (Maczey et al., 2010;Zou et al., 2011).  (van Strien et al., 2015). Central to the IBD process is the assumption that gene flow is proportionally inversely affected by the distance between populations. At a small geographical scale, the signal of geographical isolation is often overwhelmed by more immediate, sometimes anthropogenic, environmental differences such as a rural-to-urban gradient (Ziege et al., 2020) or cohort characteristics (Kitamura et al., 2018;Pearson et al., 2020). Signals of IBD at a fine scale are often weak (Crookes & Shaw, 2016;Hevroy et al., 2018), or otherwise nonexistent. Our TA B L E 3 (a) Inferred parameters and confidence interval for preresizing migration fastsimcoal2 model in Figure 3 (a) (2) constant growth or retraction; (3) exponential growth or decline; (4) instantaneous decline followed by growth or retraction (see Table S6 for model selection). Expansion coefficients >1 suggest population growth is either constant or exponentially increasing from RADseq are not well-suited for picking up signals of adaptation to environmental differences (Orsini et al., 2013). Furthermore, in modelling environmental resistance we assumed that upward and downward elevational deviations from current population ranges symmetrically impedes geneflow, while the minimum and maximum critical thermal tolerance of T. shambalaensis might be asymmetric.
To detect population-level adaptation of this Thitarodes system, additional research focused on assessing molecular markers under selection, incorporating data from physiology measurements, and taking into account differences of each population's strain of parasitic fungi as a component of "environmental" difference would be necessary.

| Glaciation and lineage history
While each population was isolated by its dispersal limits, our results also revealed a signature of glaciation history on population divergence. Specifically, our ML phylogeny found samples from the southern valleys of YZG, MZG, XG and HLG to be monophyletic.
Since these valleys were carved out from the retreat of a single large ice sheet post-LGM L (and still share the same drainage system west of the Moxi river, Figure 1b), the monophyly of the four lineages suggest that a single ancestral population diverged into four populations as the receding ice-sheet carved out four parallel valleys. This is not the case for samples from the northern valleys (HG, GBG, NMGG).
The formation of these northern valleys each resulted from an individual ice sheet retreat, and we found no signal of monophyly between populations in the northern valleys. We note that one way to provide further evidence that ancestral periglacial species diverged into the four southern valleys is to look for phylogeographic concordance among other species co-distributed with T. shambalaensis, which would suggest a convergence in response to topographic history (Marske et al., 2020;.

| Gene flow in historical landscape
Our reconstruction of fine-scale spatial-temporal environmental resistance to gene flow between populations addressed three major considerations in landscape genomics: (1) First, we constructed spatially explicit models of populations across six time periods. We reconstructed models of population distribution from moraine records. This approach allowed us to reconstruct historical scenarios as temporally-explicit as the moraine record allowed (five periods in the past 50 kya). These models assume that ancestral populations arrived at their present location in each valley following the receding ice sheets (instead of in situ "nunatak" preservation). If populations had always remained in their current positions, we would find the IBD pattern best explained by contemporary distances between populations.
(2) Second, by implementing different criteria for selecting paths along which migration occurs (Euclidean, "avoid glaciers", "least cost", "least resistance"), we contribute to the on-going discussion regarding the best practice to quantify landscape permeability to gene flow (Fant et al., 2014;MacDonald et al., 2020;McRae & Beier, 2007).
(3) Third, we applied Bayesian methods to quantify the contribution of environment and geography to genomic divergence (Bradburd et al., 2013), and converted the results into a single measure of distance between populations.
Comparing sets of pair-wise distances generated from a combination of these three considerations allows us to pinpoint historically significant factors affecting T. shambalaensis genetic structure.
Instead of incorporating results of Mantel/partial Mantel tests in a causal modelling framework to determine the set of distances that best contribute to genetic variation, we specified these distances as parameters in fastsimcoal2 coalescence to more accurately simulate the "spatial pedigree" between each individual (Bradburd & Ralph, 2019). Compared with coalescent simulations that impose presence/ absence of migration, our approach of adjusting the relative number of migrants between populations has the advantage of keeping the number of parameters constant across comparisons. However, a drawback to our approach is the need for sufficient computational resources required to enable coalescent models to converge, if at all, given as many as 50 parameters but only thousands of SNPs.
Future studies might compare our approach with models under the integrated distributional, demographic and coalescent framework (iDDC; He et al., 2013). While we mentioned in the introduction that lack of high-resolution environmental niche modeling (ENM) data might affect the application of such a framework on fine-scale landscapes, a work-around under iDDC might involve using the extent of historical ice sheet or elevational data to directly condition key demographic parameters in the model --thus circumventing the need of ENM. While such an approach would remain computationally demanding and the use of approximate Bayesian computation (ABC) model selection would be challenging, it has the advantage of reducing the number of model parameters.
Along the spatial axis, models using gene flow paths with consideration of environmental resistance and 3D topology outperform models of direct Euclidean paths (Figure 3d). This again highlights the effect of mountainous topology in redefining "distance" among populations. Along the temporal axis, our coalescent models show that contemporary T. shambalaensis genetic structure is best explained by gene flow among LGM L populations (Figure 3c). That is to say, a signature of gene flow among populations in the LGM L refugium is still preserved today. This is possible because post LGM L , populations expanded deeper into individual valleys following retreating ice sheets, but during this process, migration between populations was greatly reduced. Our study provides multiple lines of evidence for this "migrate, then expand" scenario: (1) First, while relative distances between populations fluctuated from past to present, the absolute distance between each population (measured by Euclidean distance, "least cost" distance, "least resistance" distance) increased from past to present ( Figure S8). In other words, populations have become more isolated from each other by migrating deeper into the valleys.
(2) Second, although migration could be continuously happening at varying rates throughout population history, our direct comparison between pre-resizing migration and post-resizing migration simulations (using the same number of parameters) showed higher support for the former models (Table S5). The timing of migration in the best "preresizing migration" models also corresponds to LGM L (Table 3a). In terms of resizing dynamics, both fastsimcoal2 and δaδi simulations indicate population growth as the best model for individual populations (Table 3). (3) Lastly, contemporary population statistics did not detect deviation from HWE; present-day migration, if any, was not sufficient to be detected statistically.
All these lines of evidence suggest that post-LGM L T. shambalaensis divergence along glacier valleys gradually increased the degree of isolation between populations. If it has not happened already, we suspect that given sufficient time and continued isolation,

Dobzhansky-Muller incompatibilities (intrinsic reproductive barriers)
will develop among these populations (Orr & Turelli, 2001). Similar mechanisms could explain the contribution of glacier refugia to the generation of insect diversity among the Hengduan Mountains.

ACK N OWLED G EM ENTS
ZW was supported by a graduate fellowship and a Royall Tyler Moore

AUTH O R CO NTR I B UTI O N S
Zhengyang Wang and Naomi E. Pierce conceived of the study.
Zhengyang Wang conducted the study, carried out the analysis and wrote the original draft. Both authors prepared the final manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw reads, scripts for running analysis, as well as model parame-