3 Results

3.1 Population-level genomic diversity

In the Himalayas, the high-elevation H. tibetana showed genome-wide heterozygosity between 0.03 to 0.04 (Fig. 2a) with populations in the Koshi Basin exhibiting higher heterozygosity compared with those from the Gandaki Basin. We generally observed increasing heterozygosity from the western to the eastern Himalayan populations, while the F decreased along the same geography. Populations of the low-elevation H. digitata exhibited a comparatively lower level of genome-wide heterozygosity and F (π = ~0.01, F = ~0.13). Moreover, there was no significant difference among populations from different basins.
Both HM species, H. gregoryi (high-elevation) and H. platon (low-elevation) showed similarly low levels of heterozygosity (π = ~0.001). However, populations from some of the HM basins exhibited an extremely high F while the others had a lowF . For instance, the F of most populations of H. platon were around 0.13, but the populations from the three basins located in the remote area of the Yangtze basin (ID 1634, 1280, 1205, Fig. 3) showed F approaching 1.0 (Fig. 2h). Moreover, the heterozygosity of these populations was relatively lower than that of the remaining populations. Similar to H. platon , populations ofH. gregoryi also showed fluctuating levels of F across basins: populations from some of the subbasins of Mekong (ID 2531, 2532, 2544) had a high level of F (0.5–0.7), as did populations from one subbasin of the Yangtze (ID 1684). In particular, these basins were geographically isolated, at least at the level of subbasins.

3.2 Population structure

We found four distinct clusters in the PCA plot for the high-elevationH. tibetana (Fig. 3a, b, c): (i) West and Middle Gandaki (subbasin ID 4900 and 4754); (ii) East Gandaki (ID 4789) and West Koshi (ID 4547); (iii) Middle Koshi (ID 4427); (iv) East Koshi (ID 4442). The MJ network of mtCOI sequence data, also supported four distinct clades, except that several individuals of cluster i shared mtDNA haplotypes with cluster iv, which may indicate a preservation of these old haplotypes at Middle Gandaki (ID 4754) or an invasion of individuals from East Koshi (ID 4442, Fig. S22, Fig. 3c). This same result was supported by the maximum likelihood estimation of individual ancestries using NGSadmix with the best K (K = 3), and with the taxonomic units delineated by TreeMix (Fig. S23). Analyses of PCA, admixture (with the best K = 4) and TreeMix revealed congruent structuring patterns for populations of H. digitata : populations were grouped into 5 clusters representing populations from (i) West Gandaki (ID 4900); (ii) Middle Gandaki (ID 4754); (iii) East Gandaki (ID 4789); (iv) Middle Koshi (ID 4390); (v) East Koshi (ID 4407 and 4442, Fig. 3d, e, f, Fig. S23). In addition, the MJ network of the mtCOI sequence data showed evidence of secondary contact between populations from cluster i and populations from other groups (Fig. S22).
When compared with the two species distributed in the Himalayan region, populations of H. gregoryi and H. platon showed more complex structures. Three distinct groups were recognized within the samples of H. gregoryi from (i) Mekong (including subbasin ID 2531, 2532 and 2544); (ii) West Yangtze (ID 1684); (iii) East Mekong and West Yangtze (ID 2547, 1668 and 1685, Fig. 3g, h, i). Although the optimal K could not be determined following the methods outlined in Evanno et al. (2005), we found congruence between the output of admixture and PCA: the samples from group iii formed a cluster in the PCA that simultaneously exhibits a gradient among populations of different elevations and subbasins. This cluster was assigned to an admixed group and one clade in the TreeMix but split into multiple lineages (Fig. S23). It was notable that samples within group iii were located in the border region of the three subbasins located in the same mountain area. The population structure of H. gregoryi estimated by the genomic data was partially congruent with the results of the MJ network: samples of group iii were assigned to the same haplogroup, yet, samples of groups i and ii were segregated into four haplogroups (Fig. 3i, Fig. S22). The 12 populations of low-elevation species H. platon were grouped into four distinct clusters: (i) Middle and East Yangtze (ID 1634, 1280 and 1205); (ii) West Yangtze (ID 1625); (iii) Salween (one population from subbasin ID 2949); (iv) populations from remaining basins (including Irrawaddy, Salween, Mekong and Yangtze, Fig. 3j, k, l). Similar patterns were also observed in the mtCOI MJ network, except that one individual from West Yangtze (ID 1625) fused with the near-panmictic haplogroup (iv), indicating a secondary contact (Fig. S22).
Interestingly, the PCA using thein situ and GIS environmental variables successfully separatedH. gregoryi and H. platon in the HM (Fig. S6 and S7), which indicated a different microhabitat of the two species. However,H. tibetana and H. digitata in the Himalayas were not sufficiently separated from each other in the PCA, suggesting that the two species might inhabit a similar microhabitat.

3.3 Demographic history and SDM

The demographic history of each species was investigated by PSMC analysis of each genome and as represented by the inverse instantaneous coalescent rate (IICR; Fig. 4 c, f, i, l). Although the IICR estimates became less reliable at more recent time scales, which led to a great variance of the effective population size (Ne ) between individuals, the results revealed a general trend of an expansion-decline event during or by the end of the LGM (22 ka) for all of the four target species. In addition, the Ne of low-elevation species were 10,000–100,000 times greater than those of the high-elevation distributed species. When comparing the mountain ranges, the Ne of the HM and the Himalayan species were more or less at a similar level at their respective peaks.
Estimation of the effects of climatic versus topographic variables using ensemble models produced clear results to explain the potential distribution of species. The highest percentage contribution of the ensemble model was the seasonal precipitation for the two species in the Himalayas, and the seasonal temperature, especially the annual range of temperature for the two species in the HM (Table S5). In addition, for the ensemble modelling of all four target species, the ROC and TSS scores are very high (>0.99, Table S4), indicating a very high model performance.
Coinciding with the maximum Ne during the LGM and subsequent population decline, the predicted climatically suitable area of all these four target species was greatest during the LGM and contracted to the present day (Fig. 4 a, d, g, j). For instance, during the LGM, the predicted suitable area of H. tibetana and H. digitatacomprised most of the Himalayas, north of Khasi Hills and part of the HM. However, the present-day suitable area shrunk to a small region in the middle of the Himalayas, which aligns with the known present-day distribution range. Interestingly, except for the middle of the Himalayas, the modelling identified another area on the border of the Khasi Hills and HM that has been persistent for the existence of the two species since the LGM. However, no actual presence of the two species (Fig. 4 a, d) was observed there. Similar to the species distributed in the Himalayas, species distributed in HM also experienced a contraction of climatically suitable area from the LGM to the present, especially for H. gregoryi . The suitable area of H. gregoryi during the LGM spread all over the southern and eastern part of the HM or further, as well as the eastern Himalayas, but dramatically shrunk to the central part of HM, which still vastly extends the known present-day distribution range. In addition, the suitable areas for H. gregoryi and H. platon shifted northward over time. Furthermore, the predicted elevational range for all four species during the LGM is lower than at present, implying a tendency for species to move uphill in the Holocene (Fig. 4 b, e, h, k).