Participants were recruited from the Center for Reproductive Medicine, Tongji Medical College, Huazhong University of Science and Technology. Inclusion criteria for recovered COVID-19 patients included documented recovery, age from 25 to 45 years, and having detailed medical records during hospitalization and discharge certificate. The diagnosis of COVID-19 was determined by the New Coronavirus Pneumonia Prevention and Control Program (7th edition) published by the National Health Commission of China (14). Exclusion criteria included asymptomatic cases, having received antibiotics or probiotics within 2-months, gastrointestinal diseases, and severe basic diseases. Healthy controls were individuals who were enrolled at regular physical checkups in the same hospital. Ethical approval for the study was obtained by the Ethics Committee for Clinical Research of Reproductive Medicine Center, Tongji Medical College, Huazhong University of Science and Technology. Seven male cases and seven male healthy controls were recruited in this study. All participants came from Hubei province and shared a similar moderate-fat dietary habit (mainly wheat flour, rice, pork, eggs, vegetables), in which fat provided about 30% energy. Meanwhile, participants were instructed to avoid any food and drugs affecting gastrointestinal function before sampling. The collection and storage of fecal samples adhered to the relevant standards. All fecal samples were collected immediately in sterile plastic tubes and stored at −80°C until analysis. Fecal SCFAs levels were detected using the Agilent 7,890B-7,000D GC-MS/MS platform with detailed methods described in the Supplementary Methods. Blood samples were collected from participants for analyzing serum antibody and cytokines with detailed methods described in the Supplementary Methods.

Total bacterial genomic DNA was extracted from fecal samples using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to the manufacturer's protocols. The final DNA concentrations were determined by NanoDrop 2,000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA), Thereafter, the quality of genomic DNA was checked through 1.0% agarose gel electrophoresis. The V3–V4 region of the 16S rRNA gene was amplified using primers 338F (5'-ACT CCT ACG GGA GGC AGC AG-3') and 806R (5'-GGA CTA CHV GGG TWT CTA AT-3') on the GeneAmp 9,700 PCR System (Applied Biosystems, Foster City, CA, US). The PCR reactions were conducted according to the following program: 3 min of denaturation at 95°C, 27 cycles of 30 s at 95°C, 30 s for annealing at 55°C, and 45 s for elongation at 72°C, and a final extension at 72°C for 10 min. PCR reactions were performed in triplicate 20 μl mixture containing 4 μl of 5 × FastPfu Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of each primer (5 μM), 0.4 μl of FastPfu Polymerase and 10 ng of template DNA. Amplification products were extracted from 2% agarose gels and purified by the DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, U.S.) and quantified by QuantiFluor ™ -ST (Promega, USA). Subsequently, the PCR amplification products were pooled in equimolar and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) at the Shanghai Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China).

Raw FASTQ files were demultiplexed, quality-filtered by Trimmomatic, and merged by FLASH (15, 16) with the following criteria: (i) the reads were truncated at any site receiving an average quality score <20 over a 50 bp sliding window; (ii) primers were exactly matched allowing 2 nucleotide mismatching, and reads containing ambiguous bases were removed; (iii) sequences whose overlap longer than 10 bp were merged according to their overlap sequence. And subsequent high-quality reads were clustered into operational taxonomic units (OTUs) at 97% sequence similarity via USEARCH (Version 8.1). The taxonomy of each 16S rRNA gene sequence was classified using the Ribosomal Database Project (RDP) Classifier tool (Version 2.2) against the SILVA 119 16S rRNA database with a confidence threshold of 70%. OTUs were used for alpha diversity (Shannon, Simpson), richness (Chao) with a 97% threshold. Non-metric multidimensional analysis (NMDS) and principal coordinates analysis (PCoA) were performed based on the Bray-Curtis distance matrix calculated using OTU information from each sample. The differentially enriched microbes between two groups from phylum to genus were analyzed via the Wilcoxon rank-sum test. In addition, QIIME2- DADA2 flows were also run for reexamination, and the detailed process as well as results are available in Supplementary Materials. Statistical analyses were conducted in SPSS 26.0 (Chicago, Illinois, USA). Continuous variables were expressed as the mean ± SD and analyzed with the Student t-test. Categorical variables were expressed as percentages and compared by the chi-square test. P < 0.05 were accepted as indicating a significant difference.

We recruited seven male recovered COVID-19 patients, and their clinical characteristics are summarized in Table 1. The mean length of hospital stay was 19.4 days, and the time between discharge and sampling from 78 to 106 days, with an average time of 90 days. Two patients had complications during the course of the disease; of note, gastrointestinal symptoms were reported in one of the patients. According to the New Coronavirus Pneumonia Prevention and Control Program (7th edition), the degree of disease severity was classified as mild, moderate, or severe. The CoV-IgG of all patients was positive and CoV-IgM was found in only one of seven patients, which means that most patients have entered the recovery phase entirely. Moreover, the pro-inflammatory cytokines IL-6 and TNF-α showed no significant differences between recovered patients and healthy people (Supplementary Figure 1). The comparison of the baseline characteristics of COVID-19 group and healthy controls is shown in Supplementary Table 1. There was no significant difference among the groups.

After merging and filtering, 756,359 high-quality sequences were acquired from 14 fecal samples by 16S rDNA gene sequencing with an average sequence number of 54,025 per sample. Subsequently, 397 OTUs were clustered at a 97% similarity level. The rarefaction curves became flat, and the Shannon index also approached a clear plateau, which together suggested that a near-complete diversity has been captured from each sample (Supplementary Figures 2A,B).

The compositions of gut microbiota of study subjects in the level of phylum and genus were shown in Figure 1A. Eight phyla were identified through bacteria microbiota analysis, and both the COVID-19 patients and healthy controls were dominated by four phyla, Firmicutes, Bacteroidota, Proteobacteria, Actinobacteriota, and Fusobacteriota. The Firmicutes were the predominant phylum, contributing 72.0 and 69.8% of the microbiota in healthy volunteers and COVID-19 patients, respectively. Actinobacteriota was the fourth most dominant phyla in healthy volunteers (9.7%), while a clear decrease in COVID-19 patients (1.3%). The second, third, and fifth most dominant phyla in healthy controls were Bacteroidota (13.4%), Proteobacteria (4.0%), and Fusobacteriota (0.8%); however, all of these phyla were elevated in the COVID-19 groups (Bacteroidota 16.1%, Proteobacteria 7.8%, Fusobacteriota 4.1%) (Figure 1A). The dominant genera of gut microbiota in both groups were Blautia, Bacteroides, Agathobacter, Faecalibacterium, and Escherichia-Shigella, with plenty of other sporadic genera (Figure 1B). The composition of the intestinal microbiota of each participant is shown in Supplementary Figure 3. Taken together, the gut microbiota of the two groups remains consistent in the overall microbiota composition but shows differences in proportions of some bacteria.

To investigate whether an infection has a long-term effect on the gut microbiota, we evaluated the richness and diversity of recovered patients and healthy controls using the α-diversity including the Shannon index, Chao index, and Simpson index (Figures 2A–C). Recovered COVID-19 patients showed a decreasing trend in microbial diversity and richness, there was a statistically significant difference among the groups in light of richness estimator, Chao index (p = 0.0298), and diversity estimator Simpson index (p = 0.015). But there was no significant difference in another diversity estimator, the Shannon index (p = 0.055). Other α-diversity estimators are summarized in Supplementary Table 2. Meanwhile, the number of OTUs in the recovered patient group and health group were 291 and 346, and 240 OTUs overlapped between the two groups (Figure 2D). β-diversity was calculated through NMDS and PCoA in the level of OTUs. NMDS indicated that a clear separation between recovered patients and healthy controls (R = 0.02449, p = 0.0250) (Figure 2E). PCoA also produced a similar result (Supplementary Figure 4). Taken together, our results indicated a role of SARS-CoV-2 infection in affecting the fecal microbiota structure.

Many microbes differ substantially between the two groups, characterized by the decrease in beneficial bacteria and increase in opportunistic pathogens (Supplementary Table 3). There was no significant difference in the phylum level. At the class level, recovered COVID-19 patients had a significant decrease in Coriobacteriia (p = 0.030), with an increase in Acidimicrobiia (0.031) compared with controls, partly caused by the abundant change of order Coriobacteriales and Microtrichales (Figure 3A). Additionally, at the order level, the relative abundance of Eubacteriales (0.011) and Micrococcales (0.021) increased in the recovered patients, while Oscillospirales (0.030) decreased compared with the controls (Figure 3B). At the family level, the relative abundance of Ruminococcaceae (p = 0.010) and the dramatic reduction in Coriobacteriaceae (p = 0.037) in recovered patients, with an elevation in the abundance of Eubacteriaceae (p = 0.011), Micrococcaceae (p = 0.015), Microtrichaceae (0.031) (Figure 3C). At the Genus level, recovered COVID-19 patients showed a clear decrease in Faecalibacterium (0.021), Eubacterium hallii group (0.004), Collinsella (0.037), Erysipelotrichaceae UCG-003 (0.001), NK4A214 group (0.031), and a concomitant increase in Flavonifractor (0.041), Eubacterium (0.011), Rothia (0.030), Candidatus Microthrix (0.030), especially the Erysipelatoclostridium (0.001) (Figure 3D). This is partly similar to previous results on COVID-19 inpatients (7). On the one hand, a remarkable depletion was observed in commensals, particularly in the anti-inflammatory bacteria, such as Faecalibacterium and the Eubacterium_hallii_group. On the other hand, opportunistic pathogens such as Rothia and Erysipelatoclostridium showed a higher relative abundance than controls.

To identify potential bacterial taxon biomarkers between recovered patient and healthy controls, linear discriminant analysis (LDA) effect size (LEfSe) was used to find species that differed significantly in abundance between the groups. A cladogram was used to show microbiota structures at the phylum and genus levels by an LDA score >3.5 (Figure 4). The gut microbiome of the COVID-19 group was mainly dominated by Acidimicrobiia, Microtrichales, Candidatus Microthrix, Microtrichaceae, Rothia, Micrococcales, Erysipelatoclostridium, and Micrococcaceae. These species may act as biomarkers to distinguish and analyze the outcomes of recovered patients on gut health.

Family Ruminococcaceae, genus Faecalibacterium and Eubacterium_hallii_group showed a significant decrease in recovered patients and are SCFAs-producing bacteria (17–19). Decreased production of SCFAs has been found to be a consequence in mouse models of influenza A virus infection, which was capable of affecting immune responses in the lungs and modulating disease outcomes (20). Therefore, we examined the fecal concentrations of certain SCFAs, including acetic acid (AA), propionic acid (PA), isobutyric acid (IBA), butyric acid (BA), isovaleric acid (IVA), valeric acid (VA), and hexanoic acid (HA) among the two groups (Figures 5A–G). Intriguingly, there was no significant difference in the level of any SCFA. This result indicated that the SCFA levels of recovered patients were normal, though the upstream bacteria varied. We further analyzed the relationship between SCFA levels and bacterial genera using correlation heatmap based on the spearman correlation coefficient (Figure 5H). Intriguingly, there was no significant association between altered microbes and SCFAs levels. However, other SCFA-producing bacteria, such as Alistipes, Dialister, Butyricimonas, and Phascolarctobacterium, demonstrate remarkable correlation with at least one SCFA, the relative abunance of which was not decreased according to analysis mentioned above.

Taken together, the drop in abundance of several anti-inflammatory bacteria that can secrete SCFAs does not impact SCFA levels, probably because of supplementation from other SCFA-producing bacteria, and other potential impacts of reduced anti-inflammatory bacteria should be paid attention to.