A total of 36 soil samples, located in the Qinghai Province of China, were selected from 12 field sites under three precipitation gradients (Figure 1). Samples in low precipitation (LP, mean annual precipitation < 200 mm) include LP_32, LP_34, LP_35, and LP_36. Samples in moderate precipitation (MP, mean annual precipitation is between 200 and 400 mm) include MP_2, MP_7, MP_27, and MP_29. Samples in high precipitation (HP, mean annual precipitation > 400 mm) include HP_11, HP_15, HP_17, and HP_22. All samples were from grasslands (Figure 1B) with the detailed information in Supplementary Table 1. Three adjacent samples were randomly taken from each site by inserting the coring device 10 cm into the soil surface. The coring device was cleaned with ethanol (70%, v/v) between sites to avoid cross sample contamination. Each of the soil samples was sieved through a 2 mm sieve to remove rocks and roots, mixed evenly, and placed into a clean zip-lock bag for bioinformatics and physicochemical analysis.

The soil physicochemical characteristics were determined according to established protocols (Ge et al., 2021; Shi et al., 2021). Briefly, soil pH was determined with a soil- to-water ratio of 1:2.5 (w/w) suspension using a pH meter (DELTA-320, China). Soil electrical conductivity (EC) was measured with a soil-to- water ratio of 1:5. SWC was determined by the oven-drying method at 105°C to constant weight. Soil organic matter (SOM) was estimated by the K₂Cr₂O₇ oxidation-reduction colorimetric method. The TC and TN were determined by an elemental analyzer (Vario EL III-Elementar, Germany). NH₄⁺-N and NO₃^–-N were extracted with 1 M KCl and measured by a continuous flow analyzer (SAN + +, Skalar, Holland). The soil’s total phosphorus (TP) was extracted by NaOH solution and determined by Mo-Sb colorimetric method (Wang et al., 2014).

Virus-like particles (VLPs) and microbial abundances in each soil sample were estimated using epifluorescence microscopy (EFM) (Danovaro and Serresi, 2000; Thurber et al., 2009; Han et al., 2017a). 10 ml of 0.22 μm filtered amended 1% potassium citrate (AKC) buffer [1% potassium citrate resuspension buffer amended with 10% phosphate buffered-saline (PBS) and 150 mM magnesium sulfate (MgSO₄)] was added to 3 ± 0.5 g soil (Trubl et al., 2016). Viruses and microbes were physically dispersed via 15 min of shaking at 150 rpm at room temperature. The supernatant-contained viruses and microbes were collected to new tubes. The resuspension steps above were repeated two more times to obtain approximately 30 ml supernatants and then filtered through a 0.45 μm Millex filter to remove large particles. The moderate filtrate was disposed with 5 U/ml DNase I (Thermo Fisher Scientific, Lithuania, European Union) and incubated at 37°C for 1 h to remove free DNA, then filtered through a 0.02 μm Anodisc Al₂O₃ filter membrane (25 mm diameter, Whatman; GE Healthcare, Kent, United Kingdom) supported by a 0.45 μm filter. The dried Anodisc filter membrane was stained with SYBR Green I (Invitrogen, Eugene, OR, United States) working solution (1:400) for 20 min in darkness at room temperature. The stained filter membrane was then mounted on a glass slide with ab antifade solution [50% glycerol, 50% phosphate-buffered saline (0.05 M Na₂HPO₄, 0.85% NaCl; pH 7.5), 0.1% p-phenylenediamine]. Subsequently, ten fields of view were randomly selected for observation under EFM (Nikon, Melville, NY, United States), and the abundances of VLPs and microbial cells and VMR were calculated.

The total DNA was extracted from 0.25 g of soil with the PowerLyser PowerSoil DNA isolation kit (Qiagen, Hilden, Germany), following the manufacturer’s protocol. The DNA extract was fragmented to an average size of about 400 bp using Covaris M220 (Gene Company Limited, Shanghai, China) for paired-end library construction. A paired-end library was constructed using NEXTFLEX Rapid DNA-Seq (Bioo Scientific, Austin, TX, United States). Paired-end sequencing was performed on the Illumina NovaSeq platform (Illumina Inc., San Diego, CA, United States) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using NovaSeq Reagent Kits according to the manufacturer’s instructions.¹ The sequencing depth was 20 Gbp for each sample. The DNA content of LP_34 and LP_36 samples is too little to meet the requirements of library construction.

Metagenome read data are available in the NCBI Short Read Archive under BioProject ID PRJNA782356.

Viruses and microbes in soil were enumerated using EFM. The abundance of virus-like particles (VLPs) ranged from 2.0 × 10⁷ (LP_36) to 1.0 × 10¹⁰ (HP_22) per gram of dry soil, and microbial abundance ranged from 1.0 × 10⁸ (MP_29) to 8.2 × 10⁸ (HP_15) per gram of dry soil. The virus-to-microbe ratio (VMR) also varied in different soils, ranging from 0.11 to 98.3 (Table 1). There were significant differences in VLP abundance among different precipitation gradients. The VLP abundance under HP was the highest, followed by MP and LP. Furthermore, the microbial abundance under HP was significantly higher than that under MP and LP. Interestingly, the VMR under HP was significantly higher than that under LP (Table 1).

Pearson correlations were used to evaluate the factors that might influence the soil VLP abundance, microbial abundance, and VMR. Across all samples, soil VLPs have a significant positive correlation with microbial abundance (r = 0.491, p < 0.01) (Figure 2A). The results showed significant positive correlations of VLPs and microbe counts with SWC (virus, r = 0.868, p < 0.001; microbe, r = 0.487, p < 0.01), total carbon (TC; virus, r = 0.701, p < 0.001; microbe, r = 0.399, p < 0.05), total nitrogen (TN; virus, r = 0.711, p < 0.001; microbe, r = 0.514, p < 0.01), soil organic matter (SOM; virus, r = 0.715, p < 0.001; microbe, r = 0.500, p < 0.01) and TP (virus, r = 0.652, p < 0.001; microbe, r = 0.500, p > 0.05). A significant negative correlation of VLPs was observed with C/N (r = −0.473, p < 0.01), soil pH (r = −0.677, p < 0.001), and EC (r = −0.350, p < 0.05). A significant negative correlation of microbial abundance was observed with C/N (r = −0.434, p < 0.01) and soil pH (r = −0.445, p < 0.01). Particularly, VMR was positively correlated with SWC (r = 0.371, p < 0.05) and VLP abundance (r = 0.484, p < 0.01). Furthermore, the random forest model was used to evaluate the relative importance of soil physical–chemical properties, biological factors, climate factors, and geographical location for predicting VLP abundance, which indicated that SWC was the main environmental factor affecting VLP abundance, followed by TP (p < 0.05) (Figure 2B).

A total of 54 viral orders, 152 families, and 362 genera were identified from the metagenomic data of 30 soil samples. At the family level, the average relative abundance of Siphoviridae (40.1%–54.5%) belonging to the order Caudovirales was the highest in all samples. The family Pandoraviridae (10.1%–13.9%), Myoviridae (7.2%–13.1%), and Herpesviridae (5.0%–8.2%) also accounted for a large proportion (Figure 3A). The relative abundance of Podoviridae (22.6%) in LP_35 was 5–7 times higher than that in other samples (2.7%–3.6%). Finally, other viral families like Poxviridae, Baculoviridae, Autographiviridae, Phycodnaviridae, and Adenoviridae accounting for a small proportion were also detected in each soil. The distribution of dominant viral families was greatly similar among most samples under three precipitation gradients. LEfSe analysis identified five orders and eight families of soil viruses that showed significantly different abundance among three precipitation gradients (Figure 3B). A total of 70.1% belonged to Caudovirales, and most of these belonged to Siphoviridae (50.2%). Only Imitervirales and Miniviridae were significantly enriched in MP. Four orders (Priklausovirales, Bunyavirales, Mononegavirales, and Geplafuvirales) and seven families were significantly enriched in HP.

The predicted hosts including bacteria (primarily) and archaea were identified from 8.44% of vOTUs (47 out of 557 vOTUs) via the PHISDetector tool (Figure 4 ad Supplementary Table 2). Most of the individual links (47) occurred via CRISPRs (black solid line), with three via BLAST (blue dashed line). These hosts spanned 18 genera among four phyla (Bacteria: Actinobacteria, Proteobacteria, and Tenericutes; Archaea: Euryarchaeota). Most viruses were linked to only one host at the genus level; only three viruses were, respectively, linked to two hosts but always within the same order (vOTU_4, vOTU_346, and vOTU_349). Of all the hosts, Rubellimicrobium (involving 21 vOTUs) has the most association with viruses, followed by Streptomyces (involving 6 vOTUs). The overall number of virus–host linkages (pairs) under MP (38) was significantly higher than that under LP (6) and HP (6).

To clarify the viral role in carbon cycling in soils, 557 vOTUs were further annotated for carbohydrate-active enzymes (CAZymes) by the dbCAN server based on the recognition of the CAZyme signature domain. According to results, 23.0% of vOTUs (128 of 557 vOTUs) carried 137 CAZyme genes (Figure 5A). These genes belonged to six CAZyme functional classes, most of which are affiliated to glycoside hydrolases (GHs, 42.3% of 137 viral CAZymes genes), followed by glycosyl-transferase (GTs, 28.5%) and carbohydrate-binding modules (CBMs, 20.4%) (Figure 5A). The most annotated CAZymes genes were found in soil viruses under MP, especially MP_29, followed by LP. However, the soil viruses under HP carried fewer CAZyme genes, only involving CE, GH, and GT (Figure 5B). CBM50 and GH23 were the most widely distributed, with different distributions in five samples (Supplementary Figure 1).

It has been confirmed that the abundance, diversity, and reproductive strategy of viruses in soil are affected by multiple factors, such as soil types, physical and chemical properties, cover plants, and so on (Williamson et al., 2017; Liang et al., 2021). In order to expand our understanding of the distribution of soil viruses across different precipitation gradients, the abundance of microorganisms and viruses in 36 soil samples across three precipitation gradients was quantified by fluorescence microscopy. This method only quantified free and adsorbed virus particles in soils, ignoring temperate viruses, so the true abundance of soil viruses was underestimated. Our data showed that there were significant differences in both VLPs and the microbial abundance among different precipitation gradients; both of them increased with precipitation and reached the highest at HP (Table 1), which was consistent with previous studies that soil tends to have more abundant viruses and bacteria under the condition of high moisture content and organic matter (Srinivasiah et al., 2008). In our data, high carbon and nitrogen nutrients in wet soil are beneficial to the growth of microorganisms (Supplementary Table 3), which may make it more conducive for the encounter of a virus and host, resulting in more viruses. Srinivasiah et al. (2008) found that bacterial abundance increased by 84-fold after the addition of C- and N-rich substrates. Also, VLP counts were highly correlated with multiple soil physical and chemical factors compared with microbial abundance, such as SWC, TN, TC, SOM, TP, pH, EC, and C/N (Figure 2A); the result of the random forest showed that SWC was the main environmental factor affecting VLP abundance (p < 0.05) (Figure 2B). Therefore, we supposed that viruses were more sensitive to environmental factors than microbes and SWC was a key environmental factor driving VLP abundance in soil, which indirectly reflected the impact of climate change (average annual precipitation) on VLP abundance and provided a theoretical basis for the conjecture of Williamson et al. (2017) who strongly suspected that the SWC was a key environmental parameter driving both bacterial and viral abundance in soils.

In this study, the VMR, an index to measure viral activity and study the virus–host interaction in the environment (Ogunseitan et al., 1990; Wommack and Colwell, 2000; Parikka et al., 2017), fluctuated over four orders of magnitude (from 0.11 to 98.3) and was significantly higher in HP than that in LP (Table 1). Since VMR is the ratio of VLPs to microbial abundance, its value has to do with the factors controlling both VLPs and microbial abundance. As previous results showed, both VLPs and microbial abundance have a significant positive correlation with SWC, and SWC mainly contributed to the variability in VLP abundance (Figure 2). In addition, soil VLPs have a significantly positive correlation with microbial abundance (r = 0.491, p < 0.01) (Figure 2A),as has been observed in other soils (Williamson et al., 2005; Liang et al., 2019a,2020). With regard to the correlations between viral abundance and replication strategy and microbial abundance, several hypothesized models have been previously reported. The KtW model proposed that viral predation was density and frequency dependent, while lysogeny is facilitated at low host density (Thingstad and Lignell, 1997). Moreover, a study on coral reef samples has found that both viral abundance and VMR were reduced at high host density and confirmed that lysogeny is more prevalent in an environment with high host densities. The “PtW” model was thus presented (Knowles et al., 2016). Based on our results, the viral lytic efficiency was higher in soils with higher water content, and the dynamic changes between viruses and hosts were more consistent with the KtW model (Thingstad and Lignell, 1997). While the sampling pool of this work was limited, this study indicates the interesting results of viral distribution and replication strategies likely being influenced by SWC and sheds light on the drivers of viral community dynamics in the complicated soil environment. Future effort is needed to reveal the relationship of soil properties and other factors with viral distribution and production and, in turn, the role of viruses in biogeochemical cycling in the soil of the Qinghai-Tibet Plateau.

Taxonomic diversity analysis revealed that Caudovirales (70.1%) was the major viral group in soil. The relative abundance of the Siphoviridae family was the highest in all samples (Figure 3A). This result is consistent with that in southeastern United States agricultural soil and Antarctic soil (Adriaenssens et al., 2017; Liang et al., 2019a). Many previous studies have shown that Siphoviridae was the dominant viral family in both soil and marine systems (Jin et al., 2019; Gao et al., 2020; Zheng et al., 2021). The family Pandoraviridae, Myoviridae, and Herpesviridae also accounted for a large proportion. Pandoravirus, the second largest giant virus after Mimivirus, can infect amoeba. Its DNA genomes can reach 2.5 Mb, much larger than that of other viruses (Philippe et al., 2013). Pandoraviruses were rarely detected in soil habitats; however, a most recent study by Legendre et al. (2018) found that Pandoravirus quercus can be isolated from ground soil in Marseille (France). In our results, Pandoraviridae family accounted for 10.1%–13.9% in each soil sample, which may provide convenient conditions for studying the diversity and evolution of Pandoraviridae family in soil. Importantly, the distribution of dominant viral families (top 10) was greatly similar among most samples under three precipitation gradients, indicating that the change of precipitation had little influence on the composition of soil viruses. While the viral species composition at a family level was similar among samples, different dominant virus populations were significantly enriched in soil with MP and HP as the LEfSe analysis result showed, respectively. For example, the families Imitervirales and Miniviridae were significantly enriched in soil under MP (Figure 3B). However, these results are highly dependent on the viral database that is still much incomplete, which will be further refined, as more environmental viruses are discovered and the database is improved.

In order to examine these viruses’ impacts on the microbial communities and processes, we sought to link them to their hosts via the PHISDetector tool (Zhou et al., 2020). In this study, only 8.44% of vOTUs has been assigned to hosts including bacteria and archaea that spanned 18 genera among four phyla (Bacteria: Actinobacteria, Proteobacteria, and Tenericutes; Archaea: Euryarchaeota) (Figure 4). Consistent with previous reports, the majority of hosts were annotated as bacteria. It was reported that both Proteobacteria and Actinobacteria were the dominant soil bacteria taxa (Liang et al., 2021; Wu et al., 2021; Zheng et al., 2021). In addition, compared with grassland soils in Kansas (Wu et al., 2021) and samples from Lake Michigan (Malki et al., 2015), the host range was narrow because most viruses were linked to only one host at the genus level, which was similar to the findings of Trubl et al. (2018) and ter Horst et al. (2021). Particularly, the overall number of virus–host linkages (pairs) in MP (38) was significantly higher than that in LP (6) and HP (6). The site-specific vOTUs had different predicted hosts, thus illustrating the unique assemblages of soil viruses and hosts across sites with differences in precipitation. Interestingly, only vOTU_467 from LP_35 sample located in the Qaidam Basin was assigned to archaea, probably resulting from the highest salt concentration (Supplementary Table 3). The majority of archaea phages originated in thermophilic or extremely halophilic environments (Mochizuki et al., 2010). The Qaidam Basin was originally a huge lagoon. Due to the continuous elevation of the Qinghai-Tibet Plateau, reduction of rainfall, and evaporation of water, this huge natural salty lake basin was formed. Overall, further studies of viral communities and the virus–host interaction can fuel our knowledge of viral ecology in different soil matrices and how biotic/abiotic factors play a role in structuring viral communities.

Soil is the largest carbon pool in terrestrial ecosystems (Post et al., 1982), and the Tibetan Plateau stores the most abundant soil organic carbon (Stockmann et al., 2015; Liang et al., 2019b). In the past, many studies focused on the influence of microbes on the soil carbon cycle, including bacteria and fungi, but ignored the contribution of viruses. It was confirmed that viruses affect ecosystem carbon processing via the controls of top–down (lysing dominant microbial hosts) and bottom–up (carrying AMGs) (Brum and Sullivan, 2015; Trubl et al., 2018). Virus-encoded diverse AMGs could enhance or expand the host metabolic pathways, thus opening up new ecological niches and affecting biogeochemistry (Roux et al., 2016; Gazitua et al., 2021). For example, photosystem I and II genes were obtained by the marine cyanophages from cyanobacteria, and the expression of these genes during infection promoted the photosynthetic output of host cells (Thompson et al., 2011). Thus, to reveal the potential contribution of viruses to the carbon cycle on the Tibetan Plateau, potential CAZymes in the soil viral genomes were annotated. Finally, 22.2% of vOTUs (128 vOTUs) carrying 59 CAZyme genes (categories) that spanned 137 CAZyme genes were further annotated, with most of CAZymes affiliated to polysaccharide hydrolase activities, implying that viruses may play a potential role in the decomposition of organic carbon on the Qinghai-Tibet Plateau (Figure 5A). However, compared with farmlands (10 CAZymes genes) (Bi et al., 2021) and mangroves (27 CAZyme genes) (Jin et al., 2019), soil viruses in the Qinghai-Tibet Plateau carried more diverse CAZyme genes (species), such as GH5, GH8, CE11, CE14, and so on, indicating the higher occurrence of CAZymes and non-negligible roles of viruses in the soil carbon cycle. This difference might result from the environmental specificity of CAZymes (Bi et al., 2021). Our data suggested that the GHs responsible for the breakdown of complex organic matter were the most abundant, supporting previous findings (Emerson et al., 2018). Further, we found that CAZyme genes were the most abundance in samples under MP (Figure 5B), especially in MP_29, which corresponded to the result of host prediction. However, there were the highest VLPs abundance and the lowest CAZyme genes under HP, probably resulting from virus-carrying AMGs that mainly exist in the lysogenic state. The distribution of integrase genes also supported high abundant lysogenic viruses presented in MP (Supplementary Figure 2 and Supplementary Table 4). In addition, we also identified more specific auxiliary carbohydrate metabolism genes, including GHs, glycosyl transferases, polysaccharide lyases, carbohydrate esterases, and carbohydrate-binding modules (Supplementary Figure 1), implying that more ecological functions of viruses are yet to be discovered. Together, these results suggest that viral infections contribute to soil ecosystem functioning and that further interrogation of soil viral communities will yield a more comprehensive understanding of complex functional networks and ecosystem processes in soil.