Bacterial Diversity and Community Composition Distribution in Cold-Desert Habitats of Qinghai–Tibet Plateau, China

Bacterial communities in cold-desert habitats play an important ecological role. However, the variation in bacterial diversity and community composition of the cold-desert ecosystem in Qinghai–Tibet Plateau remains unknown. To fill this scientific gape, Illumina MiSeq sequencing was performed on 15 soil samples collected from different cold-desert habitats, including human-disturbed, vegetation coverage, desert land, and sand dune. The abundance-based coverage estimator, Shannon, and Chao indices showed that the bacterial diversity and abundance of the cold-desert were high. A significant variation reported in the bacterial diversity and community composition across the study area. Proteobacteria accounted for the largest proportion (12.4–55.7%) of all sequences, followed by Actinobacteria (9.2–39.7%), Bacteroidetes (1.8–21.5%), and Chloroflexi (2.7–12.6%). Furthermore, unclassified genera dominated in human-disturbed habitats. The community profiles of GeErMu, HongLiangHe, and CuoNaHu sites were different and metagenomic biomarkers were higher (22) in CuoNaHu sites. Among the soil physicochemical variables, the total nitrogen and electric conductivity significantly influenced the bacterial community structure. In conclusion, this study provides information regarding variation in diversity and composition of bacterial communities and elucidates the association between bacterial community structures and soil physicochemical variables in cold-desert habitats of Qinghai–Tibet Plateau.


Introduction
The cold-desert ecosystem in the northeastern Qinghai-Tibet Plateau is geographically and ecologically distinct. Despite considered an extreme region compared to the surrounding areas, little is known about its microbial diversity and community composition. In addition to extreme fluctuations in annual temperature, other environmental factors, such as low nutrient availability, low-soil-moisture content, and elevated ultraviolet (UV) radiation make the area extremely arid [1][2][3][4][5]. Cold-desert habitats are proposed to reduce the potential productivity of soil or cause destruction of soil system functioning and support relatively simple ecosystems, including food-web structures comprised of cold-adapted plants and microbial taxa [6]. Soil microbes might be the most abundant and diverse types of organisms on Earth [7,8]. Among others, soil microbial diversity and community functions are greatly affected by complex environmental variations [4,9]. Besides, several studies have shown that anthropogenic activities, especially chemical contamination into the desert ecosystems, affect the diversity and richness of soil microbial communities as well as their physiological activities and functional potential of soil microbes [10,11]. Soil microbial communities might be variant owing to increased anthropogenic activities that occur in the railway track, where herbicide is regularly used to control weed growth for the safety of railway tracks [12]. A better understanding of microbial diversity and community dynamics in cold habitats could help to gain insights into the response and adaptation of bacterial communities to extreme environmental conditions. This evidence is vital to evaluate how natural and anthropogenic drivers affect the plant-microbe-soil system under the physicochemical harsh environment.
Soil microbes play a vital role in nutrient availability and maintaining soil quality. Bacteria constitute a key part of the soil biodiversity and play a crucial role in upholding soil processes, which is critical to sustaining the functioning of terrestrial ecosystems [13,14]. Therefore, it is necessary to study how soil bacterial communities are affected by specific environmental changes or disturbances. Previous work suggested that some bacterial community distribution has been well-defined on a large-scale [15,16]. Moreover, additional studies have suggested that spatial distance mainly influences the distribution of bacterial communities. For instance, over large spatial scales, soil pH is the key factor that shapes the bacterial community composition [17][18][19]. In addition, Zhang et al. [20] reported that salinity might be one of the key determinants of bacterial community differences in desert ecosystems. Similarly, bacterial community structure can be influenced by soil temperature [21] and water content (WC) [22][23][24]. Furthermore, compared with soil's physicochemical properties, microorganisms are often sensitive to fluctuations in the immediate environment and consequently may serve as sensitive indicators of soil health [25,26]. Several microbial indicators have been acknowledged to predict community responses to local environmental habitats. From this perspective, the local scale pattern of soil bacterial communities in different habitats of the cold-desert region of the Qinghai-Tibet Plateau has not been studied until now. In the second stage of the Qinghai-Tibet Plateau railway project from Germu (Golmud) to Lasa, the railway track is 1142 km long. The region contains severe desertification sites, therefore, it was expected that the substantial variation (physiological characteristics, vegetation, and anthropogenic activities) in the cold-desert region.
We hypothesized that the diverse vegetation and anthropogenic activities in the region could drive substantial changes in bacterial diversity and community compositions. To test this hypothesis, three objectives were addressed in the current study: (i) to study the bacterial diversity and community composition of the soil samples collected at a depth of 0-10 cm; (ii) to identify key environmental driving factors and understand their association with bacterial community composition, and (iii) to assess the impact of human activities and vegetation coverage on bacterial community interaction networks.

Study Site Description
The study sites were selected from a cold-desert region in the northwest of XiNing and LaSa, located on the eastern edge of the Qinghai Tibetan Plateau, China, with an elevation of about 2600-3000 m (average elevation of 2800 m). It is a continental desert feature, with a dry and windy climate and flat terrain, a unique alpine, and arid desert landscape. At present, the sandstorm and desertification disasters in the study area are very serious, which restricts the economic and social development of the region. The mean annual temperature is 2.1 • C, with the highest and lowest monthly mean temperature of~9 • C in June and −5 • C in January. The mean annual precipitation is 374 mm. The main soil type in the study area is sierozem, which is very susceptible to wind erosion.

Experimental Design and Sampling Procedure
The Qinghai-Tibet railway project (second phase) from Golmud to Lasa is 1142 km long. The region contains severe desertification sites; therefore, it was expected that the substantial variation (physiological characteristics, vegetation, and anthropogenic activities) in the cold-desert region and hold distinct soil bacterial community composition. Therefore, three sampling sites GeErMu (GEM), HongLiangHe (HLH), and CuoNaHu (CNH) were selected according to anthropogenic activities, vegetation, and desert coverage. Five different habitats were targeted at a specific distance from the railway track as an origin of sampling to ensure a higher level of anthropogenic interference ( Figure 1B). Samples were collected from the railway's track to the sand dune along such transect that had a decreasing trend of anthropogenic activities. Samples were categorized based on habitats: high-vegetation (GH) (70-80%), low-vegetation (GF) (30-50%), desert-land (DL), human-disturb (HD), and sand dune (SD), respectively (Table S1 and Figure 1). In addition, daily based meteorological data were downloaded from the website of China Meteorological Data Service Center (CMDC), i.e., http://data.cma.cn. The meteorological data used in the current study is the average of 30 years (1981-2010) ( Figure S1). The dominant vegetation in each site includes (GEM) Potentilla fruticose, Ptilagrostis dichotoma, Artemisia demissa, Carex sp., (HLH) Leymus secalinus, Calamagrostis epigeio, Leymus secalinus, Saussure sp., Stipa sp., (CNH) Hippophae thibetana, Stipa sp., Saussure sp., Potentilla bifurca, Stipa sp., Hippophae thibetana, Stipa sp., Potentilla bifurca, Saussure sp., and Tripolium sp.  (Table S1) and showing the operational taxonomic units (OTUs) found in each site (circles) and the percentages shared between different sites (lines) ( Figure S2). Shared-OTUs estimates for site pairs are shown. (B) Sampling strategy showing the individual soil sample collected (Table S1). CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu.
In June 2015, at each sampling point, approximately 500 g of soil samples were collected at depth of 0-10 cm. One core was collected from the center of the quadrat (10 m × 10 m) and two soil cores were collected from two corners 8 cm away from the end site of the quadrat and mixed thoroughly. The composite soil was placed in sterile plastic bags and immediately placed in a dry ice box for transportation to the laboratory. The composite soil samples were sieved through a 2.0 mm to remove large gravel. Moreover, the soil samples were air-dried for 1 week and subsequently used for physicochemical characterization, whereas those for bacterial community analysis were immediately stored at −80 • C for later DNA extraction.

Soil Physicochemical Analysis
The soil water content (WC) was determined gravimetrically after drying in an oven at 105 • C for 24 h. Soil pH values were measured after making 1:5 soil: water slurry with a pH meter (Göttingen, PT-10, Sartorius, Germany), and the electric conductivity (EC) of the soil was analyzed using a conductivity meter. Total nitrogen (TN) and total organic carbon (TOC) of soil samples were quantified using an automatic element analyzer (Elementar Vario-EL, Langenselbold, Germany) [27]. The total dissolved solids (TDS) were measured by a conventional conductivity meter.

DNA Extraction, PCR Amplification, and Illumina MiSeq Sequencing
The total DNA was extracted from 0.5 g of each soil sample using the PowerSoil DNA Isolation Kit (MoBio Mo Bio Laboratories, Inc., Carlsbad, CA, USA), following the manufacturer's instructions. To determine the soil bacterial community diversity and composition in each soil sample, the V4-V5 hyper-variable region of bacterial 16S rRNA gene were amplified using the bacterial universal primer set 515F/806R with the barcode [28]. Quantitative PCR reactions were analyzed using the SYBR ® Premix EX TaqTM II kit (Takara-Bio Inc., Mountain View, CA, USA), carried out in 20-µL reactions with 10-µL ddH2O, 0.5-µL (10-µM) forward and reverse primers each, 7.0-µL master mix, and 2.0-µL (2.5 ng,1-µL) of template DNA. The amplification of the 16S rRNA genes consisted of initial denaturation at 95 • C for 3 min, followed by 25 cycles of denaturation at 95 • C for 30 s, annealing of 30 s at 50 • C, and an elongation at 72 • C for 30 s, and final elongation step at 72 • C for 10 min.

Bioinformatics Analysis
Raw FASTQ files were cleaned by removing the adapter and used QIIME (1.7) to performed quality-filtered according to the following criteria: (i) the 300-bp reads were truncated at any site that achieved an average quality score of <20 over a 10-bp sliding window, and the truncated files shorter than 50-bp were removed; (ii) exact barcode matching, two nucleotide mismatch in primer matching, and reads comprising ambiguous characters were discarded; (iii) reads that could not be assembled were removed, only overlapping sequences longer than 10-bp were assembled allowing to their overlapped sequence. Open reference operational taxonomic unit (OTU) picking was done with pick_open_reference_otus.py by the default uclust method dataset (97% similarity cutoff was used, alpha issue, download from the web site of QIIME http://qiime.org/home_static/dataFiles.html), and singletons were discarded throughout OTU picking. The phylogenetic association of each sequence was analyzed by RDP Classifier (http://rdp.cme.msu.edu/) against the dataset using a confidence threshold of 97%. Single rarefaction.py in QIIME was used to generate OTUs table with even reads in each soil sample.
The alpha diversity indices including the Shannon index, Chao richness, and Simpson index were analyzed to compare the diversity and richness of bacterial communities [30]. A Bray Curtis distance tree (at a ≥97% sequence similarity level) of the bacterial community across the different habitats were calculated using PAST with the statistical computing paired group (UPGMA). The relative abundance of the dominant phyla and genera in each sample on the composition of the bacterial community, simulated abundance matrix, and grouping then standardized abundance, ggplot2 draws a stacked chart. Differences be-tween samples were visualized by Non-Metric Multidimensional scaling (NMDS) analysis based on Bray-Curtis distance (metaMDS function, vegan) [31]. Analysis of similarities (ANOSIM) was used to test whether the differences between different groups were significant, using the R package vegan (anosim) [32]. Moreover, an NMDS plot was created using the correlation between responsive variables and bacterial taxonomic composition across different habitats. PERMANOVA (999 permutations) was performed to assess the differences in community composition across the sites.
The characterization of bacterial features differentiated across the different sampling sites at different taxonomic levels was performed using the LDA effect size LEfSe method (LDA > 4.0) for biomarker discovery. To reveal the co-occurrence patterns of bacterial communities based on the relative abundances of genera were retained for network construction. The genera and environmental factors data files were organized as guided in the pipeline. All potential pairwise Spearman's rank correlations between the genera were calculated with "vegan" (psych) packages in R. A cutoff value (adjusted p < 0.05) was incorporated into the network analysis [33]. Network visualization was undertaken with "Gephi" random networks of equal size that were built as actual networks for community comparisons. The Venn diagram was drawn using the "VennDiagram" package [34].

Soil Characteristics of the Sampling Sites
The sampling sites across the different habitats significantly varied in terms of soil physicochemical characteristics (see Table 1). Soil pH was not significantly different in all habitats, ranging from 7.6 to 8.6, which shows a slightly alkaline nature, but higher pH values were recorded in HD across all three sites (GEM, CNH, and HEL). EC in the soils ranged from 26.42 to 89.63 µS/cm, it was significantly higher in G-HD and H-HD than C-HD, a similar trend was found for TDS mg/L. SAL in the soils ranged from 0.02 to 0.06 practical salinity unit (psu) of the soil samples, and it was higher in HD than in other habitats. The highest WC percentage was 11.13, 10.65, and 9.18% for GH. The TN content did not differ across the habitats in either site ranged from 0.03 to 0.06%. The only significant difference was observed in the TN for G-HD. TOC in the soil ranged from 0.31 to 1.06%, it was significantly higher in G-GH and C-HD (Table 1). Generally, the soil characteristics in the G-HD and H-HD were quite different from those in the other habitats.

Bacterial Richness, Diversity, and Taxonomic Composition
The total raw sequencing data were thoroughly processed by mothur and exceeded 384,598 quality-filtered of 16S rRNA sequences across the habitats (five habitats per site). The sequence library of each habitat, after high-quality filtering hold 17,251 to 37,877 sequences ( Table 2). The total abundance, 23,303 OTUs were annotated (identity 97%) in the complete data set, varying from 962 to 2219 across the region. The Shannon and Simpson indices were compared, and bacterial diversity was evaluated across the habitats (4.15-6.58 Shannon; 0.004-0.131 Simpson) ( Table 2). The range of Chao was 1304 to 2767, whereas the range of ACE was 1347 to 2722, which were considered as richness indices. G-DL had a very high level of the Chao index and ACE, consistently, had the leading number of OTUs. Furthermore, H-GH had a low level of Chao and ACE, and the smallest number of OTUs ( Table 2). The OTUs rarefaction curves showed that diversity was completely examined in all samples ( Figure S3 and Table 2). Besides, G-DL habitats had the highest bacterial diversity and also had high community richness. In short, our results indicated that by likening the two index values, we found that the richness index was generally positive to the diversity index. The OTUs were classified into 42 phyla, 99 classes, 219 orders, 420 families, 792 genera. Across all habitats, the explanation for a large proportion of the phyla was as follows, Proteobacteria (12.36-55.72%) was largely dominant phylum, followed by Actinobacteria  Figure 2B). Acinetobacter was the most abundant genus within the G-SD (33.85%) and G-DL (26.68%). Notably, a major portion of the abundant genera was unclassified in the G-HD (61.21%) and H-HD (50.28%), whereas C-HD displays the diverse abundant trend of genera ( Figure 2B). A similar significant difference was detected at the genus level ( Figure S4). Taken together, the results encompassed that soil bacterial community composition in the study sites were high, and there is some unique unclassified bacterial community which presents in the human-disturb habitats.

Distribution of Bacterial Community Across the Habitats
OTUs cluster analysis based on the Bray Curtis similarity index paired group process with algorithm mean was estimated (Figure 3). The dendrogram revealed five clustered groups. The members of the first group encompassed only G-SD; the second group (H-HD) clustered alone and exhibited different bacterial communities. The third group included C-GF and G-HD; the fourth group members consisted of H-DL, G-GF, H-SD, C-HD, C-SD, C-DL, C-GH, and H-GF. G-GH and H-GH have consisted of the fifth group ( Figure 3). Therefore, it was concluded that the H-HD cluster alone and had different bacterial community composition. A higher similarity was detected in the CNH sampling sites excluding C-HD.
Venn diagram revealed that the distribution of OTUs was different across the habitats (Figure 4). In total, of the 8198 OTUs, 7871 OTUs, and 7234 OTUs were detected at GEM, CNH, and HEL sites, respectively, only 451, 606, and 450 OTUs were shared across the five habitats, respectively (Figure 4). In the GEM site, the number of OTUs was unique in each habitat such as 501 (G-DL), 258 (G-SD), 129 (G-GH), 127 (G-HD), and 140 (G-GF), accounting for the total observed OTUs ( Figure 4A). In the CNH site, the number of unique OTUs was highest at C-GF (310) while the lowest were detected at G-GH (92) ( Figure 4B). However, the HLH site recorded the highest number of unique OTUs at H-DL (277) whereas it was lowest at H-GH (53) ( Figure 4C). The number of lowest unique OTUs across the habitats was found at C-GH and H-GH. Therefore, the insertion of the other habitats into the G-HD, C-HD, and H-HD seemed to have affected the community structures, and more OTUs were shared across different communities. These results indicated that the distribution of OTUs responded differently to each habitat, and all the cold-desert natural and anthropogenic habitats had significant impacts on the bacterial community distributions.  Table S3. CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu. A comparison of the bacterial community structure for different habitats based on NMDS analysis revealed that G-HD was separated from the other two C-HD and H-HD habitats. Verifying this, permutational multivariate analysis of difference shown significant variances in bacterial community composition across the habitats (ANOSIM metric, a number indicating the degree of differences between the group R 2 = 0.202, p < 0.05) ( Figure 5A).

Identification of Environmental Variables Shaping Cold-Desert Bacterial Community Structures
To assess the relationship between soil variables and bacterial communities, the constrained analysis of NMDS was performed ( Figure 5B). Overall, the first two axes NMDS1 and NMDS2 explained 69.74% of the variation in the bacterial community structure. The partial NMDS indicated that TN, TOC, and EC were the largest factors of the soil ( Figure 5B). In Figure 5B, the GH and GF samples relatively plotted closed compared to other samples. Across the selected soil variables, TN (R 2 = 0.24, p = 0.005) and EC (R 2 = 0.23, p = 0.006) were important soil variables that correlated with the bacterial community composition in the habitats. However, the other soil variables, including pH, WC, TOC, TDS, and SAL were not significantly correlated with bacterial community composition (Table S2). A co-occurrences network analysis was conducted to assess the correlation between the environmental variables and bacterial communities at the genus level across the different sites (GEM, HLH, and CNH) (Figure 6). At the GEM site, the network analysis had 94 nodes linked by 162 edges, at the HLH site, 101 nodes linked by 157 edges, and the lowest linked consist 93 nodes associated with 102 edges at the CNH site ( Figure 6A-C; Table S4). The number of negative relationships (77.16, 55.51, and 61.78%) was considerably larger than the number of positive relationships (22.84, 44.49, and 38.24%) at GEM, HLH, and CNH, respectively. The average degree (avg K) value of the CNH site was higher than the values of the other two sites. Moreover, the modularity of the HLH and CNH co-occurrence network was higher than 0.4, indicating a more modular structure compared to GEM (<0.4) ( Figure 6A-C; Table S4). Modules 4 accounted for 33.8% of network nodes with TN at the GEM site ( Figure 6A) and at the CNH site (modules 5 accounted for 22.58%). Whereas, at site HLH, modules 3 accounted for 32.67% of network nodes with TOC ( Figure 6B). The environmental variables in the different habitats were highly diverse, therefore, the results indicated that local environmental factors are vital for shaping bacterial community structure on a spatial small scale.  Table S4. CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu.

Indicator Taxa in Soil Bacterial Communities
To assess the significance of discriminatively specific bacterial taxa, the LEfSe analysis was accomplished from the phylum to genus across the three sites using the default logarithmic (LDA) value of 4 ( Figure 7A). The clade graph showed that 22 groups of bacterial taxa were enriched in the HLH site samples, 7 groups in the CNH followed by GEM (3 groups of bacterial taxa), respectively. The community compositions of the bacterial community between the three sites were significantly different, as displayed in the LEfSe cladogram ( Figure 7B). Most of the clades across the sites were unclassified. After removing the unclassified clades, the results showed that different groups at different taxonomic levels can be significantly distinguished across the sites. For instance, Gemmatimonadetes acted as a leading discriminant clade at site CNH. Phycisphaerae is the most differently abundant bacterial class at GEM. At the family level, Cyclobacteriaceae, Microbacteriaceae, Bdellovibrionaceae, Nitrosomonadaceae, Acetobacteraceae, and Rhodobiaceae were significantly enriched at the HLH site, whereas Beijerinckiaceae was dominated at CNH ( Figure 7B).

Discussion
Deserts contribute disproportionately to the terrestrial biodiversity of Earth, especially in the cold-desert soil, where they host hotspots of variation microbial diversity. With about 109-1010 bacterial cells per gram of soil, bacterial diversity varies markedly across the desert ecosystem [35]. Altogether, over the past decade, studies have revealed that certain bacterial taxa are distributed in hostile environmental conditions such as cold-desert. These bacteria adapt special protective strategies to cope with environmental stresses [36]. Comparatively very few reports are present on the bacterial diversity in different habitats across the Qinghai-Tibet Plateau desert. Most of them have focused on the bacterial communities of gas fields, salts, and/or contaminated lakes [37,38]. This study presents a comprehensive assessment of bacterial diversity and structural analyses in different habitats across the Qinghai-Tibet Plateau. Significant variations were reported in physicochemical parameters in the study area. Nevertheless, the harsh characteristics of cold-desert regions vary fundamentally from those of green mountainous regions and are expected to play a key role in maintaining and generating bacterial diversity. With continuous human interference through railway tracks and land use, the role of the desert as refugia for bacterial diversity may come under risk. Our findings consequently demonstrate that different habitats in cold-desert ecosystems harbor significantly diverse bacterial communities.

Variation in Bacterial Diversity and Community Composition
Our findings highlight that bacterial communities had a noticeable difference across the desert habitats ( Figure 5), and that bacterial richness and abundance had a dramatic variation (Table 2). Similar findings have been reported for other dry-land desert ecosystems, showing a strong effect of different habitats aspects on bacterial communities, including karst rocky desert [39] and high arctic desert [40]. In this study, bacterial abundance and richness were higher than those in previous studies in the Badain Jaran desert and Tenger desert [41], and cold desert ecosystem [42], this might be attributed to the methodological variances and region selections. In our study, human-activities and vegetation habitats might affect soil bacterial diversity. Similar findings have been reported in previous studies such as Kuwait desert soils [43] and Drass, cold-desert in western Himalaya [44]. Unexpectedly, these findings suggest that the different habitats from cold-desert in the Qinghai Tibet Plateau comprise more bacterial diversity and the desert environment consists of bacteria having a large reservoir of genetic variability. Since we found 451, 606, and 450 OTUs in common across the five habitats (GEM, CNH, and HLH sites), the five habitats likely comprise different OTUs ( Figure 4A-C), as previously reported in desert ecosystems [45,46]. The variances in bacterial composition across the habitats were significant at the OTUs level, even for the two close locations (HLH-4633 elevation and CNH-4624 elevations). Our results provide evidence that the desert environment comprises numerous different OTUs depends on the habitats rather than the geographical locations.
The bacterial communities were attributed to different phyla, and 12 phyla dominated all samples. This level of diversity recommends a larger ecosystem complexity in the cold-desert habitats than previously appreciated [41,42] (Figure 2A). The bacterial diversity in this study was dominated by desiccation and radiation-resistant bacteria that were reported in previous studies [4,26,43]. The relatively low number of bacterial communities and taxonomic evenness, coupled with the presence of unclassified phyla, show that the HD habitat may characterize a transient soil-borne inoculum rather than the indigenous community. In contrast, the specific existence of Chloroflexi across the habitats, the phylum acknowledged to be important in tundra soils [47] recommends a true soil community might occur. The occurrence of Firmicutes in other cold habitats is well reported in the literature [48]. Previous reports on Himalaya regions, Pradhan et al. [49] and Shivaji et al. [50] stated that Firmicutes as the dominant phylum contrary to our study. Cold-deserts are representative of physicochemical extreme environments, where the scarcity of water is the basic parameter affecting the bacterial community. Therefore, xerophilic bacterial communities that are adapted to comparatively high radiation and low temperatures are likely to be the leading communities in these ecosystems.
The presences of unclassified genera were surprising that indicate the anthropogenic influence on bacterial diversity and the site is rich in novel phylogenetic taxon ( Figure 2B). It has been previously stated that anthropogenic activities are significant sources of contamination [43,51,52]. Recently, the identification of unclassified taxa has motivated a new wave in research. In particular, to elucidate and quantify their contribution to biological processes and influences of essential ecological variations [4,20,26]. Conversely, the specific occurrences of rhizospheric bacteria in vegetative habitats can be attributed to their role in plant promotion in Qinghai Tibet Plateau compared to other desert habitats that lack the existence of rhizospheric bacteria [43,53].
Across the region, the samples (except for C-GF) in CNH are closely grouped compared to the samples in sites GEM and HLH. This grouping shows that the bacterial community structure of cold-desert habitats in a different location has spatial heterogeneity, which is consistent with previous studies [14,26,54]. Interestingly, the samples from GEM and HLH plotted far away and showed higher dissimilarity. Though the study area was narrow, yet the spatial separation, environmental characteristics, and/or species-sorting mechanism [55] are the key defining variables that influenced the bacterial community's composition.

Bacterial Communities in Relation to Environmental Variables
It has been previously reported that environmental factors such as WC, TOC, TN, salinity, and pH are vital in shaping the bacterial community in cold-deserts [44,56,57]. In the present study, environmental variables significantly explained the variation of bacterial communities across the region. Here, we detected that TOC, TN, and EC content were relatively the key variables in desert soil samples ( Table 1), some of these variables have earlier been implicated in correlations with bacterial diversity [26,58]. In this study, pH was significantly positively correlated with Actinobacteria and our results are in line with a previous study [59]. In the previous studies, the bacterial community composition was most affected by TN and TC contents and soil physical variables like WC, but not pH [60], while Ca, C, K, and WC influenced the structure and distribution of microbial communities in the Antarctic valley [56]. Actinobacteria are adjusted to oligotrophic environments where the hyphae allow the bacteria to restore nutrients and moisture through pores in the soil [20,61]. Cyanobacteria are positively correlated with soil variables and play an important role in the N and C economy of desert soil [62]. Of note, SAL has earlier been recommended as an important driver of bacterial community diversity [63,64]. Furthermore, the significant influence of environmental factors was reported on overall OTUs abundance in NMDS analysis. The environmental factors in the study sites differ greatly, particularly TN, TOC, and EC, which recommends that local environmental factors are the key determinants that can determine the bacterial diversity on a small special scale ( Figure 5B), as previously reported in the western Tibetan Plateau [26,65]. Samples GH and GF showed relatively higher similarity that might be due to the common ecological setup provided by vegetation coverage. Similarly, the diverse bacterial community can consume soil nutrients and drive biogeochemical cycling and launch a constant ecosystem in the local niche.

Network Interaction Differences and Potential Indicators
The responses of different network modules were distinct in different sites ( Figure 6). TN showed the maximum number of associations with the bacterial community in GEM and CNH sites while TOC showed maximum interaction in the HLH site. Likewise, TOC might support the predominant genera Arthrobacter in the dune habitat (HLH site) (Figure 2 and Figure S5). This finding might be supported by the Jiang and Yang [66] study, which showed that the north of the desert contains wide oases and organic carbon carried by dust to the desert. Overall, these findings strongly suggest that habitat-filtration may serve as an influential practice defining dune bacterial community assemblies [67]. This result chains the notion that bacterial community assembly is governed by environmentally driven functional features [18] and soil nutrients are crucial in this regard. Regarding the interspecific interactions, we found stronger negative relations than positive between environmental variables and bacterial communities across the study sites. For instance, positive relationships were found between the Actinobacteria and Proteobacteria ( Figure S4) and they are both chemoheterotrophic and aerobic bacteria [68]. Furthermore, different genera of Firmicutes and Actinobacteria can induce dry cold habitats producing endospores or forming spores, respectively; they can also produce antimicrobial substances affecting an extensive spectrum of microorganisms [69]. On the contrary, Fierer et al. [70], reported an increased abundance of antimicrobial-resistance genes in non-desert habitats than desert habitats. This finding signifies that competitive interactions are not as vital in shaping desert bacterial communities. In the notion of these conflicting statements, more emphasis on species interactions, the relative assistances of stochastic and deterministic processes, and how these differ via time and with environmental variables, is required.
Taxa within bacterial communities that quickly respond to habitat variation are considered as potential indicators [71]. Due to local ecological variation, bacterial taxa sensitive to the environmental changes in each site and habitat (Figure 7). Soils from the HLH site indicating that the number of significant differences in bacterial taxa was higher than that of GEM and CNH sites. Some bacteria are retrieved from the vegetative habitat of the HLH site, such as Rhizobium and Sphingomonas, the respective bacteria were also revealed in the plant rhizosphere of Taklimakan desert southern edge [52]. The phylum Gemmatimonadetes are found in a CNH site, which is distributed in a wide range of habitats and contains oligotrophic bacteria with good adaptation to low soil moisture content [72]. In this framework, it is worth stating that the soil moisture content was much lower in the desert and sand dune land (~4-5%, except H-DL habitat) than that in vegetation habitats (~5-11%). The presence of thermophilic bacteria, such as Spirosoma, Adhaeribacter, and Rubellimicrobium might be related to the formation of soil crust [73]. Excitingly, Deinococcus-Thermus is an acknowledged group of extremophilic bacteria that could endure in ionizing radiation environments [74], but in the current study their abundance was not sophisticated. These distinctive potential indicators deliver evidence to well respond to environmental variation occurring in the cold-desert habitats, which is in line with other studies conducted in cold environments [75].

Conclusions
Regardless of the cold and dry nature of the Qinghai-Tibet Plateau desert soils, this study found a relatively high diversity within bacterial community composition. Moreover, this study provides the first insight into bacterial communities associated with humandisturbed, vegetation coverage, desert land, and sand dune habitats across the desert soil. Higher variations in the bacterial communities among the cold-desert desert sites suggest the site-specific variation of different habitats. The bacterial community structure could be best described by TN and EC levels in the soil and much different from other soil physiochemical variables that exist in extremely cold conditions. As a whole, this work provides a comprehensive study that should be used as a foundation for future studies, likewise, unclassified genera were dominated in human-disturbed habitats, for such aspect, many unknown bacterial species appeared to be present and remain to be identified and characterized in future studies.
Supplementary Materials: The following are available online at https://www.mdpi.com/2076-260 7/9/2/262/s1, Figure S1: (A) Temperature data obtained from ordinary Inverse distance weighting (IDW) interpolation based on the 30-year average annual records (1981 to 2015). (B) Precipitation map obtained from ordinary Inverse distance weighting (IDW) interpolation based on the 30-year average annual records (1981 to 2015)., Figure S2: Venn diagram (using relevant subsets of the full OTU data set) presenting the sharing of OUTs between sites. CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu., Figure S3: Rarefaction curves of all sequences of different habitats across 3 sites sampled of cold desert, Qinghai Tibet Plateau, China. CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu., Figure S4: Hierarchical cluster analysis of bacterial communities of different habitats across 3 sites sampled of cold desert, Qinghai Tibet Plateau, China: The color intensity in each panel shows the percentage of genera in a sample, referring to the color key at the right top. The top and the left graph show the cluster results of the bacterial community at the genus-level across sites. CNH-CuoNaHu, HLH-HongLiangHe, GEM-GeErMu. Table S1: The information for the research sites., Table  S2: Results of the perMANOVA analysis of the Bray-Curtis dissimilarities for bacterial phylum community structure in pH, WC, TN, TOC, EC, TDS and SAL. Df-degrees of freedom; SS-sum of squares; MS -mean sum of squares; Pseudo-F-F value by permutation. p-values are based on 9999 permutations., Table S3: Bray-Curtis similarity indices of different habitats across 3 sites of soil sampled of cold desert, Qinghai Tibet Plateau, China., Table S4: Key topological parameters for the co-occurrence network analysis of soil bacterial communities of each site in the cold desert.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.