Baihe National Nature Reserve is a habitat for endangered wild animals such as giant pandas (Ailuropoda melanoleuca) and snub-nosed golden monkeys (Rhinopithecus roxellana) in China, but it is not included in China’s Giant Panda National Park, so there are some problems such as lax management and control at present (Liu et al., 2021). The original residents of the nature reserve depended on animal husbandry for their livelihoods, and 19.69% of the total income of the original residents came from animal husbandry, but the feeding and management methods of animal husbandry are extensive (Yuan, 2018). In the nature reserve, the grazing traces even overlap with the home ranges of the Sichuan snub-nosed monkeys. Therefore, our study was conducted in the Baihe National Nature Reserve.

In this study, soil samples were taken from the habitats of Xiapingdi monkey groups with GD and Tongjiashan monkey groups without GD (CK), and the soil samples were selected according to Figure 1. The specific sampling method was as follows: First, survey lines were set up in the area under GD and CK, and the starting and end points of the sample lines referred to the results of Ran et al. (2003). After that, five large quadrats of 20 m × 20 m were set up in each sample line, and three soil samples in each large quadrat were collected randomly; then, the litter layer on the soil surface was removed. After that, a sieve with a pore size of 2 mm was used to remove roots and stones. Lastly, the treated soils were stored at −80°C for DNA extraction.

Frozen soil (2 g) was used to extract DNA using the FastDNA^® SPIN Kit for Soil (MP Biomedicals, United States) according to the manufacturer’s protocols. The DNA concentrations and purity were determined with a NanoDrop 2000 instrument (Thermo Fisher Scientific, Wilmington, DE, United States). The three DNA extracts of the three soils from each large quadrat were mixed for the following sequencing step. The V3-V4 hypervariable regions of 16S rRNA genes were amplified by PCR using the barcoded primer set 338F (5′-ACT​CCT​ACG​GGA​GGC​AGC​A-3′) and 806R (5′-GGA​CTA​CCA​GGG​TAT​CTA​AT-3′). Purified amplicons of 16S rRNA genes were sequenced using the Illumina MiSeq platform (Majorbio Company in Shanghai) with the PE250 (paired-end sequencing 250 × 2) strategy. The data were deposited into the NCBI Sequence Read Archive (SRA) database under accession number PRJNA896787. To compare the sOTU abundance between different samples, the sequencing depth of each sample was normalized to 41723 reads. The taxonomy analyses of sOTUs were conducted using the feature-classifier plugin in QIIME2 with the Greengenes database (v. 13_8) at 97% shared similarity (Xia et al., 2022).

Shotgun metagenomic sequencing of DNA extracted from soil samples was performed simultaneously (Chen et al., 2021). In brief, approximately 1 μg of DNA was used to construct a library with a 300 bp insert size, followed by sequencing using the Zhu et al., 2021 platform (Majorbio Company in Shanghai) with the PE150 (paired end sequencing 150 × 2) strategy. The quality control method used on the raw data is described in detail in (Chen et al., 2021). There were 384443 reads of 10 soil samples. The read data were deposited into the NCBI Sequence Read Archive (SRA) database under accession number PRJNA906922. Following quality control, the clean reads were assembled using IDBAUD (Chen et al., 2021), and the assembled contigs longer than 500 bp were reserved. We obtained 7635863 contigs with an average length of 828.42 bp each from 10 metagenomic samples. The (ORFs) of the contigs were predicted using Prodigal v2.6.3 with a meta model. The coverage of each contig was calculated by mapping clean reads to the contigs using bbmap (https://sourceforge.net/projects/bbmap/) with the default parameters (Ju et al., 2019).

The ORFs were identified against the CARD database (https://card.mcmaster.ca/) using DIAMOND with an E-value ≤10^–10 (Buchfink et al., 2015), and the results that identified ≥80% were identified as ARG-like ORFs (Ma et al., 2015).

To compare coverage between different samples, the coverage of ARG-like ORFs was normalized using the data size of each sample (copies/Gb). The calculation formula was as follows: R e l a t i v e   a b u n d a n c e s = ∑ 1 n N × 150 / L G (1) where n is the total number of ARG-like ORFs in one sample, N is the number of clean reads mapped to ARG-like ORFs, 150 is the length of clean reads, L is the length of the target ARG-like ORFs, and G is the data size (Gb) of clean reads per sample (Ma, et al., 2015; Xiong, et al., 2018).

The ORFs were detected using VFanalyzer on VFDB in setB (http://www.mgc.ac.cn/VFs/) with an E-value ≤10^–10 to identify the virulence factors (VFs), and the results that identified ≥80% were identified as ARG-like ORFs. If ARG-like ORFs also had VFs, ORFs were drug resistance genes that could reflect the risk of drug-resistant pathogens (Liu et al., 2018).

The diversities and relative abundances of the soil microorganisms and their ARGs as well as the abundances of MGEs were their mean values ± standard error in GD or CK. The differences in the diversities and relative abundances of the soil microorganisms and their ARGs as well as the abundances of MGEs between GD and CK were analyzed by a t-test with a significance level of p < 0.05. The average abundances of OTUs were higher than 0.1% in 10 soil samples, and the relative abundances of all types of ARGs were used to calculate the correlation coefficients (ρ) according to Spearman between soil microorganisms and ARGs using the vegan package in R4.2.1. If ρ ≥ 0.8 and p value ≤0.01, the relationship of OTUs and ARGs was used to perform a network analysis using the igraph package in R4.2.1 and then visualized in Gephi (version 0.9.5).

Microbial diversity has a positive effect on the cycling of nutrients, such as the cycling of carbon and nitrogen (Zhang H et al., 2020). Simpson’s index gives more weight to the more abundant species in a sample. It takes into account the number of species present, as well as the abundance of each species, while the ACE diversities shows more species present in a sample (Schmidt et al., 2019). Supplementary Figure S1 shows that the Simpson’s diversities of soil microorganisms in GD were significantly higher than that of the CK (p = 0.032), while the ACE diversities in GD were significantly lower than that of the CK (p = 0.012), which indicated that GD reduced the total soil microbial species, but increased the dominance of dominant microbial species. GD reduced the total soil microbial species, which could be explained by the fact that GD reduces the diversity of plants and the content of humic substances in wildlife habitats (Zhang T et al., 2020), which results in the reduction of mineralization rates of C and N in soil and the decrease of available nutrients for microbial reproduction by GD (Wu et al., 2022). However, the relative abundances of dominant microorganisms were significantly increased by GD, they were Chloroflexi and Verrucomicrobia (p = 0.001 and p = 0.049), while GD had no significant effect on the abundances of other microorganisms (p > 0.05), thus the soil Simpson’s diversities were increased by GD, even resulted in a rise in the function of dominant microbial species in the soil of wild animal habitats. Chloroflexi and Verrucomicrobia are associated with the carbon and nitrogen contents of the soil. Chloroflexi is a photosynthetic bacterium that can synthesize carbon and nitrogen (Singleton et al., 2020), while the abundance of Verrucomicrobia is also significantly positively correlated with the total carbon and nitrogen contents in soils (Freitas et al., 2012). Thus, our results indicated that the cycling nutrients in the soils of the wildlife habitat were reduced by microorganisms under GD (Supplementary Figure S1C).

A total of 650 ORFs containing ARGs were detected in the soil microorganisms, and 291 types of ARGs were classified into 18 kinds of resistance classes in GD and CK (Figure 2A), which indicated that there was a potential risk from ARGs in the habitats of wild animals. Because ARGs in the surroundings of the habitat can affect the health of wild animal (Xu et al., 2020) and the accessibility of ARGs to animals is one of the factors that affect the risk from ARG (Zhang et al., 2022). Further analysis showed there were no significant difference in the total relative abundances of ARGs between the GD and CK (p > 0.05, Supplementary Figure S2A). However, the Simpson indexes of ARGs were increased significantly by GD (p < 0.05, Supplementary Figure S2B). Simpson’s diversity is a measure of diversity which takes into account both species and the relative abundance of the different species (Schmidt et al., 2019), thus combined with our analysis, we found two reasons for the increased the diversity of ARGs: First, there were 287 and 245 types of ARGs in GD and CK. Among them, 46 types of ARGs and triclosan and bicyclomycin were checked in GD but not in CK (Figures 2A, D), which showed that GD could increase types of ARGs. Second, the relative abundances of MLS and mupirocin in GD were higher than those in CK (p < 0.05, Figure 2B). The diversity of ARGs is higher, the risk from ARGs is greater (Yang et al., 2019). Thus, our result indicated that the risk of the ARGs was increased by increasing types and the abundances of ARGs under GD. ARGs in the surroundings of the habitat can affect the health of wild animal (Xu et al., 2020), as well as humans or other animals could promote the transfer of ARGs, and increases the relative abundance of ARGs (Li, et al., 2018; Lu, et al., 2022). Thus, we concluded that the wild animals may have a risk of ARGs in acessibility of ARGs in the habitat, especially, the ARGs of MLS, mupirocin, triclosan and bicyclomycin that only be checked in GD but not in CK. MLS, mupirocin, triclosan and bicyclomycin, can be used to treat infections caused by a variety of sensitive bacteria (Patteson, et al., 2017). We speculated that if the wild animals who lives in the habitats under GD are infected by pathogenic bacteria, the risk of treatment to wild animals will increase, for the pathogenic bacteria have increased their resistance to the four classes ARGs of MLS, mupirocin, triclosan and bicyclomycin under GD (Figure 2).

MGEs play an important role in the horizontal transfer of ARGs, the higher the abundance of MGEs is in microorganisms, the stronger the ability of ARGs to be carried by MGEs, and the greater the risk of horizontal transfer of ARGs (Lu et al., 2022). Our results showed that there were 12160 MGE-like ORFs in both CK and GD, the relative abundances of MGEs in GD were higher than those in CK (Figure 3B; p = 0.000162), which indicated that GD could increase the horizontal transfer ability of MGEs to ARGs. Gut microbiome could acquire exogenous ARGs via horizontal gene transfer (Mei et al., 2021). Table 1 also shows that there are ARG-like ORFs carried by plasmids only be checked in GD, namely, tetM, MexD, ceoB, MuxC, sul1 and catII, which indicated that if the carried ARGs are transferred to gut microbiome of wild animals under GD, the wild animals they will develop resistance to some antibiotics.

The relative abundances of MGEs carried by plasmid, prophage and virus were as follows: plasmid>prophage>virus, and the relative abundances of plasmids in the GD were significantly higher than those in the CK (Figure 3B; p = 0.000160), while there were no significant differences in the relative abundances of MGEs carried by prophages and virus between the CK and GD (Figure 3B; p > 0.05). If both ARGs and MGEs were present in one ORF, the ARGs could be moved by MGEs (Ju et al., 2019). Table 1 shows that 55 ARG-like ORFs could be carried by plasmids, and the total number of ORFs-like plasmids in GD was significantly higher than that in CK (Supplementary Figure S3; p = 0.034). Combined with Figure 3B; Table 1, we found that plasmids had the strongest ability to carry ARGs and GD increased their ability to carry ARGs in wildlife habitats. In nature, ARGs carried by plasmids are the main method by which ARGs are transferred horizontally, studies have shown that the ARGs in the microorganisms of water, activated sludge and soil are mainly carried by plasmids (Che et al., 2019; Fan et al., 2019; Li et al., 2020; Dai et., 2022), which supports our results.

If both an ARG and a virulence factor (VF) are present in one ORF, the bacteria in the ORF are called drug-resistant pathogens that reflected the highest level of risk of ARGs (Zhang et al., 2021). There were 26 ARG-like ORFs that were detected in three types of pathogenic bacteria: Klebsiella pneumoniae subsp. pneumoniae NTUH-K2044, Acinetobacter baumannii ACICU and Neisseria meningitidis MC58. The total number of the 26 types of ORFs in GD was significantly higher than that in CK (p = 0.032, Supplementary Figure S4). Among them, six types of ORFs were detected only in GD but not in CK (Supplementary Table S1), they were found to be two types of pathogenic bacteria: Klebsiella pneumoniae NTUH-K2044 and Acinetobacter baumannii ACICU. Thus, not only the total abundance but also the types of drug-resistant pathogens were increased by GD. Previous studies have shown that the wide application of animal manure could introduce a large number of ARGs and drug-resistant pathogen into the soil surroundings (Duan et al., 2017; Liu et al., 2017; Zhang Y et al., 2020), because antibiotic stress in animal manure under GD can not only increase the abundance and species of ARGs but also increase the abundance of VFs due to the adaptive evolution of pathogenic bacteria (Durrant et al., 2020). So, we speculate that once drug-resistant pathogens are transferred to wild animals, health of the wild animals are endangered. For the acquisition of antibiotic resistance in clinically relevant pathogens diminishes the effectiveness of antibiotics in treating bacterial infections (Zhu et al., 2017).

Although these findings have demonstrated that the risk of ARGs was increased by GD in the habitats of wildlife, as the host, the microbial community determines the pattern and environmental behavior of ARGs, and the transfer of ARGs is also related to microorganisms (Bing et al., 2015). The co-occurrence network correlation between ARGs and bacterial taxa is regarded as an effective way to track potential ARG hosts in various environments (Hu et al., 2017). To investigate the association between ARGs and the bacterial community, a correlation analysis was performed between the relative abundance of ARGs and bacteria. There were 271 points and 520 edges in the CK, of which 56 edges were significantly negative and 462 edges were significantly positive (|r|>0.8, p < 0.01), while there were 305 points and 1,150 edges in the GD, of which 60 edges were significantly negative and 1,090 edges were significantly positive (|r|>0.8, p < 0.01; Figures 4A, B). If there is a significant correlation between the abundance of ARGs and the abundance of microorganisms, the microorganisms can be considered potential hosts of ARGs (Zhu et al., 2018), and the network relationship between ARGs and microorganisms can explain the transfer potential of microorganisms to ARGs in the environment to some extent (Lu et al., 2022). Thus, our results indicated that the bacterial community structure determined the distribution of the resistome and that the horizontal transfer of ARGs was not sufficient to obscure their association with bacterial genomes. The species and numbers of the host bacteria of ARGs were increased by GD, which indicated that the carrying capacity of ARGs in soil microorganisms to move was increased by GD.

Although the network relationship results indicated a significant correlation between ARGs and bacteria in the habitats of wild animals, they could not reflect the ability of microorganisms to carry ARGs. If a microorganism and ARG are on one same ORF, the microorganism can carry the ARGs to transfer (Zeng et al., 2019). Figure 4 shows that there were 232 types of ARG-like ORFs that were detected in seven types of bacteria, namely, Bacteroidetes, Firmicutes, Actinobacteria, Proteobacteria, Spirochaetes, Tenericutes and unclassified bacteria, they carried 5.70%, 17.34%, 8.08%, 54.63%, 23.75%, 23.75% and 13.78% of the ARG-like ORFs in the CK, respectively. While there were 242 types of ARG-like ORFs that were detected in eight types of bacteria, namely, Actinobacteria, Bacteroidetes, Deferribacteres, Firmicutes, Proteobacteria, Spirochaetes, Tenericutes and unclassified bacteria, they carried 6.02%, 4.47%, 1.94%, 14.56%, 63.69%, 1.94%, 1.94% and 10.87% of the ARG-like ORFs in the GD, respectively. Among the ARG-like ORFs, 18 types and 16 antibiotic classes could be carried by bacteria in GD and CK, respectively, and triclosan and bicyclomycin (79 types of ARG-like ORFs were not detected in the CK) were not carried in CK. Results have showed that the bacterial hosts of ARGs were preferably distributed in Proteobacteria (Zeng et al., 2019), which was consistent with our research results that Proteobacteria was the main host of ARGs in the habitats of wild animals. Proteobacteria is one of the largest groups of bacteria in nature, and it includes many pathogenic bacteria (Rizzatti et al., 2017). In our results, the highest proportion (98.73%) of the host of the new ARGs in GD was Proteobacteria, and GD increased the ability of Proteobacteria to carry ARGs (Figures 4A, B). We concluded that the ecological risk of ARGs carried by Proteobacteria in the habitats of wild animals is the largest, which was increased by GD, as supported by the results of Supplementary Table S2 because Supplementary Table S2 shows that all three types of drug-resistant pathogens were Proteobacteria, including the two types of new drug-resistant pathogens in GD.

Figure 5 showed that the Simpson’s diversity of the MRGs was higher than that of the other ARG classes, and GD increased the Simpson’s diversity of the MRGs, as well as the MRGs accounted for the largest proportion in the ARG-like ORFs that were only checked in the GD (Figures 5A, B), which indicated that the risk of MRGs was the greatest, and GD increased this risk. MRGs was the most abundant and types and human disturbance increases the risk of MRGs in soil samples (Mei et al., 2021), which supported our result. Our result showed that the MRGs-like ORFs were the largest proportion that could be transferred by plasmids and they were increased by GD (Figure 5C; p = 0.007), the largest proportions of ARGs that the host bacteria carried were MRGs, GD increased this proportions (Figures 4C, D), as well as the MRGs had the maximum distribution degree in the network, GD increased their degree distribution, (Supplementary Table S2). Our results indicated that GD increased mobility of MRG and the storage capacity of MRGs in multidrug-resistant bacteria, which could be explained by this reason that multidrug-resistant bacteria are mainly generated by transfer of MRGs in MGEs and the accumulating MRGs in microorganisms (Kim et al., 2021). The risk of MRGs favors the development of multidrug-resistant bacteria through horizontal gene transfer (Kang et al., 2022), thus we concluded that GD increased the risk of MRGs by increasing the mobility of MRGs and the storage capacity of microorganisms to MRGs.

MRGs present in microorganisms, which could reflect the risk level of ARGs in a region (Zeng et al., 2019). Today, infection caused by multidrug-resistant bacteria is a major threat to ecosystem health (Fiona and Brion, 2013). The bad thing is that, GD increased the abundance of multidrug-resistant pathogens (Figure 5D; p = 0.033), and the multidrug-resistant pathogens-like ORFs of Klebsiella pneumoniae NTUH-K2044 only check in GD. As early as 2008, result had shown that Acinetobacter baumannii ACICU could cause a variety of bacterial infectious diseases, and it had resistance to some (Iacono et al., 2008). The propensity of this bacterium to rapidly acquire antibiotic resistance could lead to the emergence and spread of Acinetobacter baumannii strains with MRGs, and could develop resistance to latest antibiotics (Law and Tan, 2022), which supported our findings that the ARGs in Acinetobacter baumannii ACICU are MRGs, it could develop resistance to a variety of antibiotics (Supplementary Table S1). More seriously, there was one MGE in the Acinetobacter baumannii ACICU under GD, it was pGMI1000MP, which indicated that the Acinetobacter baumannii ACICU could be moved by plasmid. Acinetobacter baumannii is a pathogen responsible for several serious infections, including pneu-monia, sepsis, and meningitis (Law and Tan, 2022). Thus, we concluded that GD increased the risk of MRGs, multidrug-resistant pathogens may appear in the habitats of wild animals after GD, there is also a risk of horizontal transfer of drug-resistant pathogens to wild animals, which has a profound impact on survival of wild animals.