From 13 to 21 August 2021, E. nutans plants were collected from three distinct grasslands located in Huangyuan, Xinghai, and Chenduo counties of Qinghai Province. The samples were collected at comparable growth stages from each grassland. In Huangyuan County, the elevations of the three sites were 2,567 m, 2,610 m, and 2,662 m. In Xinghai County, the respective altitudes were 3,555 m, 3,634 m, and 3,647 m. Similarly, in Chenduo County, the sites were situated at altitudes of 4,550 m, 4,582 m, and 4,636 m. We randomly selected each grassland sample from four different locations. The gathered material was promptly sent to the nearby laboratory, where it was then cut into 2 cm pieces with a paper cutter. Once thoroughly mixed, approximately 450 ± 10 g of the cut material was placed into a polyethylene bag (38 cm × 50 cm) and sealed under vacuum. A total of 36 bags were made according to three plots times four replicates times three ensiling periods. These bags were then kept at ambient temperature, approximately 25°C, for 60 days. To assess fermentation traits and conduct metagenomic sequencing, silage samples were extracted at intervals of 7, 14, and 60 days during the ensiling process.

DNA was extracted from 20 g of E. nutans silage using a Plant Genome DNA Extraction Kit (DP305; TIANGEN, China) according to the manufacturer’s instructions. The extracted DNA was stored at −80°C before analysis. DNA quality was determined using agarose gel electrophoresis (0.5%). DNA concentration was quantified using a NanoDrop 1,000 spectrophotometer (Thermo Scientific, USA).

The metagenomic sequencing library generation was conducted using a DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s specifications (Wu et al., 2020). After library preparation, metagenomic sequencing was performed on an IlluminaSeq 6,000 platform using the 150 paired-end read strategy. For the sample metagenomic datasets, more than 10.0 Gb of reads were sequenced. The metagenomic sequences were filtered to remove low-quality sequences and adapters to obtain clean data (Chen et al., 2018). The clean data were assembled using MEGAHIT software (−-presets meta-large) (Li et al., 2015). The assembled scaffolds were interrupted at the N-junction to obtain sequence fragments called scaftigs (Nielsen et al., 2014). Filtered scaftigs (>500 bp) were used for statistical analysis and gene prediction. Gene prediction was performed as described by Wu et al. (2020).

The unigenes were compared to the Structured Antibiotic Resistance Genes (SARG) using BLASTN software to search for ARGs (identity >80%, coverage >80%; Yin et al., 2018). Each compared sequence with an identity value greater than the minimum identity value required by the database was assessed to ensure the reliability of resistance gene annotation. The abundance of ARGs in each sample and the resistance mechanism of ARGs were determined after the filtration of comparison results. The ARG and species annotation were conducted based on contigs. Gene coverage was determined using BWA-MEM default cutoff parameters and metagenomic reads mapped to assembled genes (Bogaerts-Márquez et al., 2020).

For MGE annotation, the unigenes obtained from metagenomic sequencing were directly used according to the MGE database (Ma et al., 2013). The unigenes were annotated as plasmid-like, integron-like, transposase-like, insertion sequence transposase (ist)-like, or insertion sequence (IS)-like genes if the threshold similarity for gene annotation was >90%and the length of the unigenes matched the reference sequence by at least 25 amino acids.

After quality control, the contigs were spliced using MEGAHIT (Li et al., 2015). The open reading frame (ORF) was identified using Prodigal (Hyatt et al., 2010), and then the ORFs were made on redundant using CD-HIT to obtain a non-redundant gene set (Fu et al., 2012). The abundance of each gene was calculated using Salmon software (Patro et al., 2017). The non-redundant gene set was compared to the SARG and MGE databases using BLASTX for ARG and MGE annotation, respectively, and then the annotated ARG and MGE contigs were further compared to the NCBI-NR database for taxonomic classification. The data processing and graph drawing were completed using the “dplyr” package and “ggplot 2” package (Wickham et al., 2020), respectively.

To identify clinical resistance genes in E. nutans silage, a BLAST search was conducted between the ARG reference sequence in the SARG database and the Pathosystems Resource Integration Center (PATRIC) database (Snyder et al., 2007). The threshold of comparison was that the consistency and matching length were >80%. The final matched SARG reference sequence was the clinical ARGs. Subsequently, the reference sequence-level gene abundance obtained using the ARG-OAP tool was matched with the identified clinical ARGs to determine the gene abundance of clinical ARGs in the samples (Yin et al., 2018).

Tukey’s test was used to compare the means of ARG types, ARG hosts, and MGE types in R statistical software, and the means were considered significantly different when the p-value was <0.05. The total ARG abundance in each group was visualized using box plots. Principal coordinate analysis (PCoA) was carried out using the “ape” package. The “linkET” package was used to perform a Mantel analysis of the relationship between fermentation quality and ARGs, MGEs, and bacteria (Huang et al., 2021). All additional charts not specifically mentioned were created using the “ggplot 2” package of the R software. Network analysis was performed using BH calibration and Spearman’s correlation. Genes or species with discovery rates <60% were eliminated before correlation analysis, and genes or species with low abundance (gene abundance >5 ppm and bacterial species abundance >0.5%) were removed to lower the number of false-positive results. Gephi version 0.9.6 was used to create network diagrams, which were reserved for results with correlation coefficients >0.8 and p-values <0.5.

The diversity and total abundance of ARGs in E. nutans silage were not significantly affected by altitude (p > 0.05) (Figures 1A,B). After 60 days of ensiling, the number of ARGs at low altitudes was higher than that at medium and high altitudes (p < 0.05) (Figure 1A). The results showed that multiple drugs, aminoglycosides, bacitracin, beta-lactams, and polymyxins were the major ARG types in E. nutans silage (Figure 1C). Similarly, the abundance of ARG types in E. nutans silage was not affected by altitude (p > 0.05). After 60 days of ensiling, the abundance of streptothricin was higher at low altitudes than at medium and high altitudes (p < 0.05). The results of the ARG subtypes showed that bacA, mdtK, acrF, tolC, cpxA, CRP, and acrD were the major ARG subtypes in E. nutans silage (Figure 1D). At the same altitude, the effect of fermentation time on ARG was negligible. The efflux pump was the major ARG resistance mechanism in each group, followed by regulation and antibiotic target alteration (Figure 1E). Altitude and fermentation time had no major effects on ARG resistance in E. nutans silage.

During ensiling, the abundance of Pantoea and Pantoea agglomerans gradually decreased with increasing altitude (p < 0.05) (Figure 2). The abundances of Serratia, Serratia proteamaculans, Pediococcus, and Pediococcus acidilactici were the highest at middle altitudes (p < 0.05). The abundances of Hafnia and Hafnia alvei were the highest at low altitudes (p < 0.05). At all three altitudes, the abundances of Lactiplantibacillus and Lactiplantibacillus plantarum increased with increasing fermentation time.

At low altitudes, Hafnia, Pantoea, and Enterobacter were the major carriers of some dominant ARGs (multidrug, aminoglycoside, bacitracin, beta-lactam, and polymyxin) in Elymus nutans silage (Figure 3A). At mid-altitudes, Serratia, Lelliottia, and Pantoea were the primary carriers of the dominant ARGs (multidrug, aminoglycoside, bacitracin, beta-lactam, and polymyxin) in Elymus nutans silage (Figure 3B). At high altitudes, Hafnia, Serratia, Enterobacter, and Atlantibacter were the main carriers of the dominant ARGs in Elymus nutans silage (Figure 3C). In addition, QANA01 was the only carrier of streptothricin at all three altitudes. Atlantibacter was the only carrier of florfenicol. Staphylococcus was the only carrier of fusidic acid.

As shown in Figure 4A, the MGE-based PCoA showed no remarkable separation among the different altitudes. During ensiling, altitude had no effect on the total abundance of MGEs, except that at 60 days of fermentation, the total abundance of MGEs at high altitudes was higher than that at low and medium altitudes (p < 0.05) (Figure 4B). The tnpA and IS91 were the dominant MGE subtypes in each group (Figure 4C). After ensiling for 7 days, the abundances of IS26 and tnpA1 increased at high altitudes compared to low and medium altitudes (p < 0.05). After 60 days of ensiling, the abundances of tnpAB and tnpA1-IS981 were higher at high altitudes than at low and medium altitudes (p < 0.05). Pantoea was the main carrier of the plasmid (Figure 4D). Enterobacter, Pediococcus, Hornefia, and Lacticaseibacillus were the major transposase carriers. The dominant hosts of ist were Pantoea and Buttiauxella. The dominant IS hosts were Buttiauxella, Yersinia, Enterobacter, Erwinia, and Rahnella. Enterobacter, Pantoea, and Pediococcus were the main carriers of integrases.

At low altitudes, 33 ARGs and 4 MGEs highly co-occurred with 15 bacteria (Figure 5A). Among these visual associations, tnpA2B2-ISLL6 was highly correlated with the majority of the target ARGs. Weissella confusa, Bacterium acnes, and Erwinia persicina were the major bacteria associated with ARGs. These results suggest that tnpA2B2-ISLL6 may be the major mediator that promotes ARG transmission in E. nutans silage at low altitudes. At medium altitudes, 36 ARGs and 5 MGEs showed high co-occurrence with 16 bacteria (Figure 5B). Erwinia persicna, Serratia liquefaciens, and Klebsiella pneumoniae were the major bacteria associated with ARGs. tnpA-IS683, int2, and tnpA27 were MGEs highly correlated with most target ARGs. At high altitudes, 32 ARGs and 8 MGEs highly co-occurred with 16 bacteria (Figure 5C). Lactiplantibacillus plantarum, Pantoea agglomerans, and Enterobacter cloacae were the chief bacteria associated with ARGs. tnpAB was an MGE that correlated with the majority of the target ARGs. Notably, the coexistence patterns of ARGs, MGEs, and bacterial species at high altitudes were more complex than those at low and medium altitudes, and more ARG subtypes were present at high altitudes. This suggests that high-altitude conditions may increase the risk of horizontal gene transfer (HGT) of ARGs in E. nutans silage. Procrustes analysis further showed a close correlation among the three (Figures 5D–F).

Data on the fermentation quality have been presented in a previously published paper (Su et al., 2023). The correlation between fermentation characteristics, ARGs, MGEs, and microbial composition was analyzed using the Mantel test (Figure 6A). Lactate and starch levels were highly correlated with ARG levels (r > 0.4, p < 0.01). Changes in MGEs were highly correlated with lactate, acetate, DM loss, and starch content (r > 0.4, p < 0.01). Bacterial communities were associated with the majority of the fermentation indicators, with pH and lactate levels as the key factors. Figure 6B further shows that there was a positive correlation between the dissimilarity of ARGs and the dissimilarity of environmental factors. The correlation between them increased with the increase in altitude.

Some ARGs in human pathogens have been obtained from environmental bacteria through HGT. Therefore, evaluating the clinical ARGs carried by pathogens in E. nutans silage at different altitudes is important. In the present study, the top 12 clinical ARGs in E. nutans silage collected at different altitudes were screened (Figure 7). The top 12 clinical ARGs included one beta-lactam (penA), one bacitracin (bacA), one kasugamycin (ksgA), two polymyxins (rosA and rosB), and seven multidrug genes (acrA, acrF, CRP, emrD, H-NS, mdfA, and mdtE) in E. nutans silage collected at different altitudes. The bacA and acrF were the predominant clinical ARGs in each group. In addition, the abundance of clinical ARGs was decreased with prolonged fermentation time, in particular, the abundances of acrA, CRP, ksgA, penA, and rosB. Notably, the majority of clinical ARGs (acrA, acrF, CRP, emrD, H-NS, mdfA, mdtE, and penA) with 90% amino acid homology in E. nutans silage have high sequence similarity to ARGs in human pathogens. This indicates that ARGs in E. nutans silage from different altitudes on the QTP pose a concern, with at least some potentially entering the clinical domain.