4 Discussion

To further understand the origin and evolution of mountain biodiversity, we investigated the genetic diversity and distribution of a caddisfly group in the Himalayas and the HM. For the first time in aquatic insects, our results are derived from whole genome re-sequencing of populations of four congeneric species and fine-grain environmental GIS data. And, while we interpreted our results carefully, our study provides evidence that changes in geodiversity and climate affect the evolutionary history of closely related species occupying varying ecological niches, thus supporting the MGH.

4.1 Genetic patterns of caddisfly in the Himalayas and the HM

Present-day genetic diversity patterns offer important insights on the historic interconnectivity of populations and can help inform the effect of past environmental change on the studied species. Our results indicate that in most of the cases, the high-elevation species of both mountain regions (H. tibetana and H. gregoryi ) align best with the headwater model pattern sensu Hughes et al. (2007). For instance, populations of H. tibetana (high-elevation species, Himalayas) in the West Gandaki (subbasin ID 4900) and Middle Gandaki (ID 4754) cluster together and share the same homogeneous ancestry (Fig. 3a, b, c). Instead of being connected by a stream network, these populations are connected by a series of local mountain ranges (including Tilitso Himal, Anna Dakshin, Machhapuchhara and Lamjung Kailas), as defined in the headwater model pattern. Similarly in H. gregoryi (HM): populations in the Baima Snow Mountains (occupied in North Mekong (ID 2547) and West Yangtze (ID 1668 and 1685)) originated from an admixed ancestral group and formed one haplogroup based on mtCOI (Fig. S22, Fig. 3g, h, i), in addition to a low level of inbreeding coefficient (F = ~ 0.12), indicating a good connection within this mountain range. Outside the Baima Snow Mountains, the H. gregoryi populations located in the North and South Mekong River (containing subbasin ID 2531, 2532, and 2544) were associated with a stream hierarchy model since they show non-significant genetic variance among each other. Nevertheless, the ancestral group of these populations is distinct from others, suggesting limited gene flow between this group and populations from other mountains. Therefore, the populations of high-elevation species are generally restricted by local mountain ranges and show a stronger “insularity” among the different mountain ranges, or a “Sky Island” pattern. This pattern implies some degree of overland dispersal among headwater streams. Caddisflies can disperse overland during their adult stage, although the associated dispersal distance seems to be restricted (Geismar et al., 2015). Nonetheless, high-elevation species likely disperse along local mountain ranges within or among basins, using suitable habitats (so-called “montane archipelagos”) as stepping stones, thus leading to better connectivity within these mountain ranges. Moreover, due to their strong association with specific habitats, lowland areas, especially deep and large valleys, may act as dispersal barriers for these species, thus leading to limited gene flow among mountains (Myers et al. 2001; Pauls et al. 2006). In fact, freshwater organisms, in particular headwater specialists, are often endemic to specific mountain ranges due to their low tolerance to elevated temperature or other habitat factors characterizing the surrounding lowland (Finn et al., 2007; Hughes et al., 2009; Bálint et al., 2011; Vitecek et al., 2015, 2017).
In comparison, for the low-elevation species, a significant genetic diversity was shared among populations across basins, thus providing support for the stream hierarchy model pattern. It suggests gene flow mostly occurred among populations from the same basins, indicating that dispersal within the stream network appears to be the dominant process. The divergence among dispersal strategies of high-elevation and low-elevation species is evidenced by the pattern we observed where populations of H. digitata and H. tibetana in the West and Middle Gandaki (ID 4900, 4754): although populations of these two species partly overlap, populations of the H. digitata appeared as two distinct groups associated with the two subbasins, while populations of H. tibetana grouped together (Fig. 3c, f). In addition, the population connectivity of H. digitata reflects the dendritic nature of the stream as described by Hughes et al. (2013): minimal differentiation was observed for populations within the West Gandaki subbasin (e.g., individuals from subbasin ID 4900, Fig. 3d, e, f), higher differentiation among populations from different subbasins (e.g., ID 4900 and 4754) and the highest among populations from different primary basins (e.g., Gandaki and Koshi). The stream hierarchy model pattern can also be observed in H. platon in the HM. Several populations located in the Yangtze River (subbasin ID 1625, 1634, 1280, 1205) formed distinct groups, suggesting higher connectivity among populations within subbasins (Fig. 3j, k, l). In addition, these populations exhibit low genomic heterozygosity (π < 0.001) and an extreme inbreeding coefficient (F = ~0.75 – ~0.1, Fig. 2d, h), indicating isolation and future allopatric speciation. The rest of the populations of H. platonthat extend across all four primary basins (Irrawaddy, Salween, Mekong, and Yangtze) demonstrate panmixia, accompanied by low genetic diversity and low inbreeding coefficient (Fig. 2d, h). This pattern suggests that all of these populations are highly connected and gene flow is strong within all four primary basins until recent times, which may be due to their very well-developed dispersal abilities as adults or continuous habitats (Finn et al., 2007; Baggiano et al., 2011).
Our results highlight that identifying the prevalent model of genetic structure for aquatic organisms may not always be clear-cut, especially for taxa with a broad distribution (Hughes et al., 2013; Sproul et al., 2014). In our case some species even appear to follow multiple patterns: Stream hierarchy model patterns observed in the low-elevation species could also be interpreted as isolation by distance since the genetic diversity also increased with geographic distance (or “stream distance”, the same as geographic distance in this case).
In line with previous studies that suggested that a general pattern of population structure is challenging to identify for caddisflies (Lehrian et al., 2009; Pauls et al., 2009; Engelhardt et al., 2011; Previšić et al., 2014), our study also reveals that the patterns inHimalopsyche are likely to be species-specific: high-elevation species are associated with a headwater model pattern, whereas the low-elevation species tend to be more similar to stream hierarchy model or panmixia patterns, indicating that habitat specialization is a crucial factor for shaping the population structure. Further, species in the HM show a better-connected gene flow within the mountain range than species distributed in the Himalayas (panmixia appears in the HM), implying that the effects of topography are strong for the dispersal of species with an aquatic life history in this region. Overall, these complex patterns may result from either the complexity of the dendritic stream network embedded in a high-gradient landscape or simply the geographic scale of the studies examined (Hughes et al., 2013).

4.2 Key abiotic factors in shaping the montane biodiversity of caddisflies

Ecological preferences and dispersal abilities cannot explain present-day distribution patterns without interpreting the data in the context of a species’ history (Avise 2000). In particular, topography and historical climate changes may have a strong impact on the current genetic diversity pattern of caddisflies. Our study, combined with the current literature (e.g., Favre et al., 2015; Antonelli et al., 2018; Mosbrugger et al., 2018; Muellner‐Riehl et al., 2019; Muellner-Riehl, 2019), suggests that the evolution of climate, mountain topography, as well as historical drainage rearrangement, are key factors for the evolution of aquatic biodiversity in the Tibeto-Himalayan Region. These factors have interacted with each other through time and vastly contributed to shaping current diversity patterns in these aquatic insects.

4.2.1 Genetic patterns versus climate and topographical features

The influence of climate on population diversity of caddisflies is significant, especially in the Himalayan region. In line with other organisms in previous studies (e.g., Wikramanayake et al., 2002; Price et al., 2011; Yan et al., 2013; Rana et al., 2019), we find an increased genetic diversity of our target caddisfly species from west to east in the Himalayan region at the population level, especially for the high-elevation species. There may be a link between this climatic gradient and the genetic diversity of caddisflies across the whole mountain region, leading to patterns associated with topography and longitude (Bookhagen & Burbank, 2006, 2010): species richness increased approximately threefold from the northwest to the east of the Himalayas, especially at higher elevations, which is in accordance with the fact that the lowland receives about three times more precipitation in the subtropical East than the Northwest (Anders et al., 2006; Kreft & Jetz, 2007; Bookhagen & Burbank, 2010; Rana et al. 2019, 2022). This pattern is strongly supported by studies on regional plant diversity (Yan et al., 2013; Muellner-Riehl, 2019; Rana et al., 2022). Precipitation may be the direct or indirect factor resulting in the current distribution of caddisfly species (Collier & Smith, 1997). Furthermore, since both larval and adult stages of caddisflies are reliant on flowing waters, surface water availability from rainfall and snowmelt in the east of the Himalayas may significantly benefit them (Bookhagen & Burbank, 2010).
Compared to the straightforward genetic diversity pattern of caddisflies in the Himalayas, a common distribution pattern of genetic diversity is difficult to identify for caddisflies in the HM. The complex terrain characterizing the HM region offers highly heterogeneous microclimates, making it difficult to evaluate the impacts of climate change in general (Yin et al., 2020). Yet, a variety of genetic diversity patterns of caddisflies were observed in the HM, which is strongly correlated with topographic complexity, in particular those associated with elevation and hydromorphology. Generally, species in the HM show better connectivity all over the region, but there are plenty of allopatric populations observed for both species. For the high-elevation species, the main distribution pattern is the “Sky Island”, thus the geographic barrier made by certain lower elevations. However, for the low-elevation species, the mechanism of formation of the allopatric populations is different. Unlike the alpine species, the low-elevation species generally show better connectivity within basins, especially in H. platon which shows near-panmixia across four primary basins. Furthermore, the fact that both strongly connected (the upper Yangtze: Jinsha River) and allopatric populations (the middle Yangtze: Yalong River and Dadu River) occurred in the Yangtze Basin, may derive from ancient drainage re-organizations (e.g., via a so-called “river capture”, Fig. 3m). Firstly, the upper Yangtze River used to flow southward as a tributary of the paleo-Red River until the late Eocene (Clark et al., 2004; Zheng, 2015; He et al., 2021; Zheng et al., 2021), which was not the case of its middle part (see Zhang et al., 2017). Fuelled by regional orogeny and erosion (Zheng, 2015; Wang et al., 2018), a so-called river capture event occurred, changing the course of the upper Yangtze away from the paleo-Red River and into its current course (Clark et al., 2004; Zheng, 2015; Deng et al., 2021a). Our results reveal patterns consistent with an ancient disconnection: populations located in the Yalong (ID 1634) and Dadu River (ID 1280 and 1205) showed a distinct genetic ancestry with any other populations with any K values, including those located at the upper Yangtze River (ID 1625 and 1659). Secondly, although the populations located in the Yanglong and Dadu River represented the same ancestral, they are distinct clusters in population structure, indicating a recent divergence. It is in line with the interpretation that the paleo-Yalong and paleo-Dadu Rivers drained together with the middle Yangtze to the paleo Xigeda lake before ca. 1.3 Ma, but segregate after ca. 1.3 Ma since the paleo-Dadu River was integrated into the lower Yangtze directly (Deng et al., 2021a). Until ca. 1.3 Ma, the middle Yangtze captured the upper Yangtze and drained into the lower Yangtze drainage basin (Kong et al., 2009; Deng et al., 2021a). Thus, populations ofH. platon potentially experience some gene flow along the connected river networks, as shown in our TreeMix analysis: gene flow has occurred among populations from the upper Yangtze (ID 1625) and middle Yangtze (ID 1634). However, the genetic connection seems to be very limited, since the populations located in the middle Yangtze showed extremely high inbreeding coefficient.
Interestingly, the population of H. platon located in subbasin ID 1625 shares the same ancestral source with the near-panmictic populations when K = 3. However, when increasing the value of K, it is assigned to a unique ancestry and the near-panmictic populations are inferred to have a certain level of ancestry from this unique population with relatively high inbreeding (Fig. 3, Fig. 2). Notably, this population is located at the middle of the First Bend of the Yangtze River (as shown in Fig. 3m). The First Bend is the place where upper Yangtze makes a dramatic turn to the North and starts flowing eastward to the East China Sea (Zhang et al., 2017). It is regarded as the capture point of the upper Yangtze away from the Red River, thus resulting in a key period for geo(hydro)morphology in this area (Clark et al., 2004; Li et al., 2022). The processes leading to the present-day structure and diversity patterns remain unclear, however. Very limited information is available on the genetic diversity patterns of plants or animals here. It would be worthwhile to investigate species of this specific region in future studies.

4.2.2 Distributional and demographic history of the four caddisfly species

Repeated climatic fluctuations in the Pleistocene promoted range contractions and expansions of numerous plant and animal species, which is believed to result in high diversification rates in the Tibeto-Himalayan Region (Xing & Ree, 2017; Mosbrugger et al., 2018; Ding et al., 2020). However, the effects of Pleistocene glacial cycles on the distribution and demography of species may vary across regions (Hewitt, 2000). For example, in addition to ecological preferences, life-history traits, and dispersal ability (Tonzo & Ortego, 2021), local topography and latitude likely have strong effects on rates of dispersal, colonization, and extinction of a species during climate fluctuations.
Our ensemble model results show that in comparison to the present, suitable habitats for all four target species were more extensive during the LGM but shifted from their current core distribution (Fig. 4). Nevertheless, SDMs depict similar suitable habitats for species from the same mountain region, but divergent patterns for species from the Himalayas and the HM. In the Himalayas, the core potential habitat of both species (at intermediate elevations), seems to have only shifted a little between the LGM and today, indicating minor in-situvertical displacement along the elevational gradient (Fig. 4). Dispersal between the East and the West of the range seems also to have been possible, as indicated by the results of the ancestral components. However, suitable habitats at the eastern and westernmost part of the Himalayan range were likely only available during the LGM, indicating potential local extinction or long-distance migration if species occurred there before the LGM. In comparison, for both species in the HM, a latitudinal displacement (from South at warm times to North at colder times) and a vertical displacement (move up to a higher elevation) occurred after the LGM, which was not the case in the Himalayas. Latitudinal displacement was aided by the mountain ranges long North-South extent (Fig. 5b), whereas the Himalayas are latitudinally strongly constrained by the high plateau to the North and the tropical lowland to the South (Fig. 5a).
In line with the SDM results, evidence of the population size history based on genomic data showed that all species experienced an abrupt demographic expansion during the LGM, followed by a population size decline from the end of the LGM until recent times (Fig. 4). Previous studies combining genetic data and ecological niche modelling have also detected such demographic expansions during the LGM, particularly in cold-tolerant/adapted species, for example in grasshoppers (Tonzo & Ortego, 2021), Vipers (Yousefi et al., 2015), caribou (Taylor et al., 2021), birds (Garg et al., 2020), and plants (Manthey et al., 2017). TheHimalopsyche species we investigated are cold-tolerant, and their ability to disperse overland (to some extent) may have allowed them to persist in situ regionally throughout cold snaps. In fact, this was the case for many mountain caddisfly species in Central Europe (Pauls et al., 2006; Kubow et al., 2010; Lehrian et al., 2010; Previšić et al., 2014), despite a much stronger impact of glaciation on this continent (Višnjević et al., 2020) than in the Tibeto-Himalayan Region (Mao et al., 2021). The limited formation of inland ice during the LGM in the Tibeto-Himalayan Region may have represented an opportunity for cold-adapted species to expand their range rather than being moved far away by massive ice sheets.

4.3 Implication of MGH in different mountain systems

In line with the MGH, we have shown that the historical evolution of abiotic conditions such as changing topography (including paleo-drainage reorganization) and climate are likely to have played dominant roles in shaping the genetic diversity pattern of caddisflies in the Himalayas and the HM. The three MGH boundary conditions (i.e. extensive elevational zonation, response to climatic fluctuations, and a high‐relief terrain with environmental gradients; Mosbrugger et al., 2018) are met in both the Himalayas and the HM. Yet, we observe distinct genetic diversity patterns of species in each mountainous region, indicating that although the conditions of the MGH are generally realized opportunities and constraints for species vary across different mountain systems. To explain the formation of different genetic diversity patterns in each mountain region, it seems crucial to compare the geological and climatic conditions in each region and their effects on the evolution process of species under the concept of MGH.
The uplift of the Himalayas and HM are both caused by the India-Eurasia collision since the Paleocene (65–60 Ma). The ranges reached their modern height before the Pliocene (as reviewed by Ding et al., 2022). Although Muellner-Riehl (2019) indicates that the major uplift of the Qinghai-Tibet plateau proper would not have had an impact at the level of intraspecific divergence (as it is vastly older), the distinctive topographic features of these two mountain systems that resulted from this orogenic process has influenced the interconnectivity of habitats among and within each mountainous region (Craw et al., 2016; Rahbek et al., 2019). Compared to the Himalayas, the HM provide a more extensive area in the North-South orientation, thus implying better geographic connectivity for species to shift their distribution towards the south during the cooling period (Xing & Ree, 2017). Our results have verified that caddisfly populations in the HM show better connectivity and a broader potential habitat during the LGM. For this particular scenario that is associated with Pleistocene climatic fluctuations, Flantua et al. (2018) recommended a concept called “flickering connectivity system (FCS)” to explain the rapid diversification and high species richness in mountain ecosystems. Unlike the MGH, which attempts to explain the origination and evolution of mountain biodiversity throughout the whole uplift history, the FCS emphasizes the connectivity dynamics resulting from repetitive climatic changes during the Pleistocene. As reviewed by (Muellner-Riehl, 2019), many previous phylogeographic studies have indicated that the late Pleistocene is a crucial time period for intraspecific diversification in the Tibeto-Himalayan Region. Therefore, the FCS concept is applicable to explain the formation of intraspecific diversification of caddisflies that occurred in a more recent time period (the LGM).
Caddisflies of the HM and the Himalayas may have reacted to the climatic fluctuations of the Pleistocene by shifting their geographic range within mountain areas. Consequently, the dynamically changing connectivity levels may lead to fragmentation, colonization (dispersal), intermixing, and hybridization, thus promoting diversification in mountain ecosystems (Flantua et al., 2018; Muellner-Riehl, 2019). Moreover, since topography and climate conditions are distinct in different mountain systems, the diversification of species shaped by these four processes under the concept of the FCS is unique in each mountain. The resulting interaction between topography and climate is considered the “fingerprint” of a mountain (Flantua et al., 2018). Although the “fingerprints” are not yet specified for all the mountain systems of the world, our study revealed that the “fingerprints” differ between the Himalayas and the HM. Hence, in agreement with Muellner-Riehl (2019), we think that together with climate fluctuations as one driver of diversification in mountain systems sensu MGH, the FCS is relevant for explaining the montane biodiversity patterns of caddisflies. However, given the fact that species react differently to changes in their environment, more studies are needed in the future to explore how the entire biota reacts when facing climatic challenges.