Alterations of fecal antibiotic resistome in COVID-19 patients after empirical antibiotic exposure.

As the COVID-19 pandemic spread globally, the consumption of antibiotics increased. However, no studies exist evaluating the effect of antibiotics use on the antibiotic resistance of intestinal flora in COVID-19 patients during the pandemic. To explore this issue, we collected 15 metagenomic data of fecal samples from healthy controls (HCs) with no use history of antibiotics, 23 metagenomic data of fecal samples from COVID-19 patients who received empirical antibiotics [COVID-19 (abx+)], 18 metagenomic data of fecal samples from antibiotics-naïve COVID-19 patients [COVID-19 (abx-)], and six metagenomic data of fecal samples from patients with community-acquired pneumonia [PC (abx+)] from the Sequence Read Archive database. A total of 513 antibiotic-resistant gene (ARG) subtypes of 18 ARG types were found. Antibiotic treatment resulted in a significant increase in the abundance of ARGs in intestinal flora of COVID-19 patients and markedly altered the composition of ARG profiles. Grouped comparisons of pairs of Bray-Curtis dissimilarity values demonstrated that the dissimilarity of the HC versus the COVID-19 (abx+) group was significantly higher than the dissimilarity of the HC versus the COVID-19 (abx-) group. The mexF, mexD, OXA_209, major facilitator superfamily transporter, and EmrB_QacA family major facilitator transporter genes were the discriminative ARG subtypes for the COVID-19 (abx+) group. IS621, qacEdelta, transposase, and ISCR were significantly increased in COVID-19 (abx+) group; they greatly contributed toward explaining variation in the relative abundance of ARG types. Overall, our data provide important insights into the effect of antibiotics use on the antibiotic resistance of COVID-19 patients during the COVID-19 epidemic.

Introduction

Coronavirus disease 2019 (COVID-19, or SARS-CoV-2) has rapidly spread worldwide (Auerswald et al., 2021), and it is exacting substantial medical and economic tolls upon humanity. This disease has been declared as a “Public Health Emergency of International Concern” by the World Health Organization (WHO) and is being treated with various antivirals, antibiotics, and antifungals (Miranda et al., 2020). For example, many COVID-19 patients presenting with mild disease without pneumonia or moderate disease with pneumonia are treated with antibiotics (WHO, 2020). Severely ill COVID-19 patients are treated with a variety of antibiotics to take control of the condition (Fehr and Perlman, 2015). A recent review showed that 72% (1450/2010) of hospitalized COVID-19 patients received antibiotics (Rawson et al., 2020). Current evidence indicates low rates of bacterial/fungal co-infection in patients with COVID-19, yet high prescription rates of broad-spectrum antibiotic agents have been reported (Rawson et al., 2020). In the short term, antibiotics exert a variety of roles during the pandemic, but the long-term use of antibiotics and their misuse will contribute to the development of antibiotic resistance (Miranda et al., 2020). This should be expected a marked increase in antibiotic resistance during this pandemic.

Human intestinal flora are considered an important reservoir of antibiotic-resistant genes (ARGs), which have posed immense threats to public health over the last two decades (Feng et al., 2018; Laxminarayan et al., 2013). The bacteria carrying out antibiotic resistance can persist in human intestines for several years, even after short periods of antibiotic intake (Clemente et al., 2015; Hu et al., 2013). Moreover, it is generally recognized that horizontal gene transfer (HGT) is an important evolutionary force for the spread of the ARGs among bacteria (Forsberg et al., 2014; Goossens et al., 2005). In the COVID-19 pandemic, it is unclear whether the use of empirical antibiotics will have a net-positive or net-negative impact on the emergence and development of antibiotic resistance of intestinal flora. Knowing the impact of the use of empirical antibiotics on the antibiotic resistance of intestinal flora is essential for designing antibiotic stewardship programs during the COVID-19 pandemic.

In the work described herein, we performed a metagenomic analysis on 62 fecal samples obtained from the U.S. National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database. To the best of our knowledge, this is the first study to characterize the alterations of fecal resistome in COVID-19 patients after empirical antibiotic exposure. This study aimed to report the impact of antibiotic use on antibiotic resistance during this pandemic to characterize the broad implications of the pandemic for health and social care systems. This issue must be approached with caution so that the long-term challenges of antibiotic resistance are addressed while preserving access to effective drugs.

Materials and methods

Data collection and information

To explore this issue, we selected and downloaded a total of 62 fecal metagenomic datasets from NCBI BioProject No. PRJNA624223 (Zuo et al., 2020), including 15 metagenomic data of fecal samples from healthy individuals with no history of antibiotic intake in the past 3 months, 23 metagenomic data of fecal samples from COVID-19 patients who received empirical antibiotics and without discontinued use of antibiotics during stool collection, 18 metagenomic data of fecal samples from antibiotics-naïve COVID-19 patients during stool collection, and six metagenomic data of fecal samples from community-acquired pneumonia patients who received empirical antibiotics. See Supplementary Table 1 for detailed accession number information relating to specimens. The basic information of subjects is provided in Supplementary Table 1. FastQ files for all SRA accessions were obtained from the NCBI SRA repository using prefetch and fastq-dump from the SRA toolkit (https://github.com/ncbi/sra-tools).

Bioinformatics analysis

High-quality microbial sequencing data were generated by removing low-quality and short reads and human host sequences using the “metaWRAP read_qc” module (Uritskiy et al., 2018). Metagenomic sequencing yielded 33.1 million microbial reads per sample on average, with data analysis performed using MetaPhlAn3 software, a bioinformatics pipeline developed for taxonomic classification from metagenomic sequence data (Beghini et al., 2021). ARG-like reads were determined in all samples using ARGs-OAP v2.0 (Yin et al., 2018). The default parameters of ARGs-OAP v2.0 were used, i.e., a sequence identity of 80%, an E-value cutoff of 10−7, and an alignment length of more than 25 amino acids. The Structured Antibiotic Resistance Genes (SARG) v2.0 database, provided by the ARGs-OAP v2.0 software developers, was used to search ARGs. The SARG database was created by concatenating the Antibiotic Resistance Genes Database (ARDB; http://ardb.cbcb.umd.edu/) (Liu and Pop, 2009) and the Comprehensive Antibiotic Resistance Database (CARD; https://card.mcmaster.ca/) (Jia et al., 2017), which contained 24 ARG types (e.g., tetracycline-resistant gene) and 1244 ARG subtypes (e.g., tetA and tetB).

ARG-like and 16S rRNA gene sequences were prescreened using customized Perl scripts. We calculated and normalized ARG abundance (copies of ARG per 16S rRNA) using the following equation (Li et al., 2015):where Ni(ARG-like sequence) is the number of the ARG-like sequences aligned to one specific ARG reference sequence; Li(ARGs reference sequence) is the length of the corresponding specific ARG reference sequence; N16S sequence is the number of the 16S rRNA gene sequence identified from the metagenomic reads; L16S sequence is the mean length of 16S rRNA genes sequence (1432 bp) in the Greengenes reference database (http://greengenes.secondgenome.com/downloads/database/13_5) (DeSantis et al., 2006); n is the number of the mapped ARG reference sequence belonging to the ARG type or subtype; Lreads is the sequence length of the metagenomic reads (150 nt).

MGE-like reads were identified and annotated in all samples by an aligned custom MGE database (https://github.com/KatariinaParnanen/MobileGeneticElementDatabase) (Pärnänen et al., 2018). Relatively strict criteria for MGE-like reads annotation with a sequence identity of 80%, an E-value cutoff of 10−7, and an alignment length of more than 25 amino acids were used to reduce the number of false positives. We calculated and normalized MGE abundance (copies of ARG per 16S rRNA) using the following equation:where Ni(MGE-like sequence) is the number of the MGE-like sequence aligned to one specific MGE reference sequence; Li(MGEs reference sequence) is the length of the corresponding specific MGE reference sequence; N16S sequence is the number of the 16S rRNA gene sequence identified from the metagenomic reads; L16S sequence is the mean length of 16S rRNA genes sequence (1432 bp) in the Greengenes reference database (http://greengenes.secondgenome.com/downloads/database/13_5) (DeSantis et al., 2006); n is the number of the mapped MGE reference sequence belonging to the MGE type or subtype; Lreads is the sequence length of the metagenomic reads (150 nt).

Statistical analysis

Statistical analysis and graphical plotting were performed using R software. Alpha and beta diversity, principal coordinates analysis (PcoA), Procrustes analysis, permutational multivariate analysis of variance (PERMANOVA), and partial redundancy analysis (pRDA) were computed using the Vegan R package. We used the Kruskal-Wallis rank-sum test to compare the four groups and applied the Wilcoxon rank-sum test to assess the statistical significance between the two groups. LEfSe was performed to compare the characteristics of different resistance genes from the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+). Spearman's rank correlation between the microbiota and ARG was calculated with Hmisc in R. Network visualization was conducted using the Gephi platform (https://gephi.org/). The correlations of the relative abundance of ARGs and MGEs were evaluated by correlation and best random forest model.

Result

Antibiotic resistome differences among HC group, COVID-19 (abx-) group, COVID -19 (abx+) group, and PC (abx+) group

In total, we detected 18 ARG types and 513 ARG subtypes in the fecal microbiota of the study samples. Tetracycline, multidrug, beta-lactam, macrolide-lincosamide-streptogramin (MLS), aminoglycoside, and bacitracin resistance genes were the main classes of ARGs with high relative abundance [Fig. 1 (a)]. The total relative abundance of ARGs did not show any significant difference between the HC group and COVID-19 (abx-) group. Notably, the relative abundances of total ARGs in the COVID-19 (abx+) group and PC (abx+) group were significantly higher than those in the HC group [Fig. 1(b)]. The average relative abundances of ARGs in the COVID-19 (abx+) group and PC (abx+) group were 4.77 and 5.08 times higher than the HC group, respectively. No significant differences were identified in the number of ARGs between any two of the four groups [Fig. 1(c)].

Fig. 1

(a) Comparison of antibiotic resistance gene (ARG) abundance at type level. Boxplots of (b) total relative abundance and (c) number of ARGs in four groups included in present study.

Beta diversity was analyzed using PCoA by Bray-Curtis, showing that no significant difference only between the COVID-19 (abx+) and PC (abx+) groups was observed [Fig. 2 (a)]. Note that empirical antibiotic exposure modifies not only the abundance of ARGs, but also their composition in the gut microbiome. We compared the within-group variation of the ARG composition based on Bray-Curtis distance and found that the dissimilarity among the subjects in the COVID-19 (abx+) group was statistically higher than those of other groups [Fig. 2(b)]. The magnitudes of the Bray-Curtis distance of four groups fall in the order COVID-19 (abx+)>PC (abx+)>COVID (abx-)>HC. Further grouped comparisons of pairs of Bray-Curtis dissimilarity values demonstrated that the similarity of the HC versus the COVID-19 (abx+) group was significantly lower than the similarity of the HC versus the COVID-19 (abx-) group [Fig. 2(c)]. In addition, the dissimilarity of the HC versus the COVID-19 (abx+) group was significantly higher than the dissimilarity of the HC versus the PC (abx+) group, which revealed that antibiotic treatment has a greater impact on COVID-19 patients compared to CAP patients.

Fig. 2

(a) Beta diversity was analyzed using principal coordinates analysis (PCoA) by Bray Curtis, showing segregation of the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+) group. Beta-diversity boxplots based on Bray-Curtis distances of metagenomics samples from the (b) same group (intra-beta diversity) and (c) between groups (inter-beta diversity). (d) Linear discriminant analysis effect size (LEfSe) histogram used to compare the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+) group. (e) Venn diagram showing shared and unique ARG subtypes amongst four groups. (f) Shared ARGs and abundance comparisons in the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+) group. In the stacked bar plot, the percentage (%) of a specific ARG in one of the four groups is equal to the ratio of its corresponding abundance to the sum of the abundances of the ARG in the four groups.

Next, we compared differential enrichment of ARGs (linear discriminant analysis). The discriminative ARGs were differentially colored for the different groups. The OXA_209, major facilitator superfamily transporter, EmrB_QacA family major facilitator transporter genes, mexF, and mexD were the representative ARG subtypes for the COVID-19 (abx+) group. The mean relative abundance of major facilitator superfamily transporter, EmrB_QacA family major facilitator transporter genes, mexF, and mexD increased by more than 625.79-, 179.32-, 116.89-, and 341.18-fold, respectively, in the COVID-19 (abx+) group than in the HC group. The emergence and dissemination of multi-drug resistant (MDR) bacteria pose a grave public health problem (Nagarajan et al., 2018). In addition, sul2, vanG, tetX2, penA, and VanY genes displayed discriminativity in COVID-19 (abx-) subjects. We found that a greater proportion of vanS, tetM, tet32, catD, aph(3''')-III, tetU, vatB, tetL, vanD, emrE, tet44, vanX, aac(6′)-I, mdtM, adeC, msrC, CMY_62, and ADP_ribosylatingtransferase_arr were detected in the PC (abx+) group [Fig. 2(e)].

Shared ARGs among HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+) group

Fig. 2(d) is a Venn diagram showing shared and unique ARG subtypes amongst all groups. A total of 151 ARGs were common to the four groups despite the total abundance levels. The total number of ARGs in the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group and PC (abx+) was 238, 307, 405, and 358, respectively. The total number of ARGs in the COVID-19 (abx+) group was higher than in the HC group, COVID-19 (abx-) group, and PC (abx+). We found 17 unique ARG subtypes in the HC group and identified 20, 61, and 71 unique ARG subtypes in the COVID-19 (abx-), COVID-19 (abx+), and PC (abx+) groups, respectively. More unique ARGs were identified in the COVID-19 (abx+) and PC (abx+) groups [Fig. 2(d)]. Notably, compared with the HC group, the COVID-19 (abx+) group had 193 unique ARGs, including 156 beta-lactam resistance genes, 16 MDR genes, four aminoglycoside resistance genes, five MLS resistance genes, three trimethoprim resistance genes, two tetracycline resistance genes, two chloramphenicol resistance genes, one bleomycin resistance gene, one vancomycin resistance gene, one rifamycin resistance gene, one polymyxin resistance gene, and one unclassified resistance gene. A stacked column plot further showed detailed information on the types of shared ARGs and their abundance comparison among the HC, COVID-19 (abx-), COVID-19 (abx+), and PC (abx+) groups. Among these shared ARGs, MDR genes were more abundant in the COVID-19 (abx+) and PC (abx+) groups [Fig. 2(f)].

Correlation network of cooccurring ARG subtypes and microbial taxa

Procrustes tests depicted a correlation of 0.804 between ARGs and microbial communities (Procrustes sum of squares, 0.415; P = 0.001; number of permutations, 999) [Fig. 3 (a)]. This result indicates that the ARG profile overall was significantly correlated with microbial communities. The microbial community composition shapes ARG distribution in the gut microbiota. Network analysis showed the detailed co-occurrence patterns between the specific ARG subtypes and microbial taxa. It was indicated that the potential hosts of ARGs could be tracked using non-random co-occurrence patterns in network analysis (Li et al., 2015; Pehrsson et al., 2016; Su et al., 2015). The correlation network consisting of 150 edges and 145 nodes was constructed based on the significant positive correlations (q-value <0.01, r > 0.6) between ARG subtypes and microbial taxa that occurred in at least 20% of samples [Fig. 3(b)]. The detailed cooccurrence patterns between microbial taxa and ARG subtypes are presented in Supplementary Table 2.

Fig. 3

(a) Procrustes analysis of correlations between ARGs and the microbial community. Statistical significance was verified by comparing the similarity of squared distances between matched samples in Procrustes analyses. P-value <0.05 indicated a significant correlation between antibiotic resistance and microbial community composition. (b) Network analysis showing co-occurrence pattern between microbial taxa and ARG subtypes based on Spearman's correlation analysis. An edge represents a strong (Spearman's r > 0.6) and significant (P < 0.01) correlation. The size of each node is proportional to the number of connections. The purple nodes represent microbial species, while the other nodes are colored according to ARG types. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

We speculated that 48 bacterial species were the possible hosts of 97 ARG subtypes conferring resistance to 12 kinds of antibiotics based on the co-occurrence results. Among them, E. coli was predicted to harbor 54 ARGs, 25 of which were MDR genes (TolC, acrA, acrB, acrF, emrA, emrB, emrD, emrK, mdfA, and mdtA et al.), with the balance being one kasugamycin resistance gene (kasugamycin resistance protein ksgA), one aminoglycosides resistance gene (aadA), one sulfonamide resistance gene (sul2), two fosmidomycin resistance genes (rosB and rosA), two tetracycline resistance genes (tet34 and tetA), one polymyxin resistance gene (arnA), two MLS resistance genes (macA and macB), 12 beta-lactam resistance genes (OXA-9, TEM-1, TEM-118, TEM-177, TEM-178, TEM-187, TEM-205, TEM-6, TEM-75, TEM-91, class C beta-lactamase, and metallo-beta-lactamase), and several unclassified resistance genes. Some species (Akkermansia muciniphila, Alistipes indistinctus, Clostridium leptum, Eggerthella lenta, Gordonibacter pamelaeae, Ruthenibacterium lactatiformans, and Streptococcus oralis) were also predicted to harbor four or more ARG subtypes. Other bacterial species were speculated to carry three or one ARG subtypes. For instance, K. pneumoniae carried one ARG subtypes to fosfomycin (fosA), Prevotella copri carried one ARG subtype to tetracycline (tet37), and Erysipelatoclostridium ramosum carried three ARG subtypes to vancomycin (vanH and vanD) and MLS (ermT).

MGEs in fecal microbiome among HC group, COVID-19 (abx-) group, COVID -19 (abx+) group, and PC (abx+) group

According to metagenomic analysis, we detected a total of 19 MGE types among all samples. Transposase had the highest relative abundance, followed by the plasmid. Here, we show the top 15 abundant MGEs and their differences among the four groups. Notably, the COVID-19 (abx+) group had a significantly higher abundance of integrase, IS621, ISCR, plasmid, qacEdelta, and transposase than the HC group. The PC (abx+) group had a significantly higher abundance of Tn916, IS621, plasmid, and transposase compared to the HC group (Fig. 4 ). We evaluated the relative role of MGEs and phylum-level microbial composition to explain the relative abundance of ARGs based on partial redundancy analysis (pRDA). The MGEs and microbial community showed contribution to the ARGs of 17% and 18%, respectively, and their joint effect was 21%. Our pRDA analyses showed that the explanatory proportion of MGEs and microbial community for the total ARG composition was close to equal [Fig. 5 (a)].

Fig. 4

Boxplots showing relative abundances and differences of MGEs among the HC group, COVID-19 (abx-) group, COVID-19 (abx+) group, and PC (abx+) group.

Fig. 5

(a) Total abundance of MGEs significantly correlated with the total abundance of detected ARGs based on Pearson's correlation. (b) pRDA differentiating the effect of microbial communities (at the phylum level) and MGEs on the profile of ARGs. (c) Associations of the relative abundance of ARGs and MGEs evaluated by correlation and best random-forest model. Circle size represents the variable importance (that is, proportion of explained variability estimated with out-of-bag cross-validation). Colors represent Spearman correlations. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

The normalized abundance of MGEs was significantly positively correlated with the normalized abundance of ARGs detected [Fig. 5(b)]. Subsequently, we revealed the correlation between ARG types and MGEs in the fecal microbiome. Transposase, qacEdelta, IS621, IS200, and ISCR contributed more significantly toward explaining variation in the relative abundance of ARG types. Transposase was significantly related to aminoglycoside, beta-lactam, sulfonamide, and trimethoprim resistance genes. QacEdelta was significantly associated with aminoglycoside, bacitracin, beta-lactam, fosmidomycin, kasugamycin, multidrug, sulfonamide, and unclassified resistance genes. IS621 was significantly correlated with fosmidomycin, kasugamycin, multidrug, polymyxin, trimethoprim, and unclassified resistance genes. IS200 was significantly related to beta-lactam, fosmidomycin, kasugamycin, multidrug, polymyxin, trimethoprim, and unclassified resistance genes. ISCR was significantly correlated with beta-lactam, fosmidomycin, kasugamycin, rifamycin, sulfonamide, and trimethoprim resistance genes [Fig. 5(c)].

Discussion

Up to 70% of the patients with COVID-19 receive antibiotic treatment either in the outpatient or inpatient setting (Langford et al., 2021). There are two main reasons that patients with COVID-19 potentially receive antibiotic therapy. First, COVID-19 symptoms may be similar to bacterial pneumonia. Second, patients with COVID-19 may acquire a secondary bacterial infection that requires antibiotic treatment (Knight et al., 2021). There are potential threats during the COVID-19 pandemic that could drive antibiotic resistance. Since the beginning of the COVID-19 pandemic, the dynamics of antibiotic resistance have remained uncertain. Meanwhile, the impact of empirical antibiotic exposure on fecal resistome in COVID-19 patients remains largely unknown.

Consistent with past research, tetracycline resistance genes were found to be dominant in gut flora (Hu et al., 2021; Pal et al., 2016; Feng et al., 2018). Interestingly, the total ARG abundances in the COVID (abx+) group were the highest when compared with the HC and COVID (abx-) groups. This finding confirmed the severe emergence of ARGs in the gut microbiota of COVID-19 patients after empirical antibiotic exposure. In addition to gut microbiota, a similar result was recently reported that ARGs dramatically increased in the oropharyngeal microbiome of COVID-19 patients (Ma et al., 2021). In addition, we found that the total ARG abundances in the PC (abx+) group were also significantly higher than those in the HC group. Furthermore, the results of a recent study indicated that the total relative abundance of detected ARGs was significantly higher in individuals with diseases (non-COVID-19) who received antibiotic therapy for 1 month than that in healthy individuals (Duan et al., 2020). Collectively, these results highlight the need for the rational use of antibiotics in clinical settings.

Although antibiotic resistance is a common phenomenon, antibiotic consumption is a leading cause of the emergence of antibiotic resistance as well as of the acquisition, development, and spread of the resistome (Forslund et al., 2013; Wright, 2007). As with most infections, the initial antibiotic choice in COVID-19 patients is often empirical (Huang et al., 2020). For patients, broad-spectrum empirical antibiotic therapy increases not only antibiotic drug resistance but also the risk of adverse side effects (Llor and Bjerrum, 2014). It is worth stressing that antibiotic-resistant infections pose a major and growing threat to global health, and these infections kill an estimated 700,000 people worldwide annually (O'Neill, 2016). If appropriate measures are not taken, this number is expected to rise to 10 million per year by 2050 (Roope et al., 2019). Our results emphasize the necessity for decreasing unnecessary antibiotic use in the clinical management of COVID-19.

Antibiotic use favors the selection of resistant microbes. Increasing antibiotic resistance in intestinal flora of COVID-19 (abx+) patients can be explained by bystander selection (Tedijanto et al., 2018). Although antibiotic treatment is focused on controlling the pathogenic bacteria that cause co-infection in the lungs of COVID-19 patients, currently available antibiotics usually have antibacterial activity against many bacterial species and disseminate extensively throughout the body (Sullivan et al., 2001). Thus, the bacteria that constitute the human microbiome are subject to the selection pressure exerted by antibiotic consumption (Gustafsson et al., 2003; Lindgren et al., 2009; Nyberg et al., 2007). Hence, we speculate that antibiotic resistance in the microbiota of other body sites may also increase in addition to the gut. This is worthy of further research in the future.

COVID-19 has resulted in an exponential increase in global biocide use, possibly inducing further antibiotic selection pressure and contributing to the selection and development of highly resistant micro-organisms (Getahun et al., 2020; Ruiz, 2021). Compared with the HC group, a vast majority of unique ARGs in the COVID-19 (abx+) group are against beta-lactam antibiotics. A recent study by Peng et al. also found that beta-lactam resistance genes have increased in abundance during the COVID-19 pandemic (Peng et al., 2021). Four MDR genes (mexD, mexF, major facilitator superfamily transporter, and EmrB_QacA family major facilitator transporter) were significantly more abundant among COVID-19 (abx+) patients. MDR transporters, or efflux pumps, play a prominent role in acquired clinical antibiotic resistance (Higgins, 2007; Webber and Piddock, 2003). The emergence and dissemination of MDR bacteria can impair the clinical utility of major antibiotic agents, which could pose serious consequences (Kumar et al., 2013). The worldwide advent of MDR bacteria represents a serious public health issue. Previous research based on metagenomic analysis revealed ARG exchange events between environmental bacteria and clinical pathogens (Forsberg et al., 2012). This implies that receiving antibiotic administration not only increases the COVID-19 patient's own MDR genes of intestinal flora, but may also pollute the surrounding environment. In the future, with the wide spread of MDR bacteria, more clinically available antibiotics may prove to be ineffective.

Within-individual heterogeneity of ARG composition in the COVID-19 (abx+) group was found to be significantly higher compared with the HC, COVID-19 (abx-) groups, and PC (abx+). Meanwhile, the heterogeneity of ARG composition between the HC and COVID-19 (abx+) groups was significantly higher than between the HC and COVID-19 (abx-) groups. Procrustes analysis demonstrated that the microbial community composition shapes ARG distribution in the gut microbiota. We surmise that this heterogeneity found in ARG composition may also exist in intestinal flora. This serious disorder may have a negative impact on subsequent treatment, which is not conducive to the recovery of COVID-19 patients. Compared with antibiotics-naïve COVID-19 patients, the ARG composition of COVID-19 patients receiving antibiotic administration could require more time to return to the level of a healthy population. Whether this effect persists in the long term must be confirmed in investigations with a long-term follow-up.

Infections due to antibiotic-resistant bacteria are a major public health threat. The co-occurrence pattern network implied that 48 species might be potential hosts of 97 ARG subtypes. Among them, E. coli harbored maximum numbers of ARGs, 60% of which were MDR genes. A survey of antimicrobial resistances of fecal Enterobacteriacea in healthy people, short-term hospital patients, and long-term hospital patients also found that E. coli were the most highly representative carriers of resistance in all three groups (Osterblad et al., 2000). These results are similar to those of Feng et al. (2018), but E. coli increased the resistance to aminoglycoside, sulfonamide, and kasugamycin in our results. We posit that this could imply the emergence of E. coli with stronger antibiotics resistance. Considering that E. coli can co-infect COVID-19 patients (Zhu et al., 2020), the prevalence of ARGs is of special concern.

According to previous studies, Streptococcus oralis (Byers et al., 2000), Klebsiella pneumoniae (Martin and Bachman, 2018), Klebsiella quasipneumoniae (Brisse et al., 2014), Klebsiella variicola (Maatallah et al., 2014), Streptococcus salivarius (Olson et al., 2019), Enterococcus avium (Lee et al., 2004), Prevotella copri (Posteraro et al., 2019), Alistipes indistinctus (Parker et al., 2020), Erysipelatoclostridium ramosum (Zakham et al., 2019), Eggerthella lenta (Bo et al., 2020), and Anaerotruncus colihominis (Lau et al., 2006) can be pathogens or potential pathogens that can cause diseases in humans. The fact that these ARG-carrying pathogens are able to resist antibiotics taken by humans represents a high risk to human health.

HGT is a common event in the microbial ecosystem of the human intestinal tract (McInnes et al., 2020). A concerning result is that, compared to the HC group, the abundances of IS621, qacEdelta, transposase, and ISCR were significantly increased in the COVID-19 (abx+) group. These four MGEs also contributed significantly toward explaining variations in the relative abundance of ARG types. For COVID-19 patients, the elevated abundance of these MGEs may increase the opportunity for horizontal gene transfer of ARGs in gut microbiota (Abe et al., 2020; Lu et al., 2018; Larsson et al., 2018). What is more worrying is that HGT in this niche has the potential to mediate ARG transfer from commensal organisms into potential pathogens (Broaders et al., 2013). Consequently, it is possible for E. coli, K. pneumoniae, enterococci, and other opportunistic pathogens inhabiting the human gut to acquire resistance determinants from other members of gut flora (McInnes et al., 2020). ARGs that are carried on mobile genetic elements in pathogens could give rise to an immediate threat to the successful treatment of clinical bacterial infections.

Conclusions

The results of this study indicate that the use of empirical antibiotics has a significantly negative impact on the antibiotic resistance of intestinal flora. To reduce the potential long-lasting effects on antibiotic resistance and promote access to effective antibiotics, the development of treatment guidelines to limit unnecessary antibiotic exposure and of measures to maintain conventional surveillance of antibiotic resistance must be at the forefront of research in the post-COVID-19 era.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.