Linking watershed formation with the phylogenetic distribution of a soil microscopic fungus in Yunnan Province, China

Abstract

Background

Phylogeographic studies have gained prominence in linking past geological events to the distribution patterns of biodiversity, primarily in mountainous regions. However, such studies often focus on plant taxa, neglecting the intricate biogeographical patterns of microbes, particularly soil microbial communities. This article explores the spatial distribution of the nematode-trapping fungus Arthrobotrys oligospora, a widespread microorganism, in a tectonically active region at the southeastern edge of the Qinghai-Tibetan Plateau. By analysing the genetic variation of this fungus alongside the historical structure of major river watersheds, we sought to uncover potential connections between the two. Our study involved sampling 149 strains from 116 sites across six major watersheds in the region.

Results

The resulting haplotype network revealed five distinct clusters, each corresponding closely to a specific watershed. These clusters exhibited high haplotype diversity and low nucleotide diversity, supporting the notion of watershed-based segregation. Further analysis of haplotypes shared across watersheds provided evidence for three proposed past river connections. In particular, we found numerous shared haplotypes between the Yangtze and Mekong basins, as well as between the Yangtze and the Red basins. Evidence for a Irrawaddy-Salween-Red and a Yangtze-Pearl-Red river connections were also portrayed in our mapping exercise.

Conclusions

These findings emphasize the crucial role of historical geomorphological events in shaping the biogeography of microbial biodiversity, alongside contemporary biotic and abiotic factors. Watershed perimeters emerged as effective predictors of such patterns, suggesting their suitability as analytical units for regional-scale studies. Our study also demonstrates the potential of microorganisms and phylogeographic approaches to complement traditional geological analyses, providing a more comprehensive understanding of past landscape structure and its evolution.

Background

Mountains host multiple, delicate, and unique ecosystems supporting rich biodiversity and high rates of endemism [1, 2]. The geographic isolation imposed by natural barriers such as high peaks, long ridges, and large rivers acts in favour of speciation and endemicity, limiting opportunities for the dispersion of living organisms [3, 4]. Among the major mountain ranges of the world, the distribution of fauna and flora at the eastern edge of the greater Himalayas region, a prominent example of temperate mountain biodiversity hotspot [5,6,7], has been investigated using phylogeographic approaches. Notably, significant associations between the genetic structure of current species populations and the drastic changes of the major watersheds found in the Central Yunnan Plateau have been found for several plant taxa [8,9,10,11,12,13], as well as for some frogs [14] and fishes [15,16,17]. Many of these studies refer to the reconstruction of the paleo-drainage system in southeastern Tibet proposed by Clark et al. 2004 [18], suggesting that previous to the Miocene uplift, all the upper reaches of the major rivers between the Tsangpo-Brahmaputra (west) to the Middle Yangtze (east) were connected in a dendritic pattern flowing southwards through the paleo-Red river into the South China Sea. Several river capture and reversal events occurred through the millennia, shaping the different watersheds to the current configuration [18,19,20,21]. More precisely, the authors propose that the current Yarlung-Tsangpo and upper Salween as well as the upper section of the Jinsha river (upper Yangtze) flowed into the upper Mekong river, which in turn flowed into the Red river. Moreover, the middle Yangtze was composed by several tributaries of the paleo-Red river, implying that the direction of the flow was reversed in some sections (Fig. 1a). However, there is still an open debate regarding the evolution of these rivers, the causes and timing of major drainage changes, as well as whether certain river connections ever existed [22, 23]. Some authors support Clark et al. 2004 hypothesis portraying a large paleo-Red river with sources from the current eastern Tibet and Sichuan regions but with differences in the inclusion of certain tributaries [20, 24, 25], while others refute this hypothesis and suggest a paleo-Mekong or an unnamed paleo-river with sources from eastern Tibet and including fluvial units from the Sichuan Basin flowing south-westward though the Simao, Vientiane, and Khorat basins [26,27,28,29,30]. Models including both paleo-Red and paleo-Mekong [31] or even a paleo-Salween rivers [32] sharing different sections of the current Yangtze river at different geological times, have also been put forward.

Phylogeographic approaches in this area have been almost exclusively applied to plants and a few animal taxa, while microorganisms, which represent the majority of the Earth’s biodiversity, still remain poorly explored. With the present article, we would like to contribute to the understanding of the complex links between biodiversity and geological events in the central Yunnan plateau by looking at the spatial distribution of Arthrobotrys oligospora, a widespread and abundant nematode-trapping fungus (NTF) found in the soil of both terrestrial and aquatic environments. While fungal carnivorism diverged from saprophytism about 419 Mya, it is estimated that A. oligospora evolved from its ancestor about 89 Mya [33], making it ideal for studies related to geological history. In a previous investigation, Deng et al. 2023 [34] collected soil samples from 228 sites across Yunnan province of China, identified A. oligospora strains and classified them at the phylogenetic level. By mapping the five resulting evolutionary clades (Fig. 1b), they found that the geographic location within major watersheds could predict 68% of the clades, while major climatic factors as well as Euclidean distance between sample locations showed low correlation with the phylogenetic distance of A. oligospora’s clades. These preliminary findings provided a strong support for the applicability of the vicariance hypothesis to this group of microorganisms. To complement this phylogeographic research effort, we used the same samples to evaluate the degree of genetic variability among strains and grouped them into haplotypes. We then mapped the results together with the current watershed map of Yunnan province as well as three putative historical configurations of the major rivers. We selected the paleo-Red river model proposed by Clark et al. 2004 [18], which is the most popular and widely supported one. We included here two additional models introduced previously: the one from Zhang and colleagues presented in two publications [10, 11] which supports and extends Clark et al.’s 2004 model, and the one from Wang et al. 2020 [27] which instead refutes the linkage between the Salween and Mekong with the Red river, but portrays a Yangtze river flowing into the paleo-Mekong (Fig. 1a). We aim at investigating the distinct evolutionary lineages (clades) found within each watershed and analyze those that are found in different ones in concurrence with the river change models. We hypothesize that the current distribution of clades reflects the fragmentation of their ancestral range caused by geological events modifying the watershed configuration of the region.

Materials and methods

Study region: current and putative past river network configuration

The study area covers Yunnan province of south-east China, a territory larger than 394,000 km2 which hosts the upper reaches of six major river watersheds (Fig. 1), listed here by their international name followed by their local name in parentheses: a tributary of the Irrawaddy (Dulongjiang), the upper Salween (Nujiang), the upper Mekong (Lancangjiang), the Red river (Honghe), the Yangtze river (Jinshajiang), and the upper Pearl river (Nanpanjiang). The province is characterized by a complex and varied landscape, with climate zones ranging from tropical in the south to temperate montane in the north-west, all influenced by the seasonal monsoon. Yunnan is located at the far-eastern edge of the Himalayan uplift which underwent a series of tectonic events beginning with the Indian subcontinent collision with the Eurasian plate. The surface uplift and the formation of mountains had an important impact on the configuration of water streams. Several researchers analysed the geological events in the region and proposed models of the historical drainage system.

Fig. 1

(a) Putative historical river configuration, with variations from three models [11, 18, 27]. The blue arrows indicate the flow direction of the streams. Note that this representation is just a sketch. The location of the connections as described by the respective authors are of limited precision and the geomorphology of the region (river lines, shape of the land) may have looked very different in past geological periods. (b) Yunnan province with its major river basins and the sampling sites where soil was collected, colour-coded by the evolutionary clades of A. oligospora [34]. The white points represent samples locations where this species was not detected or detected but failed DNA extraction or sequencing

Soil samples collection, treatment, NTF identification and sequencing

Because the soil sample treatment, culturing and molecular procedures, taxa identification and classification are all described in two previous works, we recommend the readers to consult the original articles [34, 35]. We include here information on molecular sequencing and focus on the additional steps done for the present analysis. After DNA extraction from the fungi’s mycelia and sequencing of amplified internal transcribed spacer (ITS), β-tubulin gene (TUB), translation elongation factor 1-alpha (TEF), and RNA polymerase II (RPB2) fragments, a total of 149 strains from 116 sites reached the standard for further analysis. These four genes are commonly used molecular markers in fungal research; they show some degree of variation between species, but are relatively conserved within species, and have a relatively constant rate of evolution. Despite that, ITS sequences are deemed suitable for single-species studies, while the other three genes have also exhibited significant discrimination capabilities among strains of our target nematophagous fungus [34, 36, 37]. The PCR products were sequenced by BioSune Biotech Company, (Limited, Shanghai, China) with the Sanger sequencing method on an ABI 3730xl sequencer. The quality of the sequences was assessed using the company’s protocol which sets a threshold of 20 for effectiveness. The obtained sequences were checked, edited and assembled using SeqMan v. 7.0 (DNASTAR, Madison, Wisconsin, USA) [38].

Phylogenetic and geographic analyses

The first step of our analyses consisted in reproducing the phylogenetic tree as in Deng et al. 2023 [34]. The sequence alignments of the four genes were generated using MAFFT v7 [39] and manually improved using BioEdit v7.7.1 (available at https://thalljiscience.github.io/). MEGA 6.0 [40] was used to link the four alignments into a multi-gene dataset. Arthrobotrys musiformis YL435 and Arthrobotrys vermicola MA95 were selected as outgroups. We used jModelTest v2.1.10 software [41, 42] to select the optimal calculating alternative model for the piecing sequence (ITS: GTR + I; TUB: SYM + I + G; TEF: GTR + I + G; RPB2: GTR + I + G). The phylogenetic tree was constructed following the Maximum Likelihood method (ML) using IQ-Tree version 1.6.5 [43]. The multi-gene dataset was partitioned, and each gene was analysed with the corresponding model. The statistical bootstrap support values (BS) were computed using rapid bootstrapping with 1000 replicates [44].

Afterwards, polymorphism of A. oligospora ITS, TUB, TEF, and RPB2 sequences were highlighted and classified into different haplotypes using DnaSP v6 [45] and then employed in downstream analyses. Insertion-deletion (InDel) mutations were coded as presence/absence characters except for gaps caused by mononucleotide repeats, which were excluded from analyses. We evaluated genetic relationships among haplotypes by the means of haplotype networks constructed for the four genes alone and combined using the median-joining network (MJN) inference method implemented in PopART [46]. MJN has the advantage to infer potential missing (unsampled) nodes which will shorten the total length of the network, decreasing complexity and yielding more readable structures [47]. The tolerance parameter epsilon (ε) regulating the generation of median vectors (additional nodes) was set to a conservative 0.

The haplotype network is a convenient method to represent the overall dataset but requires further analyses for better interpretation. To complement these results, we evaluated genetic diversity within and among watersheds using statistical metrics implemented in the DnaSP v6 software and Arlequin v3.5 [48], including number of individuals (n), number of haplotypes (Kh), haplotype diversity (Hd), nucleotide diversity (Pi), as well molecular variance (AMOVA) and the fixation index (Fst) for pair-wise comparisons, using 99,999 permutations. Neutrality tests to detect significant population expansions and bottlenecks in the evolution history of A. oligospora were performed using Tajima’s D [49] and Fu and Li’s F* statistics [50]. Moreover, population size changes were also evaluated by plotting mismatch distributions of the observed number of pair-wise differences compared to expected population growth-decline models. For this analysis, the initial effective population size (Theta initial, θ₀) and the number of generations considered (Tau, τ) were estimated assuming an infinite final effective population size (Theta final, θ₁).

To unveil more detail about haplotypes found in multiple watersheds and those found in watersheds not associated with their respective haplogroup, we mapped them over the current basin distribution of Yunnan and performed a visual assessment according to the three proposed hypotheses of historical river configuration.

Results

Our first result is the reproduction of the phylogenetic tree based on our 149 individuals of A. oligospora collected across Yunnan province. The subdivision in five main clades reported by Deng et al. was confirmed. The support values for the main branches ranged from 63 to 99%, indicating moderate to high confidence of the grouping outcome (Fig. S1).

For the following analyses, as the results for all four genes exhibited significant similarity, we focus on the combined data in this section. For reference, the Supplementary Information provides a breakdown of individual gene analyses. Among the 149 individuals, a total of 35, 31, 33, and 25 haplotypes were identified from the ITS, TUB, TEF, and RPB2 sequences, respectively; and when the genes were concatenated, the number of haplotypes reached 84. In all cases, five main clusters could be discerned. As showed in Fig. 2, these haplogroups match relatively well the six watersheds of Yunnan, with the exception of Irrawaddy and Salween forming a single group. Every haplogroup includes one or more haplotypes common to several samples but characteristic to their respective watershed. The most frequent haplotype (H3) was found in 26 individuals (17.4%) distributed in all watersheds except for the Irrawaddy, but mostly in the Mekong (14 individuals) and in the Yangtze (9 individuals). Four other haplotypes were also relatively common and widely distributed: H10 was found in 9 individuals and three watersheds, while H4, H20, and H26 were slightly less common but more widespread, found in seven to eight individuals and four watersheds, respectively. The remaining haplotypes were mostly found in single individuals (71 haplotypes/individuals), or in more than one individual but endemic to a single watershed (H13, H5, H17, H60, H79).

Fig. 2

Network of relationships among the 84 haplotypes found in A. oligospora. Missing haplotypes are represented by black and white circles. The sizes of circles are approximately proportional to sample size. Branch lengths are indicated by perpendicular bars along branches

Haplotype diversity was high in each watershed, ranging from 0.882 in the Mekong to 0.956 in the Pearl, but nucleotide diversity was relatively low (from 0.00803 in the Red to 0.01218 in the Pearl), indicating within-watershed closely related haplotypes (Table 1). However, despite a general pattern comparable to the concatenated sequences, DNA polymorphism metrics revealed slightly divergent patterns between genes (Table S1). AMOVA returned 34.52% of variation among watersheds (Salween and Irrawaddy merged) and 65.48% within watersheds (Table S2). Pair-wise analyses showed significant differences among populations (watersheds), with the most similar being Pearl and Mekong (Fst = 0.1367, p < 0.001) and Yangtze and Red (Fst = 0.13005, p < 0.001), while the greatest variation was found between Red and Salween-Irrawaddy (Fst = 0.62977, p < 0.001). This latter watershed showed greater genetic variations with all other watersheds. Figure S3 in the supplementary Information summarizes the pair-wise analyses for single-gene and the four genes combined.

Table 1 DNA polymorphism and neutrality tests of A. oligospora sequences (ITS, TUB, TEF, and RPB2 genes combined) in the six major watersheds of Yunnan province, China (Salween and Irrawaddy merged). Number of individuals (n), number of haplotypes (kh), haplotype diversity (hd), nucleotide diversity (pi), sampling variance of pi (s2), average number of nucleotide differences (k)

The neutrality tests resulted non-significant in each watershed and in the overall region (Table 1) and the observed mismatch distributions were all multimodal suggesting that A. oligospora didn’t experience population expansion or bottlenecks in its history but separated in several sub-populations (Fig. 3 and Fig. S2-S6 for single-watersheds).

Fig. 3

Population size changes analysis for 149 A. oligospora individuals (ITS, TUB, TEF, and RPB2 gene sequences) found in Yunnan province, China, based on pair-wise differences

When overlaid on the putative past river configurations, we can observe that most of the mismatching haplotypes are located in the areas where the drainage changes are thought to have occurred, especially in northwestern Yunnan (Fig. 4a and b), except for a cluster of haplotype H3 in the lower Yangtze river related to haplogroup 3 (Mekong). Those haplotypes closely related to the Yangtze and the Mekong (belonging to their haplogroups) but found in the other watershed, respectively, support both Clark et al. 2004 and Wang et al. 2020 theories. Moreover, the distribution of H3 and H10, clustered in northwestern Yunnan but also found in the current lower Jinsha (Yangtze) and all along the Lancang (Mekong), provide a significant testimonial of the relationship between these two watersheds. The presence within the Mekong watershed of haplotypes associated with the Red river (H2, H29, H40, H43) suggests a relationship between them as proposed by the paleo-drainage model of Clark et al. 2004 but rejected by Wang et al. 2020. Also, H4 highlights the communication between the upper Lancang (Mekong) and the Red watershed (Fig. 4a) while supporting the relationship Yangtze-Mekong (Fig. 4b). The hypothesis of Zhang et al. 2011 propounding an historical connection Yangtze-Pearl-Red in the eastern part of what is now Yunnan province is also geographically pertinent when observing their mismatching and shared haplotypes (Fig. 4c). In particular, H4 which is found in these three watersheds, and H1 and H66 which are shared by two watersheds are almost exclusively found along the putative past river lines.

Fig. 4

Hypotheses of historical river connections and geographic distribution of haplotypes not matching the current watershed configuration or found in multiple watersheds according to the inferred haplotype network (Fig. 2). Only haplotypes relevant to each hypothesis are represented in each map. (a) Dendritic pattern proposed by Clark et al. 2004 [18]; (b) Yangtze connected to Mekong hypothesis by Wang et al. 2020 [27]; (c)Yangtze-Pearl-Red rivers connection by Zhang et al. 2011 [11]. The cluster of H3 and H10 illustrated in (b) is also valid for (a) but not indicated for better visualization

Discussion

The phylogeographic analysis presented in this work revealed that the genetic origins of a widely distributed microorganism species align closely with the boundaries of major watersheds. The haplotype network, which depicts genetic relationships, could be divided into five haplogroups that correspond remarkably well to the watershed structure. The metrics describing genetic variations within and between watersheds are in line with the linkages shown in the haplotype network and confirm the strong watershed pattern in the genetic evolution of A. oligospora, previously put forward by Deng et al. 2023 [34]. However, within each haplogroup there are few haplotypes found in a different watershed or even in multiple watersheds, a “mismatch” that prompted us to test for potential drainage reorganisations. In particular, haplogroup 5 is geographically associated with the Salween-Irrawaddy watershed but it includes haplotypes from the Yangtze, Mekong and Pearl rivers. The former two suggest genetic relationships between the watersheds as proposed by Clark et al. 2004 [18]. This may be confirmed by looking at haplogroup 1 (Yangtze) and haplogroup 3 (Mekong), which both include branches connecting to haplotypes from the Salween (H3, H39). Moreover, they show similarities with the Red river (i.e., H3, H10, and H83) while the haplogroup associated with the Red river has connections with haplotypes found in the current Yangtze (H20, H61, H78), Mekong (H12, H20, H29, H40, H43, H52), and Salween (H20). Most interestingly, the shared haplotypes H3, H4, H5, H10, H20, and H26 highlight these putative past relationships between the four watersheds. The haplotype network also supports the model of Wang and colleagues [27], pointing to a connection between the upper Jinsha river (Yangtze) and the Mekong, as shown by the mixture of haplotypes found in their respective haplogroups. The third hypothesis suggesting a Yangtze-Pearl-Red watercourse [11] can also find support in the network. H1 and H33 in haplogroup 4 (Pearl) are located in the Red river watershed while haplogroup 1 (Yangtze) includes haplotypes from the Pearl watershed (H66 and two samples of H4).

These findings have important implications for our general understanding of microbial biogeography. The existence of biogeographic patterns for microorganisms, in the sense that their spatial distribution is determined by both contemporary environmental factors and major geological forces limiting dispersal is becoming more evident in recent years [51,52,53,54,55]. Authors at the front line of this research direction argue that, although the rules and patterns known for macroorganisms may be useful to study microorganisms, different approaches should be considered to face the challenges presented by this important component of the living world [53, 56, 57]. Clearly, variables related to geological history should be better included in the assessments. For example, in their study over an elevational gradient of 3000 m in Galongla mountain in China, Hu et al. 2020 [58] found that the inclusion of geological effects, such as parent rock and weathering, increased the explanation potential of soil bacterial communities distribution by 35.9%. However, most studies and reviews on nematophagous fungi focused mainly on the role of biotic and abiotic factors driving their spatial distribution with little consideration for geographical distance and isolation [59,60,61]. Previous analyses on the genetic distribution of A. oligospora at the national scale of China found historical population divergences between two clusters, a pattern that was mainly attributed to local environmental adaptations, while also reporting recent hybridization by long distance dispersal [62, 63]. Interestingly, the authors observed higher genetic diversity with several unique alleles and genotypes in samples located in southwestern China confirming our results of DNA polymorphism. The authors recognized that the harsh relief of this region may be the main factor leading to local adaptations, however, geographic distance was a poor explanatory factor in their study. We suspect that the sampling strategy employed in those studies was not adequate to test for the vicariance hypothesis. For this reason, the effect of topography should be accounted when calculating distances, for example by using geographic distance instead of Euclidean distance when compiling matrices, or using a stratified sampling approach that considers relief. In the present study and in Deng et al. 2023 [34], major river basins were used as spatial analytical units and the results revealed a well-marked distributional pattern, suggesting that vicariance was a leading factor promoting biological evolution and that high mountain ranges could effectively isolate populations, limit dispersal and lead to distinct genetic drift. Therefore, phylogeographic analyses designed with adequate sampling effort and stratification can reveal restricted geographic distributions of genetic traits of microorganisms otherwise believed to be widespread and with unlimited dispersal. These findings strongly support the inclusion of watersheds as biogeographical units to improve the delimitation of the existing and commonly used biological provinces, a solution that may not only bring microbes into global biogeography discussions but also suit the vast majority of living organisms [64]. Watersheds are relatively contained systems characterized by specific geomorphological and environmental features interacting with each other and with neighbouring units at different degrees of magnitude. Their natural organisation in subsystems (river branches) allows for multi-scale analyses that do not rely on arbitrary unit subdivisions such as regular grids. Hence, they represent a suitable and meaningful eco-geographical partitioning system with potential applications not only for ecological research but also for their integrated environmental and socio-economic management.

The haplotype network analysis and mapping exercises reported here corroborate previous efforts to reconstruct past geological events, in particular drainage reorganisation. Taken alone, the present study may be speculative; the mismatched distribution of A. oligospora’s haplotypes may be explained by other opportunities for gene flow such as effective dispersal by any kind of vector, decreasing the differences among watersheds. However, the fact that our results targeting a microscopic fungus showed similar outcomes to analyses performed on other organisms [13,14,15,16,17] while also being coherent with conclusions coming from other research fields (i.e., geology, sedimentology, and geomorphology), make up a significant volume of evidence supporting the three reconstructions of past drainage pattern. Because we cannot refute any of these three hypotheses and we are unable to provide specific time estimation for A. oligospora’s diversification (because of the absence of fossil records for this particular group of fungi and the choice of the gene sequences) we may posit that all the proposed watershed connections existed at different time periods. In fact, while Wang and colleagues state that their model supports a paleo-drainage configuration that existed from the late Early Cretaceous to the Late Cretaceous (120 − 66 Ma) and was disrupted by surface uplift that occurred at that period, the other two models point to the more recent geological events of the Neogene (23 − 2.4 Ma). Therefore, integrating solid molecular clock estimations in future studies targeting genetic signatures of soil microbes may enhance our knowledge of past drainage systems.

This article aimed at discussing new perspectives on microbial biogeography, in particular the suitability of using watersheds as spatial analytical units in ecology. Although the nematode-trapping fungus targeted in this study seems to be ubiquitous in temperate regions, it shows a strong gene isolation potentially explained by vicariance events catalysed by surface uplift. This finding suggests that watersheds should be considered for developing a global biogeography theory that could finally include microorganisms, as well as planning natural resources management and conservation. Moreover, by leveraging the unique features of soil microbial communities and adopting creative investigation approaches such as phylogeographic analyses, researchers can unravel the intricate connections between drainage configurations, geological events, and the evolution of microbial life, providing valuable insights into Earth’s past and informing sustainable management practices for the present and the future.