Niche Differentiation and Co-Occurrence Network of Fungal Communities Associated with Host A liations in an Extremely Arid Desert Ecosystem

Yiling Zuo Hebei University College of Life Sciences https://orcid.org/0000-0001-6919-2408 Xia Li Hebei University College of Life Sciences Jingya Yang Hebei University College of Life Sciences Jiaqiang Liu Hebei University College of Life Sciences Xueli He (  xlh3615@126.com ) Hebei University College of Life Sciences https://orcid.org/0000-0002-2783-3390 Lili Zhao Hebei University College of Life Sciences


Background
Interorganism associations between plants and microbial communities represent symbiotic relationships, with these taxa having coevolved for centuries [1]. The microbial communities can be referred to as the host's second or extended genome and are therefore critical for the health of their multicellular hosts [2,3]. Fungal endophytes, as an important component of microbial diversity, play essential roles in regulating multiple ecosystem functions and enhancing ecosystem stability [4][5][6]. These assemblages present extensive host and ecological ranges and are even distributed in habitats that are stressful to plants, such as desert regions [7][8][9]. The establishment of particular fungus-plant mutualistic relationships can confer thermotolerance, drought resistance, and a multitude of functional capabilities that enhance the survival, primary productivity and community structure of plants [10][11][12]. Arid desert is a unique habitat type with extreme conditions (e.g., low water availability, high salinity, high irradiance, and low nutrient availability) that pose challenges for the present biota, and the reservoir of fungal biodiversity adapted to extreme arid desert conditions is incompletely described. Thus, more insights into the endophytic fungal diversity and distribution patterns in desert environments are needed to fully elucidate plant-fungus interactions.
Fungal microorganisms tend to establish strong links with their hosts and show high host speci city; the interruption of these associations may result in the loss of host adaptability [13][14][15]. As a consequence of the reliance of endophytic fungi on the plant's inner environment for growth, plant species identity therefore presumably plays a particularly decisive role in fungal colonization of internal tissues [16,17].
Several studies have documented the signi cant effects and selective pressure of host plants on endophytic fungal species and community composition [18][19][20][21][22]. For example, Yao et al. [19] found that plant identity had a greater effect on endophytic fungal community composition than epiphytic fungi in the leaves of six mangrove species in South China, and Martínez-García et al. [20] showed differences in the arbuscular mycorrhizal fungal communities under different shrubs in a semiarid environment. Plantfungus interactions have developed as a response to drought in desert ecosystems [21,22]. In arid desert areas, different shrub species generate islands of fertility and grow on resource-rich patches, which leads to the heterogeneous distribution of substrate resources and reveals the potential of the host genotype to induce the establishment of speci c microbial communities [15,20,23]. Considering the spatial differences in fungal communities and varied speci city among host species, it is therefore essential to investigate fungal community composition to select representative plant species in desert ecosystems that are distributed in the same area and that are exposed to the same fungal pool [24][25][26].
The above-and belowground tissues of plants potentially represent distinct microbial habitats, such as the soil rhizosphere and the root, stem and leaf endospheres, which may host contrasting fungal communities by providing speci c biotic and abiotic conditions for the resident microorganisms [3,27,28]. Disentangling the niche differentiation of fungal communities in different plant tissues is very useful for describing patterns of plant-microbe interactions [29]. Recently, most attention has been dedicated to niche differentiation of microbes at the rhizospheric soil-root interface [20,30], in the leaf and root endospheres [31][32][33], and even in different organs of plants [15,27]. For example, Martínez-García et al. [20] demonstrated that rhizospheric soil arbuscular mycorrhizal fungal communities showed great differences from the fungal communities in root habitats in desert shrubs. In addition, Guevara-Araya et al. [33] observed signi cant differences in the community composition of endophytic fungi between leaf and root tissues of Aristolochia chilensis in an arid ecosystem, and Durand et al. [15] revealed contrasting aboveground (leaf and stem) and belowground (soil and root) fungal communities collected from poplar at a Hg phytomanagement site using ITS ribosomal RNA (rRNA) gene pyrosequencing. Despite the evidence of the specialization of microbiota to their respective niches, the taxonomic overlap of these microbes within different tissues has been less documented. Bai et al. [34] established Arabidopsis leafand root-derived microbiota culture collections and found large overlap in the genome-encoded functional capabilities of leaf-and root-derived bacteria, which suggests the possibility of reciprocal relocation between root and leaf microbiota members. Nevertheless, research involving the overlap of fungal communities in different plant tissue niches is quite limited [35], especially in extremely arid desert ecosystems. Understanding niche differentiation and fungal interactions above and below ground is essential to explain how these symbioses confer bene ts and ecological adaptations to plants in desert areas [36].
Multispecies assemblages usually do not live in isolation but instead form interactive relationships through mutualistic or competitive connections [37]. These interactive connections ultimately result in the construction of complex interspecies networks that maintain the structure of an ecological community and consequently the functions of the ecosystem [38,39]. Describing the architecture of these networks may help to decipher underlying species coexistence and ecosystem stability [40]. Co-occurrence network analysis has proven to be a useful tool for analysing microbial relationships and has recently been used to study complex ecological systems, such as bacterial or fungal communities in various environments [41][42][43]. The characteristics of keystone species, member interactions and community organization can be predicted by network analysis to determine the importance of operational taxonomic units (OTUs) and their associations in the community [44,45]. Recently, the taxon co-occurrence patterns of fungal networks have received increased levels of attention. Durand et al. [15] compared the fungal network structures in different tissue habitats of poplar and indicated that 80% of the soil fungal species were negatively correlated revealing the dominance of a pattern of mutual exclusion. However, there remains a lack of studies assessing the links across different niche habitats, which could provide new insights into the way fungal interactions are organized in above-and belowground environments.
The Anxi Extreme Arid Desert of the National Nature Reserve in Gansu Province, Northwest China, is located at the intersection of an Asian temperate desert, an extremely arid desert and a typical desert region, and has sparse native vegetation composed of shrubs [46]. Ephedra przewalskii, Reaumuria soongorica, Sympegma regelii, Nitraria sphaerocarpa and Salsola passerina are typical endemic species of extreme arid deserts in Northwest China. These xerophytic species, characterized by exceptionally high drought resistance, are challenged by several forms of abiotic stress and highly adapted to their environment. The long-term evolution of desert ecosystems has resulted in the development of special characteristics and unique associations of endophytes with these host plants. It is speculated that bene cial endophytes confer the ability of the host to tolerate drought and allow these plants to become pioneers in extremely arid environments [47,48]. To the best of our knowledge, no reports on the distribution patterns of the endophytic fungi associated with these prevalent plant species in extreme arid deserts are available.
Our previous investigations predominately focused on the spatial patterns of desert plant root-associated fungi, such as arbuscular mycorrhizal fungi and dark septate endophytes [49][50][51][52][53], and demonstrated the geographic and host a liations of fungi that are symbiotic with ecologically important plants in arid desert. However, there are still important gaps in our understanding of the occurrence and distribution pattern of fungal endophytes in different tissues of desert plants. In the present study, the genetic diversity of the endophytic fungal microbiota and its relationships with the plant communities in an extremely arid desert in Northwest China were determined using Illumina MiSeq sequencing techniques.
We hypothesize that plant species in uence the structure of the endophytic fungal communities in desert areas. We also characterized the endophytic fungal community composition in different tissue niches, such as in the leaf, stem and root tissues, and analysed the niche differentiation and taxonomic overlap in the above-and belowground environments. Co-occurrence network analysis was also used to investigate the composition and assembly of the fungal communities in these samples. We addressed the following questions: 1) How variable are the fungal communities associated with different shrubs within the same area in the extremely arid desert? 2) Do the endosphere fungal communities vary among the different plant tissue niches? 3) How do the aboveground fungal communities relate to the belowground fungal communities?

Results
Characterization of Illumina sequencing data After removing sequences of low quality, potential chimeras and singletons, the remaining non-chimeric fungal internal transcribed spacer 2 (ITS2) sequences (2,783,159 in total) were clustered into 504 operational taxonomic units (OTUs) at a 97% sequence similarity level. Of these 504 OTUs, 421 (2,692,818 reads) were identi ed as fungal OTUs. We then excluded the OTUs with < 10 reads from all the samples, and the number of sequences per sample was normalized to the smallest sample size (39,948 reads), ultimately resulting in a normalized dataset composed of 337 fungal OTUs (1,677,816 reads). The identi ed fungal OTUs included 280 Ascomycota, 54 Basidiomycota, and 3 Rozellomycota, representing 3 phyla, 20 classes, 49 orders, 102 families, 163 genera and 219 species. A total of 124 relatively abundant OTUs (> 1000 reads) accounted for 99.9% of the reads for endophytic fungi (Additional le 1: Table S1). The DNA sequences of Dothideomycetes (29.84-81.92%) were the most abundant at the class level, and those of Pleosporales (29.83-74.76%) were the most abundant at the order level, with varying relative abundances among all samples (Fig. 1a, b). The members of Sordariomycetes (18.27-46.34%) were mainly distributed in root tissue niches, while members of Tremellomycetes (2.21-30.59%) were mainly distributed in leaf and stem niches (Fig. 1a). Rarefaction curves for the Sobs index at the OTU level across all samples approached an asymptote and showed that the overall fungal diversity was well represented (Additional le 2: Figure S1).

OTU richness and alpha diversity
Sobs index values were calculated to describe the observed OTU richness. Alpha diversity, including the Shannon diversity index and evenness, was analysed based on OTU richness (Fig. 2). The OTU richness of endophytic fungi in the stems, leaves and roots ranged between 19.0-52.3, 20.3-30.7 and 24.7-104.7 (mean values), respectively, and the root endophytic fungal communities were much more diverse than the aboveground communities (Fig. 2a). The Kruskal-Wallis test revealed that plant species had a signi cant effect on the OTU richness of the stem endophytic fungi (χ 2 = 10.939, P = 0.027) but not on that of the leaves (χ 2 = 3.308, P = 0.347) or roots (χ 2 = 5.957, P = 0.114). For example, the OTU richness of the stem endophytes in N. sphaerocarpa was signi cantly higher than that in E. przewalskii, R. soongorica and S. regelii. In addition, the integrated OTU richness was highly dependent on plant compartment effects (χ 2 = 6.612, P = 0.037). The OTU richness of the stem niche in N. sphaerocarpa was comparable with that of the leaf samples from S. regelii, N. sphaerocarpa and S. passerina and root samples from S. regelii and N. sphaerocarpa (Fig. 2a). In terms of the Shannon diversity and evenness estimates, we also found signi cant variation in the stem tissues among the 5 plant species (P < 0.05) (Fig. 2b, c).

Community composition of endophytic fungi
The composition and overlap of endophytic fungi among the plant species and tissue niches were documented (Additional le 2: Figure S2). Among the identi ed fungi, R. soongorica represented the highest percentages of unique OTUs in the leaf (29.08%) and root (36.30%) tissues, and only 7 OTUs (1.16 and 3.57%) were shared by the 5 plant species (Additional le 2: Figure S2b, c). In addition, N. sphaerocarpa and S. passerina exhibited the highest proportion of unique OTUs (26.97 and 28.29%) in the stem samples, with 16 shared OTUs (5.26%) occurring among plant species (Additional le 2: Figure  S2a). Furthermore, 81 OTUs (12.05%) were shared by the leaf, stem and root tissue niches, with the root niche exhibiting the highest percentage of unique OTUs (41.54%) of endophytic fungi (Additional le 2: Figure S2d). Notably, the OTUs of leaves (10.36%) and stems (10.48%) shared greater proportion of taxa with roots than with each other (2.46%).
A heatmap was constructed based on the top 50 most abundant OTUs to reveal the varied occurrence of relatively abundant fungal OTUs among the different plant species and niche habitats (Fig. 3). Among the endophytic fungi distributed in stem tissues (Fig. 3a), four OTUs (Ascomycota OTU480 and OTU212, Sirobasidiaceae OTU173 and Hazslinszkyomyces lycii OTU398) were distributed dominantly in N. sphaerocarpa, and four OTUs (Fusarium proliferatum OTU400, Amphobotrys ricini OTU123, Ascomycota OTU144 and Pleosporales OTU120) were present mainly in S. regelii, while Sporidiobolus pararoseus OTU2 was exclusively found in E. przewalskii. The distribution of leaf fungal OTUs was signi cantly biased and showed preferences for particular plant species (Fig. 3b). For example, the fungal OTUs (Thyrostroma OTU11, Acremonium chrysogenum OTU167, Sporormia grandispora OTU166 and Dothiora OTU9) occurred in only R. soongorica; six OTUs (Fusarium OTU132, Periconia byssoides OTU133, Sporormiella septenaria OTU142, Cystobasidiomycetes OTU104, Ascomycota OTU144 and Dematiopleospora OTU476) were mainly found in S. regelii; Dothideales OTU21 and Ascomycota OTU473 occurred in N. sphaerocarpa; and three OTUs (Apiotrichum montevideense OTU128, Udeniomyces kanasensis OTU390, Neocamarosporium OTU393) were present only in S. passerina. The detectable OTUs in roots with high relative abundance were generally shared by multiple hosts, except Sordariales OTU436 and Preussia polymorpha OTU442, which were restricted to the root samples of R. soongorica, and Pleosporales OTU 273, which was exclusively observed in the root of N. sphaerocarpa (Fig. 3c). More fungal OTUs with lower abundances were strongly associated with roots of different plant species (Fig. 3c). Furthermore, the relative abundances of some abundant OTUs in different tissue niches signi cantly differed among plant species (P < 0.05), with the most fungal richness differences occurring in the stem niche (Additional le 2: Figure S3a, b, c). The heatmap also displayed the biased occurrence of relatively abundant endophytic OTUs in above-(stem and leaf) and belowground (root) tissues ( Fig. 3d), which corresponds with the observed signi cant difference between the clustered above tissue OTUs (stem and leaf) and the root OTUs. Plant tissue niche also exhibited signi cant effects on relatively abundant OTUs; for example, Thyrostroma OTU10, Aporospora OTU497 and OTU445, and Sordariales OTU278 signi cantly differed among plant tissues (P < 0.001) (Additional le 2: Figure S3d).
Nonmetric multidimensional scaling (NMDS) ordination was used to visualize the differences in endophytic fungal communities among plant species and tissue niches (Fig. 4), and signi cant dissimilarity among all factors was con rmed with the ANOSIM test. The endophytic fungal community composition of the different niches signi cantly differed among the ve plant species (stem: F = 0.7911, P = 0.001; leaf: F = 0.4090, P = 0.005; root: F = 0.2756, P = 0.024) and presented characteristics of permutation dispersion with higher heterogeneity of the OTU distribution in different hosts ( Fig. 4a, b, c). Moreover, the Bray-Curtis method revealed a signi cant effect of the plant compartment (F = 0.2556; P = 0.001), and the stem and leaf niches were clearly distinguished from the root tissues but did not cluster completely according to their respective niches (Fig. 4d).

Linear discriminant analysis of effect size
We used linear discriminant analysis (LDA) of effect size (LEfSe) to identify the taxa that most likely explain the differences among plant species and tissue niches. Indicator fungal species in the communities with LDA scores of 4 or greater were identi ed as the best biomarker for each category (Fig. 5). The signi cant host effect observed for the stem fungal communities was explained by 27 indicator fungal species: 5 for E. przewalskii, 12 for N. sphaerocarpa, 6 for R. soongorica, 1 for S. passerina, and 3 for S. regelii (Fig. 5a). In contrast, 8 and 10 indicator taxa were signi cantly enriched in the leaf and root fungal communities associated with the different plant species, namely, unclassi ed Ascomycota and Pleosporaceae in the leaves of N. sphaerocarpa; Filobasidiales in the leaves of S. passerina; Cystobasidiomycetes, Bulleribasidiaceae and Tremellales in the roots of E. przewalskii; Agaricomycetes in the roots of N. sphaerocarpa; Trichomeriaceae and Comoclathris in the roots of R. soongorica; and Pyronemataceae, Lasiobolidium, unclassi ed Ascobolaceae, and Trichophaeopsis in the roots of S. regelii (Fig. 5b, c). Among the plant tissue niches, root samples presented the most abundant species with signi cant differences (Fig. 5d), evidencing that the strong dominance of root fungal microbiota (38 indicator taxa) explains the niche differentiation of fungal communities in desert regions.

Co-occurrence networks of endophytic fungi
Co-occurrence analysis resulted in biased patterns of fungal clustering for the different niche habitats and provided information on their corresponding network properties in terms of network connectivity, the clustering coe cient and network density (Fig. 7a, b, c; Additional le 3: Table S2). The root habitat showed the highest network connectivity, exempli ed by the highest numbers of nodes (120) and edges (1385) (Fig. 7a, b, c; Additional le 3: Table S2). The root network also had signi cantly higher mean values of network density (0.194 versus 0.098 and 0.100) and centralization (0.247 versus 0.181 and 0.134) than the stem and leaf networks, while the clustering coe cient for the leaf habitat was higher than those for the stem and root networks (Additional le 3: Table S2). Notably, the analysis of the relationships between fungal OTUs in the whole dataset revealed a high percentage of co-presence rather than mutual exclusion (Additional le 3: Table S2). OTUs from Pleosporales were the most abundant in all three networks, accounting for 46.91% of the total OTUs in stems, 51.90% of the total OTUs in leaves and 39.17% of the total OTUs in roots (data not shown). According to the number of connections established with the rest of the network, 10 dominant keystone OTUs were arbitrarily identi ed here for each habitat based on the number of degrees for each node (Additional le 3: Table S3). We found that the taxa with a high number of degrees but relatively low abundance served as important nodes in terms of maintaining the function and structure of the microbial community (Additional le 3: Table S3). MCODE analysis returned several signi cant clusters of fungal groups (6 in stems, 6 in leaves and 7 in roots) that were particularly closely connected in the network (Additional le 2: Figure S4), indicating the typical small-world characteristic of the three networks in this arid desert.
To take into account the interactions between the above-and belowground fungal communities, an integrated co-occurrence network was constructed to interpret the relationships among the fungal communities across the different niches (Fig. 7d). The construction of this network led to the identi cation of 112 nodes and 441 edges, with the highest values of characteristic path length and network heterogeneity (Additional le 3: Table S2). In particular, the full network contained all OTUs exhibiting positive connections and showing a 100% co-presence percentage (Additional le 3: Table S2). The fungal nodes in the network tended to be shared by all three niche habitats (51.79%), which con rmed the results in terms of the relationships and overlap among the different plant tissues. The nodes shared by leaves-roots and stems-roots accounted for 14.29 and 11.61% of the total, respectively, while those shared by the leaves and stems accounted for 0.89% of the total. Additionally, 21.43% of the nodes appeared to be speci c to the root niche habitat. Although the fungal network tended to contain the most OTUs belonging to Pleosporales (47 nodes), eight Hypocreales members were identi ed as dominant keystone OTUs (Additional le 3: Table S3). MCODE analysis returned eleven signi cant clusters with network scores ranging between 3 and 18, and the Hypocreales groups produced the rank 1 cluster, which had 20 nodes and 171 edges (Additional le 2: Figure S4).

Discussion
In the present study, a total of 337 OTUs of endophytic fungi were identi ed at a 97% sequence similarity level. This level of endophytic fungal diversity associated with shrubs growing under extreme desert conditions in Northwest China was signi cantly lower than that associated with plants in subtropical forests, farmland plantations and even wetland ecosystems [19,37,41]. Considering the unique habitat provided by arid areas, it seems that the typical harsh environmental conditions of extreme desert lands truly reduce the number of plant endophytes. Rarefaction analyses and richness estimators indicated that much of the total diversity detectable with Illumina-based sequencing was observed. Although a great number of endophytic fungi have been isolated from desert plants worldwide [6,12,[54][55][56][57], this study is the rst to describe plant-fungal endophyte associations in an extremely arid desert ecosystem with the simultaneous consideration and comparison of above-and belowground environments.
Thanks to their various life forms and ecological plasticity, fungal assemblages display strong adaptation to extreme environments, and their role is particularly crucial for the recycling of organic matter and favouring nutrient uptake under stressful conditions [58]. Our study showed that Ascomycota were dominant among the desert plants, which is in accordance with the ndings of some previous studies [57,59,60]. Most Ascomycota members are saprophytic and therefore the main decomposers in habitats, as they can decompose large amounts of refractory organic matter, thus playing an important role in desert nutrient cycling [61]. In this study, the DNA sequences of Dothideomycetes (29.84-81.92%) were the most abundant at the class level, and those of Pleosporales (29.83-74.76%)were the most abundant at the order level. Members of Dothideomycetes have been reported to be abundant in endophytic fungal communities in other plant taxa, such as poplars and grapevines (Vitis vinifera) [62,63]. Pleosporales have been ubiquitously detected in arid and semiarid environments [64], which is consistent with our results. Moreover, we also found that members of Sordariomycetes (18.27-46.34%) were distributed mainly in root tissues, while Tremellomycetes (1.21-30.59%) were distributed mainly in leaves and stems. This differentiation of the fungi in above-and belowground tissues could possibly be the result of resource utilization, which drives the rapid evolutionary radiation of microbes competing for resources, causing them to diverge and adapt to different niches to reduce competition [18,65].
We found that fungal OTU richness was highly dependent on the plant compartment (χ 2 = 6.612, P = 0.037) and that the root endophytic fungal communities showed relatively higher OTU richness than those in the aboveground tissues. Studies have demonstrated that plant roots can be highly colonized by abundant fungal endophytes and that the deeply and extensively distributed underground root systems of desert plants provide more hosts and substrate for infection by endophytic fungi than aboveground tissues [66]. Although 16 OTUs in the stem tissues were shared by 5 hosts, plant species signi cantly affected the OTU richness of the stem endophytic fungi (χ 2 = 10.939, P = 0.027), which indicates the high variability in endophytic fungi in the stem niche among the 5 plant species. This high variability may be related to the ability of endophytic fungi to penetrate plant tissues and ultimately successfully colonize the endodermis, the pericycle and even the xylem vessels [3,18,67]. In particular, it is worth noting that the leaves (10.36%) and stems (10.48%) showed a greater percentage of shared taxa with those in the root habitat than with each other (2.46%) in our study, which indicates the existence of taxonomic overlap in the endophytic fungal assemblages of the root and aboveground compartments. Zarraonaindia et al. [68] investigated the grapevine (Vitis vinifera)-associated microbiota and found that the microbiotas of the leaves, owers, and grapes shared greater proportions of taxa with the belowground communities (albeit with the soil microbiome) than with each other, suggesting that soil served as a common reservoir for plant microbiotas. Soil harbours an extraordinarily rich diversity of microorganisms and largely determines the initial inoculum of endophytic root-associated microbiota [34,69,70]. In arid deserts, the abundant root endophytes may also provide an inoculum source for the colonization of aboveground tissues, thereby leading to taxonomic overlap between aboveground and root tissues.
As expected, plant identity signi cantly affected the community composition of endophytic fungi, as suggested by our results from the heatmap of fungal relative abundance, the ANOSIM of fungal community structures and the LEfSe analysis. Although host plants in the same arid habitat seem to receive similar fungal propagules [20,34], differences in host chemistry and nutrient contents can affect endophytic fungal colonization and promote differences in the dominant members of the fungal community [19,71,72]. Another explanation for the differentiation of the fungal community among plants can be attributed to host/fungus preferences. Our host/fungus preference analysis demonstrated that four plant species (except R. soongorica) showed signi cant preferences for endophytic fungi and that 70 of 124 endophytic fungal OTUs displayed signi cant host preferences. The ve extremely endemic plants in this study belong to different families, and previous studies have demonstrated host preferences of endophytic fungi at the plant family level [73,74]. To a certain extent, the enrichment and depletion of speci c endophytic microbiomes are not passive processes but rather depend on the active selection of microbial consortia by the plant host [3,69,75]. Furthermore, the ability of endophytic fungi to successfully colonize a host plant requires the presence of speci c traits, such as the production of cell wall-degrading enzymes, and intricate interplay between the fungi and the host plant's innate immune system [3,66,67]. Some endophytic fungi can produce enzymes that dissolve a variety of plant cell walls, allowing them to enter and colonize a wide range of hosts, but a fungal species presenting host speci city may only produce a speci c enzyme that permits the endophyte to infect certain host plant cells [76,77].
Compared with plant species, the structure of the endophytic fungal communities was more obviously differentiated in association with the different tissue niches in our study. The clustering heatmap and the Bray-Curtis method (F = 0.2556, P = 0.001) revealed the signi cant occurrence of endophytic fungal communities in the above-(stem and leaf) and belowground (root) compartments, which is in agreement with previous studies by Tardif et al. [78] and Xiong et al. [18], who demonstrated signi cant plant compartment effects among plant endophytic microbes. The signi cant difference between the shoot and root niches observed here may be due to the in uence of fundamentally different conditions between the above-and belowground environments in the desert ecosystem [61,79]. The light intensity, humidity and temperature of the habitats in which the different microbes survive and the chemical activity of the different tissues against fungal infections can absolutely lead to differentiation in the endophytic fungal communities and niche speci city [80,81]. Although a proportion of fungal OTUs were shared by all the plant tissues, the niche speci city of fungi is highlighted by the existence of speci c proportions of OTUs that were exclusively found in the different plant niches. Each of the plant ecological niches provides relevant biotic and abiotic gradients, such as in the availability of soluble organic compounds [3,19,69], and the different ecological processes can also contribute to the high degree of compartment speci city among the fungal endophyte communities [33,82]. Here, the tissue/fungus preference analysis showed that all tissue niches presented signi cant preferences for endophytic fungi and that 88 of 124 fungal OTUs revealed signi cant tissue preferences. In addition, 44 of 372 tissue niche and endophytic fungus pairs were observed to show signi cantly strong preferences. These results suggest that individual tissue niches as distinct substrates could re ect preferences for dominant taxa [83]; simultaneously, different endophytic fungi also vary in their capacity for utilizing or surviving on a speci c substrate [84,85]. For instance, the deposition of large amounts of carbon in root tissue, such as root exudates, could directly drive the colonization of root-associated fungi [3,86]. Furthermore, some leaf endophytic fungi have been shown to harbour corresponding receptors to ensure survival in the phyllosphere environment [3,87].
Co-occurrence analysis revealed that the root niche showed the highest network connectivity and the highest number of interactions (120 nodes and 1385 edges), suggesting that the root endophytic fungal network is more stable and functional than the stem and leaf networks in this arid desert. Studies have indicated that the interspecies relationships in fungal communities could be strengthened by upregulated drought or wetting disturbance [37]. Thus, the strong fungal interactions in the root habitat are likely one of the factors explaining the reasonable adaptation of desert plants to this particular soil type. In our study, MCODE analysis returned several signi cant clusters of fungal groups (6 in stems, 6 in leaves and 7 in roots) in the fungal networks, indicating the typical small-world characteristic of the networks in this arid desert. The links between OTUs in the different modules in the network can also contribute to more intensi ed interactions among fungal OTUs at a high organizational level of the network [37]. The edges in co-occurrence networks could represent positive or negative correlations and imply a biologically or biochemically meaningful relationship between fungal species [15,88]. In our study, the relationships between fungal OTUs in the whole dataset revealed a high percentage of co-presence. Studies have indicated that microbes can establish positive or negative correlations based on their environmental preferences [15]. In arid deserts, fungal assemblages may persist in the community by establishing reciprocal symbioses with each other to resist stressful conditions [15,19,89]. Additionally, the disconnected microhabitats and the presence of a high number of dormant cells in the bulk sand can also explain the identi cation of a dominance of co-presence interactions in the arid desert [11,90].
Keystone species are represented by the centralized nodes in the network, have a large in uence on 'information' transfer throughout the community, and act to stabilize the microbial community [91,92]. Interestingly, we found that taxa represented by nodes with a high number of degrees but with relatively low abundance were important in maintaining the function and structure of the microbial community (Table S3). This nding is in agreement with those from a study by Mandakovic et al. [89], who compared microbial co-occurrence networks representing soil bacterial communities along a western slope of the Andes in the Atacama Desert and revealed that highly abundant OTUs were not as important as some of the relatively rare but highly connected species in terms of network dynamics and stability. Xue et al. [93] demonstrated that the removal of keystone taxa from microbial co-occurrence networks signi cantly decreased the number of edges and average connectivity of the network. Although all of the networks shared similar node OTU names, no overlapping keystone OTUs were found in the networks among the different plant compartments (Table S3). This indicates that role shifts of fungal OTUs occurred in the different plant tissue niches [37,94]. However, the keystone species originated from nearly the same taxa at a higher phylogenetic level, for example, members of the order Pleosporales were keystone species in all networks, suggesting that the functional roles of nodes in the networks could be conserved at a higher phylogenetic level [37]. In addition to those of Pleosporales, members of Hypocreales, Chaetothyriales (Cyphellophora) and Eurotiales (Penicillium) played important central roles in connecting the above-and belowground fungal networks, indicating the potential functional features of these species in the studied desert ecosystem. Abundant studies have focused on the entomopathogenic effect of Hypocreales strains and their potential use as biocontrol agents of insects [95][96][97]; thus, the important connections of this fungus in the co-occurrence network may be related to the biological defence mechanisms of desert plants. Penicillium, a well-known Ascomycota taxon that secretes antibiotics, is also known for its tolerance to low water potential and could develop strong relationships with other members in the module through antibiotics or secondary metabolites [37,98]. Overall, the occurrence and assembly of the fungal communities in the studied desert habitat were found to be highly correlated with host species and tissue niches, and the identi cation of the dominant taxa and keystone species could provide essential information for developing strategies to manipulate desert plant-associated microbes for desert vegetation restoration [18,99].

Conclusions
A total of 337 operational taxonomic units (OTUs) of endophytic fungi were identi ed in our study. The taxa of Pleosporale were mainly dominated in desert plants and played irreplaceable role as keystone species in maintaining the stability of fungal networks. The occurrence and assembly of fungal communities in desert habitats were highly correlated with host species and tissue niches. Compared with aboveground tissues, root-associated fungi represented the most abundant reservoir of biodiversity in desert habitat, and the aboveground stems and leaves showed higher taxonomic overlap with underground root tissues than with each other. There existed strong and positive inter-species connections in desert fungal community, and the root fungal co-occurrence network harbored the highest connectivity and established higher number of interspecies interactions. Our study rstly reported the endophytic fungi communities of desert plants in an extremely arid desert with the simultaneous consideration and comparison of above-and belowground niches, and suggested the central importance of root-associated fungi to desert plant communities and ecological functions. Understanding the complex host-microbe interactions of desert plants could provide the basis for the exploitation of the plant-fungi associations in manipulation of shrubs' microbiome for desert ecological restoration.

Study site and sampling
The study was conducted in the Anxi Extreme Arid Desert of the National Nature Reserve in Gansu Province (40°06′-40°08′N, 96°17′-96°36′E), Northwest China. This study site is located at the intersection of an Asian temperate desert, an extremely arid desert and a typical desert region companying with a typical continental climate. The average annual temperature is 7.8-10 °C, and the average annual precipitation is no more than 52.0 mm, which is far less than the average annual evaporation of 2754.9 mm [48]. The soils are dark chestnut and sandy chestnut [100]. The vegetation is dominated by typical super-xerophytes, including Ephedra przewalskii, Reaumuria soongorica, Sympegma regelii, Nitraria sphaerocarpa and Salsola passerina, which is abundant and exhibits a heterogeneous distribution.
Sampling was carried out in July 2019, consisting of collecting leaf, stem, and root samples from three replicated plots. Within each sampling plot, ve healthy plants were randomly selected, with each plant individual at least 100 m away from another individual. Stem samples were collected from the branches of each plant (0.1 to 0.6 cm diameter) at the half-crown level [3]. In addition, all leaves from the sampled offshoot were collected to represent the leaf compartment (except in the case of Ephedra przewalskii). Before collecting the root samples, the upper layer of soil (approximately 1-5 mm) was removed to clear away litter, and the roots were subsequently obtained from a depth of 30 cm under the canopy of each plant. All samples were collected over a 1-day period to reduce any heterogeneity imparted by the climatic conditions and immediately placed in sterile plastic bags, labelled, and transported to the laboratory in an ice box. The ve subsamples collected from each plot were equally merged into one, and a total of 42 plant tissue samples were analysed in this study. All the samples were immediately surface sterilized and stored at -80 °C until used for DNA extraction.

Molecular analysis
The plant materials (leaves, stems and roots) were washed with deionized water and surface sterilized by consecutive immersion for 5 min in 70% ethanol, 2 min in 5% sodium hypochlorite, and 30 s in 70% ethanol, followed by three rinses in sterile distilled water. Subsequently, the treated samples were freeze dried using liquid nitrogen and prepared for DNA extraction. The genomic DNA of the endophytic fungi from the plant samples was extracted using the Invisorb Spin Plant Mini Kit according to the manufacturer's protocol (Stratec Biomedical AG, Birkenfeld, Germany). The nal DNA concentration and purity were determined by a NanoDrop 2000 UV-vis spectrophotometer (Thermo Scienti c, Wilmington, USA), and DNA quality was checked with 1% agarose gel electrophoresis.

Bioinformatics processing
Raw FASTQ les were quality ltered with Trimmomatic and merged with FLASH to obtain valid and highquality sequences on the basis of the following criteria: (i) The reads were truncated at any site receiving an average quality score < 20 over a 50 bp sliding window. (ii) Sequences whose overlap was longer than 10 bp were merged according to their overlap with a mismatch of no more than 2 bp. The sequences from each sample were separated according to barcodes (exactly matching) and primers (allowing 2 nucleotide mismatches), and reads containing ambiguous bases were removed. Operational taxonomic units (OTUs) were clustered according to a 97% similarity cut-off using UPARSE (version 7.0.1090, http://www.drive5.com/uparse/) with a novel 'greedy' algorithm that performs chimera ltering and OTU clustering simultaneously after dereplication and discarding all singletons [104]. The taxonomy of each representative sequence was analysed with the RDP Classi er algorithm (http://rdp.cme.msu.edu/) against the UNITE (version 8.0, https://unite.ut.ee/) database using a con dence threshold of 70% [105,106]. To eliminate the effects of the different numbers of sequences among the samples on the identi ed fungal community, the number of sequences per sample was normalized to the smallest sample size using the subsample command in MOTHUR. Subsequently, rarefaction curves were assembled, and the alpha diversity indices of OTU richness, Shannon diversity and evenness were calculated. The relative abundance of speci c fungal taxa on the basis of OTU, order and class was de ned as the number reads of a particular taxon as a percentage of the number of all reads in a sample [19].

Statistical analysis
All statistical analyses were implemented in R version 3.5.1. The Shapiro-Wilk test and Bartlett's test were employed to check the normality and homoscedasticity of the data, respectively [15]. As the endophytic fungal OTU richness data did not satisfy the normality of distribution and homogeneity of variance, the effects of plant species and tissue niche on the alpha diversity estimates for the fungal assemblage were examined using the Kruskal-Wallis (KW) test [19], followed by Welch's tests for paired comparisons between samples [107]. The numbers of OTUs that were shared between plants and tissues were visualized using Venn diagrams. The relative abundances of abundant endophytic fungal OTUs in the different plant species and tissues were depicted using the pheatmap function in the pheatmap package in R; then, differential abundance analysis was performed with the KW H test (P < 0.05), and STAMP software (v.2.1.3) was used [108]. Nonmetric multidimensional scaling (NMDS) was used to visualize the community composition dissimilarities of the endophytic fungi based on the Bray-Curtis method using the metaMDS command in the vegan package in R. To statistically support the clustering of fungal communities shown in the NMDS plot, the nonparametric ANOSIM test (an analogue of univariate ANOVA) was used to examine signi cant differences based on 999 permutations [15,109]. Linear discriminant analysis (LDA) effect size (LEfSe), which can take into account statistical signi cance and biological relevance, was conducted to identify the OTUs differentially represented between the different species and tissues by using the nonparametric factorial KW sum-rank test and LDA [110,111]. High LDA scores re ect a signi cantly higher abundance of certain taxa than that of others.
According to Toju et al. [112], we evaluated the host/tissue-fungus preference based on the interaction specialization index (d′) using the 'dfun' function in the R bipartite package [113]. The d' index measures how strongly a fungus deviates from a random choice among plant partners that are available at a study site [114]. Considering the di culty of estimating the host preferences of rare fungi, the d′ estimates of the abundant endophytic fungal OTUs (> 1000 reads) were used [19,112,114]. Firstly, the sample data matrix was binarized to a sample-level matrix (present-absent data) and then converted into a specieslevel matrix, in which rows denote plant species, columns represent fungal OTUs, and cell entries indicate the number of samples from which respective combinations of plants and fungi were observed.
Subsequently, the randomized species-level matrices were obtained based on 1000 permutations by shu ing the plant species labels in the sample fungal OTU matrix. The d′ value of each fungal OTU [115] was standardized as follows: standardized d′ = [d′observed − Mean (d′randomized)]/SD(d′randomized), where d′observed is the d′ estimate of the original data, and Mean (d′randomized) and SD (d′randomized) are the mean and standard deviation of the d′ values of the randomized data matrices. The standardized d' value was also calculated for each plant species based on the original and randomized data matrices mentioned above. Additionally, the two-dimensional preferences (2DP) seen in a pair of a plant species (i) and a fungal OTU (j) were quanti ed based on the species-level original and randomized matrices to evaluate to what extent each pair of each plant-fungus association was observed (counts) more or less frequently than would be expected by chance. 2DP (i, j) = [N observed (i, j) -Mean (N randomized (i, j))]/SD (N randomized (i, j)), where N observed (i, j) is the number of the samples from which a pair of a plant-fungal OTU association was observed in the original data, and Mean (N randomized (i, j)) and SD (N randomized (i, j)) are the mean and the standard deviation of the number of samples for the focal plant-fungal OTU pair across randomized matrices. The P value obtained from the preference analysis was adjusted based on the false discovery rate (FDR) [116]. Similarly, tissue × fungal OTU matrix were carried out, and the standardized d' value of both each niche tissues and each fungal species was also calculated, as well as the2DP estimators.
Network analysis was performed to identify the microbial co-occurrence patterns across the three plant tissues associated with the ve plants, respectively. Alternatively, we explored the interspecies relationship between above-and belowground niches by constructing a total co-occurrence network pattern using the overall database. Before network construction, we removed very rare OTUs and the abundant endophytic fungal OTUs (> 1000 reads) were used [15,88]. Spearman's rank correlations were calculated with all possible OTU pairs and robust correlations with correlation coe cient above 0.5 and statistical signi cance (P < 0.01) were noted [117]. The nodes of the network represent individual OTUs that have at least one signi cant co-occurrence relationship with other OTUs. The networks were explored and visualized with Cytoscape v3.5.1, and NetworkAnalyzer tool was used to calculate network topology parameters [118]. Structural attributes of networks such as the number of nodes and edges, graph density, network connectivity, the clustering coe cient, network centralization, characteristic path length and network heterogeneity were determined [119]. Modular structure and groups of highly interconnected nodes were analyzed using the MCODE application with standard parameters [118,120]. A module is a group of nodes connected more densely to each other than to other nodes outside the group. Possible keystone species were those that demonstrated the higher degree values [121].

Additional les
Additional le 1: Table S1. Molecular identi cation of endophytic fungi (> 1000 reads) at a 97% sequence  Table S2. Structural attributes of networks obtained through network analysis for different plant niches (stem, leaf and root) and the total networks. Table S3

Availability of data and materials
The raw data sequences wered eposited in the NCBI Sequence ReadArchive (SRA) under the Bioproject number PRJNA681429 and Bio Sample accession numbers SAMN16953820 to SAMN16953861.
Authors' contributions XLH and YLZ designed the experiment; YLZ and JQL collected the samples; XL, JYY and LLZ performed the laboratory work; XLH, YLZ and XL analyzed the data and wrote the manuscript. All authors read and approved the nal manuscript.
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Figure 2
Operational taxonomic unit (OTU) richness and alpha diversity estimates of the fungal communities. The black line inside each box represents the median value. Kruskal-Wallis test revealed that plant species identity had a signi cant effect on the OTU richness of stem endophytic fungi (χ2 = 10.939, P = 0.027), but not on that of leaf (χ2 = 3.308, P = 0.347) and root compartment (χ2 = 5.957, P = 0.114). Alternatively, OUT richness was highly dependent on plant compartment (χ2 = 6.612, P = 0.037). Paired comparisons were performed by Welch's t tests for between samples. RS, Reaumuria soongorica; SR, Sympegma regelii; NS, Nitraria sphaerocarpa; SP, Salsola passerina; EP, Ephedra przewalskii. S, stem; L, leaf; R, root. *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001.       Analysis of host-tissue/fungus preferences. The standardized d' estimate of preferences for fungal operational taxonomic units (OTUs) are shown for each hosts and tissue niches (column). Likewise, the standardized d' estimate of preferences for plant species and tissue niches are indicated for each of the observed fungal OTUs (row). Each cell in the matrix indicates a two-dimensional preference (2DP) estimate, which evaluates to what extent each pair of each plant/tissue-fungus association was observed (counts) more or less frequently than would be expected by chance. The P values were adjusted based on the false discovery rate (FDR). RS, Reaumuria soongorica; SR, Sympegma regelii; NS, Nitraria sphaerocarpa; SP, Salsola passerina; EP, Ephedra przewalskii.

Figure 6
Analysis of host-tissue/fungus preferences. The standardized d' estimate of preferences for fungal operational taxonomic units (OTUs) are shown for each hosts and tissue niches (column). Likewise, the Page 38/39 standardized d' estimate of preferences for plant species and tissue niches are indicated for each of the observed fungal OTUs (row). Each cell in the matrix indicates a two-dimensional preference (2DP) estimate, which evaluates to what extent each pair of each plant/tissue-fungus association was observed (counts) more or less frequently than would be expected by chance. The P values were adjusted based on the false discovery rate (FDR). RS, Reaumuria soongorica; SR, Sympegma regelii; NS, Nitraria sphaerocarpa; SP, Salsola passerina; EP, Ephedra przewalskii.

Figure 7
Co-occurrence networks of microbial taxa in the fungal communities. Nodes represent fungal OTUs, whereas edges represent signi cant interactive correlations between pairs of OTUs. Node color represents the order of fungal OTUs and the size of nodes corresponds to the relative abundances of speci c fungus. Red edges indicate positive relationships, and green edges indicate negative relationships. The thickness of each edge is proportional to the P values.