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).