Molecular Phenomics of Gut Ecosystem in COVID-19 Patients&nbsp;

doi:10.21203/rs.3.rs-208797/v1

Download PDF

Research

Molecular Phenomics of Gut Ecosystem in COVID-19 Patients

https://doi.org/10.21203/rs.3.rs-208797/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Gut ecosystem has profound effects on host physiology and health. Gastrointestinal (GI) symptoms were frequently observed in patients with COVID-19. Compared with other organs, gut antiviral response can result in more complicated immune responses because of the interactions between the gut microbiota and host immunity. However, there are still large knowledge gaps in the impact of COVID-19 on gut molecular profiles and commensal microbiome, hindering our comprehensive understanding of the pathogenesis of SARS-CoV-2 and the treatment of COVID-19.

Results: We performed longitudinal stool multi-omics profiling to systemically investigate the molecular phenomics alterations of gut ecosystem in COVID-19. Gut proteomes of COVID-19 were characterized by disturbed immune, proteolysis and redox homeostasis. The expression and glycosylation of proteins involved in neutrophil degranulation and migration were suppressed, while those of proteases were upregulated. The variable domains of Ig heavy chains were downregulated and the overall glycosylation of IgA heavy chain constant regions, IgGFc-binding protein, and J chain were suppressed with glycan-specific variations. There was a reduction of beneficial gut bacteria and an enrichment of bacteria derived deleterious metabolites potentially associated with multiple types of diseases (such as ethyl glucuronide). The reduction of Ig heave chain variable domains may contribute to the increase of some Bacteroidetes species. Many bacteria ceramide lipids with a C17-sphingoid based were downregulated in COVID-19. In many cases, the gut phenome did not restore two months after symptom onset.

Conclusions: Our study indicates widely disturbed gut molecular profiles which may play a role in the development of symptoms in COVID-19. Our findings also emphasis the need for ongoing investigation of the long-term gut molecular and microbial alterations during COVID-19 recovery process. Considering the gut ecosystem as a potential target could offer a valuable approach in managing the disease.

General Microbiology

Metaproteomics

Metabolomics

Glycoproteomics

Gut microbiome

COVID-19

The ongoing coronavirus disease 2019 (COVID-19) outbreak caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) continues to pose a threat to human beings. In addition to fever and cough, gastrointestinal (GI) symptoms were frequently observed in patients with COVID-19 [1–3]. We have provided the first clinical and molecular evidence for GI infection of SARS-CoV-2 with viral RNA detected in both GI tissues (esophagus, duodenum and rectum) and stool samples [1, 4]. Furthermore, viable virus have been isolated in patient stool samples [5]. Our studies revealed that 11.6% of patients have GI symptoms on admission, and 49.5% of patients developed GI symptoms during hospitalization [1]. Importantly, 11.6% of patients did not exhibit any imaging features of COVID-19 pneumonia but only showed GI symptoms. However, little is known about why and how SARS-CoV-2 affects GI tract. Current evidence suggests that SARS-CoV-2 uses the host receptor angiotensin converting enzyme 2 (ACE2) for cell entry and the transmembrane serine protease 2 (TMPRSS2) for Spike protein (S) priming [6].

Compared with other organs, gut antiviral response can result in more complicated immune responses because of the interactions between the gut microbiota and host immunity. It is known that a normal gut ecosystem plays an important role in maintaining host health and immune homeostasis. Gut microbiota may be associated with the COVID-19 susceptibility, severity, treatment and outcomes in many different ways. Metagenomics has revealed that COVID-19 patients had significant altered fecal microbiome compared with controls [7, 8]. A recent study revealed that severe COVID-19 infection is associated with systemic release of bacterial products [9]. Previous research has demonstrated the gut bacterial translocation to the lung in post-stroke pneumonia due to increased gut barrier permeability [10]. Gut microbiota can also regulate the intestinal infection of certain virus such as norovirus [11]. Furthermore, gut microbiome perturbation can alters immunity to vaccines [12], thus may affect the efficacy of SARS-CoV-2 vaccines.

Currently, there are still large knowledge gaps in the impact of COVID-19 on gut molecular profiles and commensal microbiome, hindering our understanding of the pathogenesis of SARS-CoV-2 and the treatment of COVID-19. Here, we investigate the comprehensive molecular phenomic signatures, including metaproteome, glycoproteome, metabolome, and lipidome, to deepen our comprehension of the gut molecular features related to SARS-CoV-2 infection.

We collected a total of 53 stool samples from 13 COVID-19 patients with a range of one to nine longitudinal time-points that occurred 1–94 days post symptom onset (Fig. 1A, Table S1). Stool samples from 21 healthy subjects severed as controls. Positive RT-PCR results for SARS-COV-2 were observed in stool and/or perianal swab samples even 3 months post symptom onset in a patient with diabetes (patient 1). Furthermore, SARS-CoV-2 viral RNA can persist in stool and/or perianal swab samples long after respiratory samples have tested negative in patients 1, 4, 6, 7, 8, and 9, highlighting the susceptibility of GI tract to SARS-CoV-2 infection.

Multi-omics profiling was performed on each sample to study the alterations of molecular phenomics of the gut ecosystem in COVID-19. Principal component analysis (PCA) showed partially or largely separated multi-omics profiles between COVID-19 patients and controls (Fig. 1B). Quality control (QC) samples of different omics data closely clustered together indicating reproducible MS measurements. Samples from the same patients tended to cluster together indicating greater individual similarity even in the presence of the SARS-COV-2 perturbation.

Disturbed host proteome

Using untargeted metaproteomics, a total of 16279 protein groups (including 268 human proteins) with at least one unique peptide (Table S2.1) and a total of 435632 peptides (Table S2.2) were identified. No SARS-COV-2 protein was detected due to the sensitivity limitation of metaproteomics. Metaproteomics revealed suppressed expression of host proteins involved in immune regulation in COVID-19 (q < 0.05, Fig. 2A, Table S2.3), including IGHV3-64D (immunoglobulin (Ig) heavy variable 3-64D), IGHV3-74 (Ig heavy variable 3–74), and IGLL1 (Ig lambda-like polypeptide 1). Meanwhile, two members of the carcinoembryonic antigen-related cell adhesion molecule (CEACAM) family that belong to the immunoglobulin superfamily (CEACAM5 and CEACAM6), were also down-regulated in COVID-19. Significant reduction of CEACAM 5 has been reported in intestinal epithelial cells (IECs) from inflammatory bowel disease (IBD) patients and this defect correlated with the inability of IBD IECs to activate CD8⁺ T cells [13–15], the main T cell population that can kill virus-infected cells. Hence, the reduced CEACAM 5 in COVID-19 may impair cytotoxic CD8⁺ T cell response in GI tract. Serpin B6, an inhibitor of chymotrypsin-like proteases decreased in COVID-19, while proteases CELA3A (chymotrypsin-like elastase family member 3A), CTRC (Chymotrypsin-C), and MEP1A (Meprin A subunit alpha) increased in COVID-19 (Fig. 2A, Tables S2.3 and S3). Serpin B6 can prevent cathepsin G-dependent neutrophil death [16] and protect dendritic cells from cytotoxic T lymphocyte induced apoptosis [17]. Serpin deficiency may lead to high inflammation via cathepsin G and gasdermin D. In contrast, intestinal alkaline phosphatase (ALPI), which inhibits host inflammatory responses by detoxifying gut bacterial lipopolysaccharide [18], and PLA2G2A (phospholipase A2, membrane associated), which participates in host antimicrobial defense and inflammatory response, were upregulated in COVID-19. These differentiating host proteins did not exhibit significant associations to other typical confounding variables (e.g. age, gender, and medicine, Table S2.3), suggesting they represented COVID-19-associated gut pathologies.

We then investigated the longitudinal changes of these altered proteins in patients 5, 11, 12 and 13, who had serial stools displaying positive to negative stool SARS-CoV-2 infection. Overall, all patients showed considerable protein abundance variations, indicating an unstable gut proteome during the disease course of COVID-19 (Fig. 2B). In many cases, the protein abundance did not restore to the normal levels even several weeks after symptom onset. For patient 5, IGHV3-64D, IGLL1, CEACAM5, and CELA3A showed similar trajectory changes, whose abundance reached a high level 38–40 days after symptom onset but reduced dramatically on day 49. For patient 11, IGHV3-64D and CELA3A elevated on day 37 but decreased sharply on day 40. The longitudinal proteome changes of patient 13 were characterized by a steep fall of IGHV3-64D on days 37 and 39 and a gradually increase of CELA3A, ALPI, and PLA2G2A.

Glycosylation insight into mucosal immunological pathogenesis

To further investigate the phenomics alterations of COVID-19, we studied the protein glycosylation which plays a key role immunological regulation [19] by HILIC based enrichment (Tables S3 and S4). We first analyzed the intact N-glycopeptides using pGlyco 2.0 because this search engine improves the identification accuracy by comprehensive quality control at all three levels of glycans, peptides, and glycopeptides [20]. In total, 4960 glycopeptide-spectrum matches (GPSMs) derived from 54 human proteins were identified with a 1% GPSM FDR (combing peptide FDR and glycan FDR) (Table S4.1). Only 3 microbial N-glycopeptides were identified probably because of the lower abundance of microbial glycoproteins compared with the dominant human glycoproteins (such as IgA) and the incomplete microbial N-glycan database of pGlyco 2.0 (currently only for human and mouse glycans). With the less stringent criteria (1% peptide FDR and 5% glycan FDR), we retrieved significantly more glycosylation features with 8423 GPSMs corresponding to 486 distinct site-specific N-glycans on 177 glycosylation sites from 83 human glycoproteins (Table S4.1). Frequency is important to understand the pathology of the different post-translational modifications (PTMs). The glycosylation frequency (merged from different glycan types (Fig. 3A, Table S4.2) or calculated separately (Fig. 3B, Table S4.3)) of major N-glycosylated sites of proteins involved in neutrophil degranulation (including ANPEP, AZU1, MGAM, CEACAM6, CEACAM8, LCN2, OLFM4, and SERPINA1) and neutrophil migration (GP2) were decreased by up to 81.6% in COVID-19. The N-glycosylation of mucins was dominated by Hex1Fuc1, the frequency of which reduced by approximately 25% (Fig. 3B). The N-glycosylation of proteases was dominated by Hex1Fuc1 and Hex3HexNAc2Fuc1 and the frequency of both glycans was reduced in COVID-19.

In contrast to the above proteins, Ig related proteins including IGHA2, FCGBP, and JCHAIN exhibited greater glycosylation heterogeneity. On the glycosite N131 of IGHA2, the frequency of glycan Hex3HexNAc4 decreased by 63.3% but that of analogue Hex3HexNAc5 (with an additional HexNAc) increased by 46.9% in COVID-19 (Fig. 3B). On the same glycosite, the frequency of glycan Hex3HexNAc3 was comparable between two groups but the analogue Hex3HexNAc4 (with an additional HexNAc) was only detected in COVID-19. On the glycosite N205 of IGHA2, the relative frequency of glycans increased as the number of HexNAc increased. These results suggest the N-glycosylation alterations of gut IGHA2 are characterized by the conjugation of more complex glycans through the attachment of more HexNAc. The glycan specific alteration was also observed in JCHAIN (N71), where glycan Hex3HexNAc3Fuc, with an additional Fuc compare to its counterpart, exhibited higher frequency in COVID-19. On the other hand, the frequency of the same glycan on different sites can be quite different. For instance, while Hex3HexNAc2Fuc was only detected in COVID-19 on N1063 of FCGBP, the frequency of this glycan was decreased in COVID-19 on N1317. Taken together, the overall N-glycosylation of IGHA2, FCGBP, and JCHAIN was suppressed with glycan-specific and site-specific variations.

We also extended our analysis to O-glycosylation and performed intensity based label-free quantification (Tables S4.4 and S4.5). Similar to N-glycoproteome, O-glycoproteome also revealed increased glycosylation of proteases and reduced glycosylation of IGHA2, FCGBP, ANPEP, and GP2. As shown in Fig. 3C, the relative abundance of O-HexNAcylated peptides QQLQS₂₀₅KNECGILADPK from FCGBP, PSTPPTPS₁₁₁PSTPPTPSPSCCHPR from IGHA1, and SVTWSESGQNVT₄₉AR from IGHA2 were significantly decreased in COVID-19 (q < 0.05), while the corresponding protein abundance did not change.

Reduction of beneficial gut bacteria and potential host-bacteria interactions

We used the metaproteomics approach, which is more accurate than sequencing methods for biomass estimates, to investigate microbial community structure and activity [21]. The relative abundance of 34 bacterial taxa were significantly changed between healthy subjects and patients with COVID-19, most of which were from the Firmicutes phylum (20 out of 34, 58.8%) followed by the Bacteroidetes phylum (10 out of 34, 29.4%) (Table 5). Strikingly, the relative abundance of all 20 altered members in the Firmicutes phylum significantly decreased in COVID-19 (p < 0.05), the majority of which were butyrate-producers [22] belonging to the Lachnospiraceae family, such as genera Lachnoclostridium, Ruminococcus, Butyrivibrio, and Dorea, and species Blautia hansenii, Ruminococcus lactaris, and Tyzzerella nexilis (Fig. 4A). There was also a significant depletion of butyrate-producing genus Eubacterium in COVID-19, which also carry out bile acid and cholesterol transformations in the gut, contributing to gut and hepatic homeostasis through modulation of bile acid metabolism [23]. In addition, a recent study has found that several species of the phylum Firmicutes (such as genera Clostridium, Ruminococcus, and Eubacterium) were positively associated with memory scores, while species from the phylum Bacteroidetes mainly presented negative associations with memory scores [24]. Taken together, these data suggest a significant reduction of beneficial gut bacteria in COVID-19.

The relative abundance of all altered members in the Bacteroidetes phylum significantly increased in COVID-19 (p < 0.05), such as Bacteroides coprophilus, Bacteroides coprocola, Bacteroides graminisolvens, Bacteroides uniformis, and Bacteroides stercoris (Fig. 4A). Importantly, it has been shown Bacteroidetes and Firmicutes bacteria mainly down-regulate and up-regulate ACE2 expression in the murine gut, respectively [25]. Therefore, the enrichment of Bacteroidetes and the reduction of Firmicutes may potentially inhibit SARS-CoV-2 entry by down-regulating intestinal ACE2 expression.

Association analysis of altered host proteins and bacteria revealed potential host-microbiome interactions. Overall, bacteria groups increased in COVID-19 including B. coprophilus and B. coprocola exhibited negative correlations with host proteins, while those increased in COVID-19 such as Ruminococcus and Fusobacteria exhibited positive correlations (Fig. 4B). An exception was CEACAM6, which was positively associated with B. coprophilus. A recent study has shown CEACAM6 is critical for pathogen enterotoxigenic Escherichia coli adhesion [26]. The reduction of host proteins such as IGHV3-73 and IGHV3-64D may potentially contribute to the enrichment of Bacteroidetes phylum because of the reduced anti-bacteria Igs.

Functional alteration of gut microbiome

Gene ontology (GO) analysis of metaproteomics data revealed that 9 biological processes of microbial proteome exhibited significant difference between healthy subjects and patients with COVID-19 (q < 0.05) (Fig. 4C). Among them, CTP biosynthetic process, GTP biosynthetic process, and UTP biosynthetic process reduced in COVID-19, while de novo' AMP biosynthetic process increased in COVID-19. Untargeted metabolomics revealed that nucleobase (guanine), nucleosides (adenosine, guanosine, 2'-deoxyadenosine, and inosine) and nucleotides (adenosine 5'-monophosphate (AMP), thymidine 5'-monophosphate (TMP), 2'-deoxyguanosine 5'-monophosphate (dGMP)) decreased in COVID-19 (q < 0.05) (Fig. 5A), while cyclic AMP (cAMP), methylated purines (1-methyladenine, 6-dimethyladenine) and methylated pyrimidine (5-methylcytosine, 1,3-dimethylxanthine) increased in COVID-19. Association analysis of microbial and metabolomics data revealed that adenosine was positively associated with class Clostridia and order Clostridiales, and guanine and guanosine were positively associated with genus Butyrivibrio (Fig. 4C). In contrast, 1-methyladenine was negatively correlated with genus Dorea, order Clostridiales and class Clostridia, and 5-methylcytosine and 6-dimethyladenine were negatively correlated with genus Ruminococcus. Consistent with the metabolomics findings, GO analysis indicated that there was a 1.8-fold increase in the protein abundance of DNA methylation process in the COVID-19 group, although this difference only reached a relaxed statistical significance threshold (raw p = 0.02). On the other hand, the process of tRNA aminoacylation (lysyl-tRNA aminoacylation, isoleucyl-tRNA aminoacylation) an essential step of protein synthesis, increased in COVID-19.

Enrichment of bacterial related deleterious metabolites

Using untargeted metabolomics, we identified 96 fecal metabolites significantly differed between control subjects and COVID-19 patients, mainly including nucleosides, nucleotides, bile acids, carboxylic acids, dipeptides, tripeptides, and acylated amino acids (Table S6). Notably, we detected an enrichment of several gut microbiome-related deleterious metabolites in COVID-19 (Figs. 5B and S1), including phenylacetyl glutamine (q = 0.01), which promotes cardiovascular disease such as platelet thrombosis [27], and salsolinol (q = 0.003), which is a potential gut bacterial neurotoxin contributing to the development of neurodegenerative diseases [28, 29]. The reduction of Firmicutes phylum (such as class Clostridia, order Clostridiales, and genus Dorea) may be at least partially responsible for the increment of phenylacetyl glutamine because they were correlated inversely with each other (Figs. 5C). Longitudinal analysis indicated that the phenylacetyl glutamine level was sustained at high levels in sever patient 5 and in patient 11 who exhibited significant GI symptoms two months after symptom onset. In contrast, this metabolite was kept at a steady and normal level in patient 12 throughout the course of disease and restored to a normal level in patient 13 after one month following symptom onset (Fig. 5B).

We also observed elevated levels of uric acid (q = 0.002) in COVID-19, a uremic toxin playing an important role in several kidney diseases such as lithiasis, gout nephropathy, and preeclampsia. One third of endogenous uric acid is extrarenally excreted via the gut lumen, where it undergoes uricolysis by gut microbiota [30, 31]. Increased fecal uric acid was positively associated with several Bacteroides species (Figs. 5C). Interestingly, although all COVID-19 patients involved in this study were non-drinkers, a significantly higher abundance of ethyl glucuronide in COVID-19 (q = 0.01), a metabolite of ethanol formed by glucuronidation, was observed in the COVID-19 group, which indicates a higher susceptibility of ethanol toxicity. Recent studies have demonstrated that certain gut bacteria (such as Klebsiella pneumoniae) contribute to endogenous ethanol production and promote the development of non-alcoholic fatty liver disease [32–34]. Furthermore, gut microbial (such as E. coli and Clostridum sordellii) β-glucuronidases could hydrolyze ethyl glucuronide, which may increase the retention of ethanol in the body by enterohepatic circulation [35]. For patient 11, both ethyl glucuronide and uric acid climbed sharply on day 35 of disease onset, when the discriminative proteins IGHV3-64D, CELA3A, ALPI, and PLA2G2A also dramatically increased (Fig. 5B).

Bile acids are critical for lipid absorption, antibacterial defense and immune regulation [36]. Gut microbiome mediates the primary-to-secondary bile acid conversion. Primary bile acids (chenodeoxycholic acid and muricholic acid), two glycine conjugates (glycochenodeoxycholic acid and glycocholic acid), and secondary bile acids (ursodeoxycholic acid and hyodeoxycholic acid) were decreased (q < 0.05) in fecal samples from participants with COVID-19, compared with control samples (Figs. 5A). Furthermore, a newly discovered conjugated bile acid phenylalanocholic acid³⁷ was also decreased in COVID-19. Hyodeoxycholic acid exhibited significant positive associations with many taxonomic groups including class Clostridia, order Clostridiales, family Eubacteriaceae, and genera Butyrivibrio, Dorea, and Eubacterium (Figs. 5C).

Alerted microbial lipidome profiles

A total of 4,124 lipid features covering 5 lipid categories (sphingolipid, phospholipid, neutral lipid, glycoglycerolipid, fatty acyl and other lipid subclasses) and 67 lipid subclasses (Table S7) were identified based on diagnostic fragment ions along with associated acyl chain fragment information. The most commonly identified lipid species in the fecal lipidome belonged to the ceramide (Cer) subclass with 923 identifications, followed by the triacylglycerol (TG) and monohexosylceramides (Hex1Cer) subclasses with 467 and 349 identifications, respectively (Fig. 6A). Other frequently identified lipid species included the diradylglycerol (DG, 265 identifications), phosphatidylcholine (PC, 245), monogalactosyldiacylglycerol (MGDG, 174), dihexosylceramide (Hex2Cer, 166), phosphatidylethanolamine (PE, 154), OAcyl-(gamma-hydroxy) fatty acid (OAHFA, 135), and sphingomyelin (SM, 122) subclasses. Among the top 30 identified lipid species, Hex1Cer, SPH, and Cer, all of which belong to the sphingolipid category, underwent the greatest amount of change, with 24.9%, 17.2%, 15.9% significantly increased (q < 0.05) in the COVID-19 group compared to the control group, whereas only 0.9%, 3.5%, and 3.3% significantly decreased (q < 0.05) in the same comparison, respectively (Fig. 6A, Table S7). Within the phospholipid category, the proportions of upregulated lipids were much greater than those of downregulated lipids for the PC (5.7 vs. 1.6%), PE (9.1 vs. 0.7%), and cardiolipin (CL, 6.4 vs. 0.9%) species, while the proportions of upregulated lipids were lower than those of downregulated lipids for lysophosphatidylglycerol (LPG, 0 vs. 11.4%), lysophosphatidylserine (LPS, 0 vs. 11.1%), lysophosphatidylethanolamine (LPE, 0 vs. 3.9%), phosphatidylglycerol (PG, 2.9 vs. 4.9%). Within the neutral lipid category, the proportions of upregulated lipids were comparable to or lower than those of downregulated lipids for the DG (2.6 vs. 6.8%) and TG (5.8 vs. 4.5%) species.

Gut bacterial sphingolipids like Cer, although less well characterized than their mammalian counterparts, are increasingly understood to play important roles in microbial-host interactions [38–40]. The sphingoid backbones and attached fatty acyl chains of bacterial sphingolipids are often odd-chain length, hydroxylated or methylated, while the sphingoid bases for mammals are predominantly even chained and linear backbones⁴¹. We found many Cer lipids downregulated in the COVID-19 group have a C17-sphingoid base, which probably derived from gut bacteria (based on the odd number of carbon atoms). Specifically, five C17-Cer lipids with trihydroxy sphingoid bases, including Cer(t17:0/17:0 + O), Cer(t17:0/23:0 + O), Cer(t17:0/24:0 + O), Cer (t17:1/16:0), and Cer (t17:1/23:0 + O) were significantly reduced in COVID-19 (q < 0.05), but no C17-Cer lipids with trihydroxy bases were increased in COVID-19 (Fig. 6B). In addition, a total of 6 Cer lipids with monohydroxy sphingoid bases, including Cer(m17:1/24:1), Cer(m17:1/20:0), Cer(m17:1/26:0), Cer(m17:1/15:0 + O), and 2 Cer(m17:1/16:0 + O) isomers, were significantly reduced in COVID-19 (q < 0.05), all of which have a C17-sphingoid base (Fig. 6C). Unlike C17-Cer species with monohydroxy or trihydroxy sphingoid bases, C17-Cer species with dihydroxy sphingoid bases did not shown significant difference in COVID-19. The chain-length-dependent alteration was also observed for fatty acyl lipid (acyl carnitine (AcCa)) and neutral lipid (acylGlcSitosterol, acylGlcStigmasterol, acylGlcCampesterol, (AcHexSiE, AcHexStE, AcHexZyE, AcHexCmE) subclasses, for which a total of 7 lipids were significantly upregulated in COVID-19 (Fig. 6D). All of them have a C18 acyl chain regardless of the number of double bonds (Fig. 6F).

In addition to chain length, the degree of unsaturation also influenced the behavior of certain lipid species. Highly unsaturated Cer lipids with trihydroxy bases carrying 6 or 7 double bonds (Cer(t40:7) and Cer(t40:6)) were downregulated in COVID-19, while those upregulated Cer lipids have no more than 3 double bonds (Fig. 6B). Similarly, many highly unsaturated DG lipids carrying 5–8 double bonds, such as DG(29:8), DG(28:7), DG(38:6), DG(27:5), and DG(32:5) were downregulated in COVID-19, while those upregulated DG lipids in COVID-19 only have 2 or 3 double bonds (Fig. 6E).

Increased lipid peroxidation and disturbed redox homeostasis in host and microbiome

We also observed proteomics level evidence of altered lipid features. The increased frequency of protein modification by reactive lipid peroxidation products 4-hydroxynonenal (HNE) and 4-oxononenal (ONE) suggests oxidative stress in COVID-19 (Fig. 6E). Indeed, human superoxide dismutase (SOD1), the major antioxidant enzyme for superoxide removal and the first line of defense against oxidative stress, was significantly downregulated in COVID-19 (Fig. 2A). Meanwhile, NADH peroxidase, which reduces peroxides, of several bacterial species and genera belonging to order Clostridiales were also down regulated (Fig. 4D). These results indicate a redox homeostasis disruption for both host and gut bacteria.

GI tract is susceptible to SARS-COV-2 infection due to the high expression of ACE2 receptor. GI symptoms are frequently observed in patients with COVID-19. Gut's immune responses to SARS-CoV-2 necessitate greater attention because they can alter the commensal microbiome and the crosstalk between microbiota and extra-intestinal organ immunity. However, little is known about the importance of the enteric SARS-CoV-2 for the development of COVID-19-associated pathologies. Increasing evidence has shown that COVID-19 can promote cardiovascular disorders such as myocardial injury, acute coronary syndrome, and thromboembolism [42], neurologic symptoms such as myalgias, encephalopathy, and dizziness [43], and kidney manifestations such as proteinuria and dipstick hematuria [44]. Our study revealed an enrichment of gut bacteria related deleterious metabolites including phenylacetylglutamine (capable of causing cardiovascular diseases), neurotoxin salsolinol, and uremic toxin uric acid. In addition to metabolites, we observed a larger number of altered host and bacterial lipids (predominated by sphingolipids such as ceramide and hexosylceramide). Sphingolipids produced by gut bacteria can enter host metabolic pathways and impact host ceramide level [40]. Our study may provide an alternative microbiome-based molecular mechanism to explain how the gut ecosystem may play a role in the development of symptoms in COVID-19 and impact the host metabolome and lipidome.

The anti-viral response may impose an immunological off-target effects on gut microbiome in COVID-19 patients. Indeed, we observed disturbed mucosal immunological defense. The reduction of host Igs such as IGHV3-73 and IGHV3-64D may potentially contribute to the enrichment of Bacteroidetes phylum because of the reduced anti-bacteria Igs. The suppressed expression of proteins involved in neutrophil degranulation and migration can also impair the gut anti-bacteria defense system. Furthermore, there is an increased risk of colonic mucosal damage and therefore greater risk of viral and bacterial infection in COVID-19 because of the increased intestinal protease and glycosylation (indicating potential higher activity) and suppressed mucin glycosylation (important for mucin protection function). A major limitation of our study of is the limited sample size and further larger scale studies are needed. Nevertheless, our study has demonstrated widely disturbed gut molecular profiles which may play a role in the development of symptoms in COVID-19. Considering the gut ecosystem as a potential target could offer a valuable approach in managing the disease.

Using metaproteomics, metabolomics, glycoproteomics, and lipidomics, our study has demonstrated widely disturbed gut molecular profiles and microbial structure in COVID-19 characterized by disturbed immune, proteolysis and redox homeostasis. Our findings suggest that considering the gut ecosystem as a potential target could offer a valuable approach in managing the disease.

Subject details

This study involved 13 patients with COVID-19 and 21 healthy controls. The associated clinical data and metadata are provided in Table S1 and Fig. 1A. SARS-CoV-2 infection was confirmed by two consecutive real-time reverse transcription PCR (RT-PCR) tests. These patients were classified into three groups according to the severity of their symptoms: mild (mild clinical symptoms without pneumonia manifestations in CT imaging (7 patients)); moderate (respiratory symptoms, fever, and imaging features of COVID-19 pneumonia (5 patients); severe (respiratory distress (respiratory rate ≥ 30 breaths/min), oxygen saturation ≤ 93% and arterial oxygen tension (PaO₂)/fractional inspired oxygen (FiO₂) ratio ≤ 300 mm Hg (1 patient)). We collected 53 stool samples from COVID-19 patients with a range of one to nine longitudinal time-points that occurred 1–94 days post symptom onset (Fig. 1A, Table S1). Stool samples from 21 healthy subjects severed as controls. A total of 73 stool samples were subjected to multi-omics analysis.

Metaproteomics sample preparation

Fecal sample (~ 150 mg) was lysed by boiling for 5 min at 95°C in 800 µl of lysis buffer (6 M Guanidinium hydrochloride (GdmCl), 10 mM tris(2-carboxyethyl)phosphine (TCEP), 40 mM chloroacetamide, 100 mM Tris pH 8.5) in Eppendorf protein LoBind tubes. The lysate was then sonicated for 15 min using a waterbath sonicator. The crude protein extract was centrifuged at 16,000 g for 5 min with the clarified lysate subjected to ultrafiltration (cut off 30 kD, Millipore), diluted 1:10 with dilution buffer (10% (v/v) acetonitrile (ACN), 25 mM Tris pH 8.5) containing 1 µg sequencing grade trypsin (1/50, w/w), and digested overnight at 37°C. The digest was acidified to an end-concentration of 1% trifluoroacetic acid (TFA) and debris were removed after centrifugation at 16,000 g for 5 min. Finally, the peptides were desalted on StageTips assembled by Empore C18 disk, dried using a SpeedVac centrifuge at 45°C, and suspended in 2% ACN and 0.1% formic acid (FA).

Metaproteomics data acquisition

Peptides were trapped onto an Acclaim PepMap 100 C18 column (75 µm × 20 mm, 3 µm, 100 Å, Thermo Scientific) at a flow rate of 8 µL/min and separated using an Acclaim PepMap 100 C18 column (75 µm × 250 mm, 2 µm, 100 Å, Thermo Scientific) at 300 nL/min on a Thermo Scientific Dionex UltiMate 3000 RSLCnano LC system. Mobile phase solvents were 0.1% formic acid in water (A) and 0.1% formic acid in 80% acetonitrile (B). The separation gradient was as follows: 15% B rising to 30% B at 130 min, rising to 98% B at 137 min, and keeping at 98% B for 27 min. Finally, the separation column was equilibrated using 15% B for 6 min. The trap column was switched online with the separation column at 3 min and switched back to load position at 168 min. Data were acquired on an Orbitrap Fusion Lumos Tribrid mass spectrometer with a Nanospray Flex ion source in positive ionization mode with a spray voltage of + 2600 V using Xcalibur software (Thermo Scientific, San Jose, CA, USA). The ion transfer tube temperature was 300°C, the vaporized temperature was 325°C, the sheath gas flow was 40 units, the auxiliary gas flow was 15 arbitrary units, and the sweep gas was 1 unit. Full scan MS spectra was acquired in the 400 − 1,600 m/z range with an AGC target of 5 × 10⁴, a maximum injection time of 50 ms, and a resolution of 60K at m/z 200. MS/MS spectra were acquired using higher-energy collisional dissociation (HCD) with a normalized collision energy (NCE) of 30% and a resolution of 15K with an AGC target of 5× 10⁴ and a maximum ion injection time of 22 ms.

Glycopeptide enrichment

Glycopeptides were enriched using hydrophilic interaction liquid chromatography (HILIC) cartridges packed with the C18 plug followed by microcrystalline cellulose resins [45]. The resin was washed with 300 µL of 0.1% TFA and initialized using 300 µL of 0.1% TFA in 80% acetonitrile. After loading ~ 200 µg of peptides in 300 µL 80% acetonitrile/0.1% TFA), the resin was washed with 80% acetonitrile/0.1% TFA three times to remove non-specific peptides. Then glycopeptides were eluted by 300 µL of H₂O, followed by 200 µL of 80% acetonitrile. Peptides were dried using a SpeedVac centrifuge at 45°C, and suspended in 2% ACN and 0.1% formic acid (FA).

Glycoproteomics data acquisition

Peptides were trapped onto an Acclaim PepMap 100 C18 column (75 µm × 20 mm, 3 µm, 100 Å, Thermo Scientific) at a flow rate of 8 µL/min and separated using an Acclaim PepMap 100 C18 column (75 µm × 250 mm, 2 µm, 100 Å, Thermo Scientific) at 300 nL/min. Mobile phase solvents were 0.1% formic acid in water (A) and 0.1% formic acid in 50% acetonitrile and 40% isopropanol (B). The separation gradient was as follows: 5% B at 0–15 min, 20%-30% B at 90–100 min, and 98% B at 107 min and kept for 20 min. Data were acquired on an Orbitrap Fusion Lumos Tribrid mass spectrometer with a Nanospray Flex ion source in positive ionization mode with a spray voltage of + 2600 V using Xcalibur software (Thermo Scientific, San Jose, CA, USA). The ion transfer tube temperature was 300°C, the vaporized temperature was 325°C, the sheath gas flow was 40 units, the auxiliary gas flow was 15 arbitrary units, and the sweep gas was 1 unit. Full scan MS spectra was acquired in the 400 − 1,600 m/z range with a maximum injection time of 50 ms and a resolution of 60K at m/z 200. MS/MS spectra were acquired using HCD with stepped NCE at 20%, 30%, and 40% to generate fragment ions of both glycan and peptide of a glycopeptide in a single spectrum MS/MS spectra. The resolution of HCD was 15K with a maximum ion injection time of 22 ms.

Metabolomics sample preparation

Samples were prepared from ~ 150 mg stool samples by protein precipitation with the addition of 1 ml of 80% methanol. Samples were centrifuged (10 min, 16,000g, 4°C) and the supernatants were evaporated to dryness under nitrogen, reconstituted in 150 µl of 0.1% formic acid in 5% acetonitrile, and kept at -80°C until analysis.

Metabolomics data acquisition

Metabolic extracts were separated on a Thermo Scientific Dionex UltiMate 3000 Rapid Separation LC (RSLC) using an ACQUITY UPLC HSS T3 analytical column (2.1 × 150 mm, 1.8µm, 100 Å, Waters) protected by an ACQUITY UPLC HSS T3 VanGuard pre-column (2.1 × 5 mm, 1.8 µm, 100 Å, Waters). Mobile phase solvents for positive ionization mode were 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B); mobile phase solvents for negative ionization mode were 0.01% formic acid in water (A) and acetonitrile (B). The following gradient elution was used: 0–3 min, 95% A; 5–13 min, 80%-30% A; 15–18 min, 2% A; 18.1–22 min, 5% A. The flow rate was 0.3 mL/min, the injection volume was 2 µL and the column oven was set at 35°C. Data were acquired on an Orbitrap Fusion Lumos Tribrid mass spectrometer fitted with a HESI source in both positive and negative ionization modes with an independent run for each polarity and a spray voltage of + 3500 V and − 3500 V, respectively (Thermo Scientific, San Jose, CA, USA). The ion transfer tube temperature was 300°C, the vaporized temperature was 350°C, the sheath gas flow was 40 units, the auxiliary gas flow was 15 arbitrary units, and the sweep gas was 1 unit. Metabolite profiling was profiled in full scan mode using a mass range of m/z 100–1000 with a resolution of 120K at m/z 200, an AGC target of 5 × 10⁴, and a maximum injection time of 50 ms. For metabolite identification, data dependent MS/MS data were acquired on quality control samples (QC) containing equally volumes of all samples used in this study. In-depth MS/MS was performed using nine staggered gas-phase fractionations (sGPFs) to allow more homogeneous selection of precursor ions in low, medium, and high m/z ranges [46]. This was achieved in nine separated LC-MS runs: (run 1) 100–110, 200–210, 300–310, 400–410, 500–510, 600–610, 700–710, 800–810; (run 2) 110–120, 210–220, 310–320, 410–420, 510–520, 610–620, 710–720, 810–820; (run 3) 120–130, 220–230, 320–330, 420–430, 520–530, 620–630, 720–730, 820–830; (run 4) 130–140, 230–240, 330–340, 430–440, 530–540, 630–640, 730–740, 830–840; (run 5) 140–150, 240–250, 340–350, 440–450, 540–550, 640–650, 740–750, 840–850; (run 6) 150–160, 250–260, 350–360, 450–460, 550–560, 650–660, 750–760, 850–860; (run 6) 160–170, 260–270, 360–370, 460–470, 560–570, 660–670, 760–770, 860–870; (run 7) 170–180, 270–280, 370–380, 470–480, 570–580, 670–680, 770–780, 870–880; (run 8) 180–190, 280–290, 380–390, 480–490, 580–590, 680–690, 780–790, 880–890; (run 9) 190–200, 290–300, 390–400, 490–500, 590–600, 690–700, 790–800, 890–900. Each sGPF LC-MS run was performed twice. Quadrupole isolation window was 1.4 m/z and dynamic exclusion was enabled for 10s. The stepped NCE at 10%, 25%, and 40% was employed to obtain information-rich MS/MS spectra. The run order was the blank first (0.1% formic acid in 5% acetonitrile), pooled QC samples for DDA-MS/MS, and a pooled QC every 12 randomized clinical samples.

Lipidomics sample preparation

Extraction of lipids started with the addition of 1 mL of methanol to 150 mg of fecal samples and the tube was vigorously shaken with a vortex for 30 s [47]. Subsequently, 5 mL of methyl tertbutyl ether was added, vortexed for another 30 s, and shaken for 20 min at 200 rpm at room temperature. Next, phase separation was induced by adding 3 mL of ultrapure water with 2.5% trichloroacetic acid (w/v) and centrifugation for 5 min at 3000 rpm. Thereafter, 1 mL of the upper layer (consisting of methyl tert-butyl ether) was transferred and evaporated to dryness at 37°C under a gentle stream of nitrogen. The residue was sequentially resuspended in 250 µL of chloroform and 650 µL of methanol.

Lipidomics data acquisition

Lipid extracts (2 µL) were separated on a Thermo Scientific Dionex UltiMate 3000 Rapid Separation LC (RSLC) using an ACQUITY UPLC HSS T3 analytical column (2.1 × 150 mm, 1.8µm, 100 Å, Waters) protected by an ACQUITY UPLC HSS T3 VanGuard pre-column (2.1 × 5 mm, 1.8 µm, 100 Å, Waters). Mobile phase solvents A and B were ACN:H2O (6:4 v/v) and isopropanol: ACN (9:1 v/v), respectively, both contained 10 mM ammonium acetate and 0.1% acetic acid. The separation was performed at 55°C with a flow rate of 0.35 mL/min using the following gradient: 0–3.0 min, 30%-35% A; 5.0–14.0 min, 65%-98% A; 18.0-18.1 min, 98%-30% A; 18.1–22.0 min, 30% A. Data were acquired on an Orbitrap Fusion Lumos Tribrid mass spectrometer fitted with a HESI source in both positive and negative ionization modes with an independent run for each polarity and a spray voltage of + 3500 V and − 3500 V, respectively (Thermo Scientific, San Jose, CA, USA). The ion transfer tube temperature was 300°C, the vaporized temperature was 350°C, the sheath gas flow was 40 units, the auxiliary gas flow was 15 arbitrary units, and the sweep gas was 1 unit. Lipid profiling was profiled in DDA mode using a full MS scan range of m/z 150–2000 (resolution was 60K at m/z 200) with top ranked precursor ions subjected to DDA-MS/MS using a maximum injection time of 22 ms. The stepped normalized collision energy (NCE) at 25, 30, and 35 was employed to obtain information-rich MS/MS spectra with a resolution of 15K at m/z 200. Quadrupole isolation window was 1.6 m/z and dynamic exclusion was enabled for 10s. To promote lipid identification, in-depth DDA MS/MS of QC sample was performed using the following four sGPFs which was performed in four separated runs (Yan and Yan, 2015): (run 1) 150–250, 550–650, 950–1050, 1350–1450, 1750–1850; (run 2) 250–350, 650–750, 1050–1150, 1450–1550, 1850–1950; (run 3) 350–450, 750–850, 1150–1250, 1550–1650; (run 4) 450–550, 850–950, 1250–1350, 1650–1750.

Metaproteomics data analysis

Peptide identifications were performed using the search engine PEAKS DB combined with PEAKS de novo sequencing [48] (De Novo ALC(%) threshold was 15). False discovery rate (FDR) was set to 1% using the decoy fusion approach. Raw files were refined by precursor ion mass correction and resolving chimeric MS/MS spectra. The precursor mass tolerance was set to 15 ppm and the fragment mass tolerance to 0.03 Da. Enzyme specificity was set to trypsin and up to three missed cleavage sites were allowed. The maximum number of variable posttranslational modifications per peptide was three, including acetylation of protein N-terminus, carbamidomethylation of Cys, oxidation of Met, deamidation of Asn and Gln as well as Pyro-glu from Gln. PEAKS PTM search tool [49] was used to search for peptides with unspecified modifications (313 built-in post-translational modifications), and the SPIDER [50] search tool was used for exploring novel peptides that are homologous to peptides in the protein database.

Database search was performed using a comprehensive meta-database containing human, microbial, and dietary organism sequences [51]. The gut microbial protein database was generated by combining the following parts: (1) the integrated gene catalog of 1,267 human fecal metagenomes [52]; (2) the 1,520 reference genomes of > 6,000 cultivated human fecal bacteria isolates [53]; (3) the genomes of 215 human fecal bacteria isolates [54]; (4) all Archaea, Bacteria, and Fungi sequences in NCBI RefSeq (Release 90) and UniProtKB (Release 2017_06). The microbial database was appended by the SARS-COV-2 protein sequences [55], a UniProt human reference proteome (downloaded on 2017_06), and a food database of common dietary organisms. A total number of 130,975,891 non-redundant sequences were obtained after dereplicating at 100% amino acid identity using USEARCH v11.0.667 (–fastx_uniques) [56]. Proteins identified by at least one unique peptides (1% false discovery rate (FDR) using the decoy fusion approach) was considered for further analysis. Label-free quantification of protein groups was performed based on the number of peptide spectrum matches (PSM).

Taxonomy and functional analysis of gut microbiota

Taxonomy and functional analysis of peptides was performed with UniPept (version 4.3.7) [57] based on the lowest common ancestor (LCA) algorithm using the following parameters: Equate I and L, Advanced missing cleavage handling. Peptide functional annotations were performed using Gene Ontology (GO) terms and Enzyme Commission (EC) numbers. The relative abundance of microbial taxonomic and functional groups were determined using the normalized number of corresponding peptides. Functions of the unannotated microbial proteins were predicted using protein-protein BLAST (BlastP, https://blast.ncbi.nlm.nih.gov/Blast.cgi) against the non-redundant protein sequences (nr) with an E-value threshold of 1e-10.

Glycoproteomics data analysis

High-confidence identification of intact N-glycopeptides was performed by pGlyco 2.0 [20]. Sequences of proteins identified in the above metaproteomics analysis as well as the SARS-COV-2 protein sequences were used in glycopeptide identification. The precursor mass tolerance was set to 10 ppm and the fragment mass tolerance to 20 ppm. For N-glycopeptide analysis, a FDR of 5% at the glycan level and a FDR of 1% at the peptide level was used. Match between run and intensity based label-free quantification was not supported by pGlyco. To calculate the frequency of N-glycosylation at a specific site, any sample containing a glycopeptide bearing an N-glycan at that site was considered positive regard less of the variations of glycopeptide sequences and glycans. To calculate the frequency of a specific N-glycan at a specific site, any sample containing a glycopeptide bearing that glycan at that site was considered positive regard less of glycopeptide sequence variations. The frequency was defined as the number of positive samples relative to the total number of samples in each group (18 for the control and 49 for the COVID-19 groups, respectively).

Global identification of both O-glycopeptides and N-glycopeptides were performed using PEAKS. Acetylation (protein N-term), carbamidomethylation, deamidation (NQ), pyro-glu from Q, oxidation (M), HexNAcylation (ST), Hex1HexNAc1, Hex1HexNAc2, Hex2HexNAc2, Hex3HexNAc2, Hex1HexNAc(3), Hex(2)HexNAc(3), Hex(3)HexNAc(3), Hex(4)HexNAc(2), Hex(6)HexNAc(2), Hex(7)HexNAc(2), Hex(8)HexNAc(2), Hex(9)HexNAc(2), and dHexHex(3)HexNAc(2) were set as variable posttranslational modifications. The precursor mass tolerance was set to 10 ppm and the fragment mass tolerance to 0.02 Da. The PEAKS output file gave the monosaccharide composition of the attached glycan. Label-free quantification was performed using match between runs with a mass error tolerance of 20 ppm and a retention time shift tolerance of 1 min.

Metabolomics data analysis

Metabolomics features were extracted, aligned, identified and quantified using Compound Discoverer (v3.1, Thermo Fisher Scientific). The analysis employed the following major steps and parameters: retention time alignment (alignment model = adaptive curve, mass tolerance = 5 ppm, maximum shift = 2 min), unknown compound detection (mass tolerance = 5 ppm, intensity threshold = 30%, S/N threshold = 3, minimum peak intensity = 1 × 10⁶, adducts ions = [M + H]⁺¹, [M + H-H2O]⁺¹, [M + H-NH3]⁺¹, [M + K]⁺¹, [M + Na]⁺¹, [M + NH4]⁺¹, [2M + H]⁺¹, [2M + K]⁺¹, [2M + Na]⁺¹, [2M + NH4]⁺¹, [M + 2H]⁺², [M-H]⁻¹, [M-2H]⁻², [M-H + HAc]⁻¹, [M-H-H2O]⁻¹, [M-H + FA]⁻¹, [M + Cl]⁻¹, [2M-H]⁻¹, [2M-H + HAc]⁻¹), compound grouping (mass tolerance = 5 ppm, RT tolerance = 0.2 min), prediction of elemental compositions (mass tolerance = 5 ppm, maximum element counts = 90 × C, 190 × H, 10 × N, 15 × O, 5 × S and 3 × P), filling gaps across all samples (mass tolerance = 5 ppm,, S/N threshold = 1.5), chemical background subtraction (using blank samples), identifying compounds by searching ChemSpider (by formula or mass, https://www.chemspider.com/), mzVault and mzCloud (by MS and MS/MS data, precursor mass tolerance = 10 ppm, fragment mass tolerance = 10 ppm, match factor threshold = 60, https://www.mzcloud.org), and QC-based batch normalization (regression model = Cubic Spline). The mzCloud and mzVault match were performed base on Similarity Forward method and HighChem-HighRes search algorithm, respectively. Extracted ion chromatogram (EIC) and MS/MS spectra of all metabolites of interests were manually inspected.

Lipidomics data analysis

Raw data files were processed using the LipidSearch software (version 4.1) (Thermo Fisher Scientific) to identify and quantify lipid molecular species. Peak detection was performed as follows: Recalc Isotope, on; RT interval (min), 0.01. Lipid identification was as follows: Search type, Product; Exp type, LC-MS; Precursor tolerance, 10 ppm; Product tolerance, 10 ppm; Intensity threshold, 1.0%; Target class, ALL lipid classes; Ion adducts (positive ion mode) of + H, +NH₄, +Na, +H-H₂O, and + 2 H; Ion adducts (negative Ion mode) of − H, +HCOO, +CH₃COO, +Cl, and − 2H; Top rank filter, On; Main node filter, Main isomer peak; m-Score threshold, 5.0; FA priority, On. ID Quality Filter, Check A, B, C, D (A: lipid class and FA are identified, B: lipid class and some FA were identified, C: Lipid class or FA were identified, D: Lipid identified by other fragment ions (H₂O loss, and other non-specific neutral losses). Quantitation was performed using a m/z tolerance of −/+ 5.0 ppm and a RT range of − 0.5/+ 0.5 min. Peak alignment was performed using the following parameters: Alignment Method, Max; RT Tolerance, 0.25 min; Calculate unassigned peak area, On; Filter type, New filter; Top rank filter, On; Main node filter, Main isomer peak; m-Score threshold, 5.0; ID quality filter, A, B and C.

Statistical analysis

The raw quantification data matrix of different omics was imported to MetaboAnalyst [58] for further processing and analysis. Data filtering was performed using interquantile range (IQR) to remove baseline noises. Missing values were imputed using KNN. Quantile normalization and pareto scaling were employed. Unsupervised multivariate data analysis was performed using principal component analysis (PCA) and hierarchical cluster analysis (HCA). Significantly differentiated omics features between COVID-19 and control groups (should present in at least 50% of samples) were detected using Wilcoxon’s rank sum test (FDR adjusted p value (q) < 0.05). Microbial taxonomic and functional groups were normalized by total abundance. Statistical significance of microbial taxonomic groups was calculated using Mann-Whitney U test (p value < 0.05). Protein-microbiome and metabolite-microbiome correlations were determined using Pearson correlation (q < 0.05) using R. The association between differentiating omics features with categorical confounding variables (gender and medicine) were determined using Wilcoxon’s rank sum test. The association between differentiating omics features with continuous confounding variables (age and Body Mass Index (BMI)) were determined using Pearson correlation.

Ethics approval and consent to participate

This study was approved by the Ethics Committee of The Fifth Affiliated Hospital, Sun Yat-sen University (K161-1).

Consent for publication

Not applicable.

Availability of data and materials

The multi-omics data generated for this manuscript have been deposited in ProteomeXchange Consortium (https://www.iprox.org/) under the following identifier: Metaproteomics (IPX0002453001), Metabolomics (IPX0002453002), Lipidomics (IPX0002453003), and Glycoproteomics (IPX0002453004).

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China (31900070), Guangdong Basic and Applied Basic Research Foundation (2019A1515011771), Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai, SML2020SP003), the National Key Research and Development Program of China (2020YFC082400), Science and Technology Development Fund of Macau SAR (FDCT0017/2020/A), the Task-Force Project on the Prevention and Control of Novel Coronavirus of Guangdong Province (20201113), the Three Major Constructions of Sun Yat-sen University (the Task-Force Project on the Prevention and Control of Novel Coronavirus of Sun Yat-sen University), the Emergency Task-Force of SARS-CoV-2 research of Guangzhou Regenerative Medicine and Health Guangdong Laboratory, the Emergency Task-Force Project on the Prevention and Control of Novel Coronavirus of Zhuhai 2020.

Authors' contributions

Conceptualization, Z.Y.; Data analysis, Z.Y.; Methodology, Z.Y.; Sample coordination and preparation, F.H., Z.Y., and F.X.; Sample collection, Z.F., S.H., Z.G., H.S., Z.Z., and H.Z.; Clinical data collection, Z.Y., H.Z., L.L., and J.L.; Resources, F.X., H.H., D.L., and P.P.; Clinical laboratory tests, G.J., and K.L.; Writing - Original Draft, Z.Y.; Writing - Review & Editing, Z.Y., H.S., X.L., X.F., and R.Y.; Work supervised by Z.Y., H.S., X.L., X.F., and R.Y.

Lin, L., Jiang, X., Zhang, Z., Huang, S., Zhang, Z., Fang, Z., Gu, Z., Gao, L., Shi, H., Mai, L., et al. (2020). Gastrointestinal symptoms of 95 cases with SARS-CoV-2 infection. Gut 69, 997-1001.
Jin, X., Lian, J.S., Hu, J.H., Gao, J., Zheng, L., Zhang, Y.M., Hao, S.R., Jia, H.Y., Cai, H., Zhang, X.L., et al. (2020). Epidemiological, clinical and virological characteristics of 74 cases of coronavirus-infected disease 2019 (COVID-19) with gastrointestinal symptoms. Gut 69, 1002-1009.
Song, Y., Liu, P., Shi, X.L., Chu, Y.L., Zhang, J., Xia, J., Gao, X.Z., Qu, T., and Wang, M.Y. (2020). SARS-CoV-2 induced diarrhoea as onset symptom in patient with COVID-19. Gut 69, 1143-1144.
Xiao, F., Tang, M., Zheng, X., Liu, Y., Li, X., and Shan, H. (2020). Evidence for Gastrointestinal Infection of SARS-CoV-2. Gastroenterology 158, 1831-1833.
Zhang, Y., Chen, C., Zhu, S., Shu, C., Wang, D., Song, J., Song, Y., Zhen, W., Feng, Z., Wu, G., et al. (2020). Isolation of 2019-nCoV from a stool specimen of a laboratory-confirmed case of the coronavirus disease 2019 (COVID-19). China CDC Weekly 2, 123-124.
Hoffmann, M., Kleine-Weber, H., Schroeder, S., Krüger, N., Herrler, T., Erichsen, S., Schiergens, T.S., Herrler, G., Wu, N.H., Nitsche, A., et al. (2020). SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor. Cell 181, 271-280.
Zuo, T., Zhan, H., Zhang, F., Liu, Q., Tso, E., Lui, G., Chen, N., Li, A., Lu, W., Chan, F., et al. (2020). Alterations in fecal fungal microbiome of patients with COVID-19 during time of hospitalization until discharge. Gastroenterology S0016-5085, 34852-34856.
Zuo, T., Zhang, F., Lui, G., Yeoh, Y.K., Li, A., Zhan, H., Wan, Y., Chung, A., Cheung, C.P., Chen, N., et al. (2020). Alterations in gut microbiota of patients with COVID-19 during time of hospitalization. Gastroenterology S0016-5085, 34701-34706.
Arunachalam, P.S., Wimmers, F., Mok, C., Perera, R., Scott, M., Hagan, T., Sigal, N., Feng, Y., Bristow, L., Tak-Yin Tsang, O., et al. (2020). Systems biological assessment of immunity to mild versus severe COVID-19 infection in humans. Science 369, 1210-1220.
0 Stanley, D., Mason, L.J., Mackin, K.E., Srikhanta, Y.N., Lyras, D., Prakash, M.D., Nurgali, K., Venegas, A., Hill, M.D., Moore, R.J., et al. (2016). Translocation and dissemination of commensal bacteria in post-stroke infection. Nat. Med. 22, 1277-1284.
Grau, K.R., Zhu, S., Peterson, S.T., Helm, E.W., Philip, D., Phillips, M., Hernandez, A., Turula, H., Frasse, P., Graziano, V.R., et al. (2020). The intestinal regionalization of acute norovirus infection is regulated by the microbiota via bile acid-mediated priming of type III interferon. Nat. Microbiol. 5, 84-92.
Hagan, T., Cortese, M., Rouphael, N., Boudreau, C., Linde, C., Maddur, M.S., Das, J., Wang, H., Guthmiller, J., Zheng, N.Y., et al. (2019). Antibiotics-driven gut microbiome perturbation alters immunity to vaccines in humans. Cell 178, 1313-1328.
Toy, L.S., Yio, X.Y., Lin, A., Honig, S., and Mayer, L. (1997). Defective expression of gp180, a novel CD8 ligand on intestinal epithelial cells, in inflammatory bowel disease. J. Clin. Invest. 100, 2062-2071.
Roda, G., Dahan, S., Mezzanotte, L., Caponi, A., Roth-Walter, F., Pinn, D., and Mayer, L. Defect in CEACAM family member expression in Crohn's disease IECs is regulated by the transcription factor SOX9. Inflamm Bowel Dis. 2009 15(12):1775-83.
Roda, G., Jianyu, X., Park, M.S., DeMarte, L., Hovhannisyan, Z., Couri, R., Stanners, C.P., Yeretssian, G., and Mayer, L. (2014). Characterizing CEACAM5 interaction with CD8α and CD1d in intestinal homeostasis. Mucosal Immunol. 7, 615-624.
Burgener, S.S., Leborgne, N., Snipas, S.J., Salvesen, G.S., Bird, P.I., and Benarafa, C. (2019). Cathepsin G Inhibition by Serpinb1 and Serpinb6 Prevents Programmed Necrosis in Neutrophils and Monocytes and Reduces GSDMD-Driven Inflammation. Cell Rep. 27, 3646-3656.
Toes, R.E., Toebes, M., Schumacher, T.N., Bladergroen, B.A., Ossendorp, F., Kummer, J.A., Melief, C.J., and Offringa, R. (2001). Expression of the serpin serine protease inhibitor 6 protects dendritic cells from cytotoxic T lymphocyte-induced apoptosis: differential modulation by T helper type 1 and type 2 cells. J. Exp. Med. 194, 657-667.
Bates, J.M., Akerlund, J., Mittge, E., and Guillemin, K. (2007). Intestinal alkaline phosphatase detoxifies lipopolysaccharide and prevents inflammation in zebrafish in response to the gut microbiota. Cell Host Microbe. 2, 371-382.
Zauner, G., Selman, M.H., Bondt, A., Rombouts, Y., Blank, D., Deelder, A.M., and Wuhrer, M. (2013). Glycoproteomic analysis of antibodies. Mol. Cell. Proteomics 12, 856-865.
Liu, M.Q., Zeng, W.F., Fang, P., Cao, W.Q., Liu, C., Yan, G.Q., Zhang, Y., Peng, C., Wu, J.Q., Zhang, X.J., et al. (2017). pGlyco 2.0 enables precision N-glycoproteomics with comprehensive quality control and one-step mass spectrometry for intact glycopeptide identification. Nat. Commun. 8, 438.
Kleiner, M., Thorson, E., Sharp, C.E., Dong, X., Liu, D., Li, C., and Strous, M. (2017). Assessing species biomass contributions in microbial communities via metaproteomics. Commun. 8, 1558.
Vital, M., Karch, A., and Pieper, D.H. (2017). Colonic butyrate-producing communities in humans: an overview using omics data. mSystems 2, e00130-17.
Udayappan, S., Manneras-Holm, L., Chaplin-Scott, A., Belzer, C., Herrema, H., Dallinga-Thie, G.M., Duncan, S.H., Stroes, E., Groen, A.K., Flint, H.J., et al. (2016). Oral treatment with Eubacterium hallii improves insulin sensitivity in db/db mice. NPJ biofilms microbiomes 2, 16009.
Arnoriaga-Rodríguez, M., Mayneris-Perxachs, J., Burokas, A., Contreras-Rodríguez, O., Blasco, G., Coll, C., Biarnés, C., Miranda-Olivos, R., Latorre, J., Moreno-Navarrete, J.M., et al. (2020). Obesity impairs short-term and working memory through gut microbial metabolism of aromatic amino acids. Cell Metab. 32, 548-560.
Geva-Zatorsky, N., Sefik, E., Kua, L., Pasman, L., Tan, T.G., Ortiz-Lopez, A., Yanortsang, T.B., Yang, L., Jupp, R., Mathis, D., et al. (2017). Mining the human gut microbiota for immunomodulatory organisms. Cell168, 928-943.
Sheikh, A., Tumala, B., Vickers, T.J., Alvarado, D., Ciorba, M.A., Bhuiyan, T.R., Qadri, F., Singer, B.B., and Fleckenstein, J.M. (2020). CEACAMs serve as toxin-stimulated receptors for enterotoxigenic Escherichia coli. Proc Natl Acad Sci U S A. 17, 29055-29062.
Nemet, I., Saha, P.P., Gupta, N., Zhu, W., Romano, K.A., Skye, S.M., Cajka, T., Mohan, M.L., Li, L., Wu, Y., et al. (2020). A Cardiovascular Disease-Linked Gut Microbial Metabolite Acts via Adrenergic Receptors. Cell 180, 862-877.
Villageliú, D.N., Borts, D.J., and Lyte, M. (2018). Production of the neurotoxin salsolinol by a gut-associated bacterium and its modulation by alcohol. Front. Microbiol. 9, 3092.
Kurnik-Łucka, M., Panula, P., Bugajski, A., and Gil, K. (2018). Salsolinol: an unintelligible and double-faced molecule-lessons learned from in vivo and in vitro experiments. Neurotoxicity Res. 33, 485-514.
Sorensen, L.B. (1965). Role of the intestinal tract in the elimination of uric acid. Arthritis Rheum. 8, 694-706.
Hosomi, A., Nakanishi, T., Fujita, T., and Tamai, I. (2012). Extra-renal elimination of uric acid via intestinal efflux transporter BCRP/ABCG2. PloS one. 7, e30456.
Zhu, L., Baker, S.S., Gill, C., Liu, W., Alkhouri, R., Baker, R.D., and Gill, S.R. (2013). Characterization of gut microbiomes in nonalcoholic steatohepatitis (NASH) patients: a connection between endogenous alcohol and NASH. Hepatology 57, 601-609.
Michail, S., Lin, M., Frey, M.R., Fanter, R., Paliy, O., Hilbush, B., and Reo, N.V. (2015). Altered gut microbial energy and metabolism in children with non-alcoholic fatty liver disease. FEMS Microbiol. Ecol. 91, 1-9.
Yuan, J., Chen, C., Cui, J., Lu, J., Yan, C., Wei, X., Zhao, X., Li, N., Li, S., Xue, G., et al. (2019). Fatty liver disease caused by high-alcohol-producing Klebsiella pneumoniae. Cell Metab. 30, 675-688.
Baranowski, S., Serr, A., Thierauf, A., Weinmann, W., Grosse Perdekamp, M., Wurst, F.M., and Halter C.C. (2008). In vitro study of bacterial degradation of ethyl glucuronide and ethyl sulphate. Int. J. Legal Med. 122, 389-393.
Shapiro, H., Kolodziejczyk, A.A., Halstuch, D. and Elinav, E. (2018). Bile acids in glucose metabolism in health and disease. J. Exp. Med. 215, 383-396.
Quinn, R.A., Melnik, A.V., Vrbanac, A., Fu, T., Patras, K.A., Christy, M.P., Bodai, Z., Belda-Ferre, P., Tripathi, A., Chung, L.K., et al. (2020). Global chemical effects of the microbiome include new bile-acid conjugations. Nature 579, 123-129.
An, D., Oh, S.F., Olszak, T., Neves, J.F., Avci, F.Y., Erturk-Hasdemir, D., Lu, X., Zeissig, S., Blumberg, R.S., and Kasper, D.L. (2014). Sphingolipids from a symbiotic microbe regulate homeostasis of host intestinal natural killer T cells. Cell 156, 123-133.
Brown, E.M., Ke, X., Hitchcock, D., Jeanfavre, S., Avila-Pacheco, J., Nakata, T., Arthur, T.D., Fornelos, N., Heim, C., Franzosa, E.A., et al. (2019). Bacteroides-Derived Sphingolipids Are Critical for Maintaining Intestinal Homeostasis and Symbiosis. Cell Host Microbe 25, 668-680.
Johnson, E.L., Heaver, S.L., Waters, J.L., Kim, B.I., Bretin, A., Goodman, A.L., Gewirtz, A.T., Worgall, T.S., and Ley, R.E. (2020). Sphingolipids produced by gut bacteria enter host metabolic pathways impacting ceramide levels. Nat. Commun. 11, 2471.
Stoffel, W., Dittmar, K., and Wilmes, R. (1975). Sphingolipid metabolism in Bacteroideaceae. Hoppe-Seyler’s Z. Physiol. Chem. 356, 715-726.
Madjid, M., Safavi-Naeini, P., Solomon, S.D. and Vardeny, O. (2020). Potential effects of coronaviruses on the cardiovascular system: a review. JAMA Cardiol. 5, 831-840.
Liotta, E.M., Batra, A., Clark, J.R., Shlobin, N.A., Hoffman, S.C., Orban, Z.S., and Koralnik, I.J. (2020). Frequent neurologic manifestations and encephalopathy‐associated morbidity in Covid‐19. Ann. Clin. Transl. Neurol. 7, 2221-2230.
Hong, D., Long, L., Wang, A.Y., Lei, Y., Tang, Y., Zhao, J.W., Song, X., He, Y., Wen, E., Zheng, L., et al. (2020). Kidney manifestations of mild, moderate and severe coronavirus disease 2019: a retrospective cohort study. Clin. kidney J. 13, 340-346.
Ohta, Y., Kameda, K., Matsumoto, M., and Kawasaki, N. (2017). Rapid glycopeptide enrichment using cellulose hydrophilic interaction/reversed-phase stageTips. Mass Spectrom (Tokyo). 6, A0061.
Yan, Z., and Yan, R. (2015). Improved data-dependent acquisition for untargeted metabolomics using gas-phase fractionation with staggered mass range. Anal Chem. 87, 2861-2868.
Van Meulebroek, L., De Paepe, E., Vercruysse, V., Pomian, B., Bos, S., Lapauw, B., and Vanhaecke, L. (2017). Holistic lipidomics of the human gut phenotype using validated ultra-high-performance liquid chromatography coupled to hybrid orbitrap mass spectrometry. Anal. Chem. 89, 12502-12510.
Zhang, J., Xin, L., Shan, B., Chen, W., Xie, M., Yuen, D., Zhang, W., Zhang, Z., Lajoie, G.A., and Ma, B. (2012). PEAKS DB: de novo sequencing assisted database search for sensitive and accurate peptide identification. Mol. Cell. Proteomics 11, M111.010587.
Han, X., He, L., Xin, L., Shan, B., and Ma, B. (2011). PeaksPTM: mass spectrometry-based identification of peptides with unspecified modifications. J. Proteome Res. 10, 2930-2936.
Han, Y., Ma, B., and Zhang, K. (2005). SPIDER: software for protein identification from sequence tags with de novo sequencing error. J. Bioinform. Comput. Biol. 3, 697-716.
Yan, Z., He, F., Xiao, F., He, H., Li, D., Cong, L., Lin, L., Zhu, H., Wu, Y., Yan, R., et al. (2020). A semi-tryptic peptide centric metaproteomic mining approach and its potential utility in capturing signatures of gut microbial proteolysis. Microbiome 9.
Li, J., Jia, H., Cai, X., Zhong, H., Feng, Q., Sunagawa, S., Arumugam, M., Kultima, J. R., Prifti, E., Nielsen, T., et al. (2014). An integrated catalog of reference genes in the human gut microbiome. Nat. Biotechnol. 32, 834-841.
Zou, Y., Xue, W., Luo, G., Deng, Z., Qin, P., Guo, R., Sun, H., Xia, Y., Liang, S., Dai, Y., et al. (2019). 1,520 reference genomes from cultivated human gut bacteria enable functional microbiome analyses. Nat. Biotechnol. 37, 179-185.
Browne, H.P., Forster, S.C., Anonye, B.O., Kumar, N., Neville, B.A., Stares, M.D., Goulding, D., and Lawley, T.D. (2016). Culturing of 'unculturable' human microbiota reveals novel taxa and extensive sporulation. Nature 533, 543-546.
Gordon, D.E., Jang, G.M., Bouhaddou, M., Xu, J., Obernier, K., White, K.M., O'Meara, M.J., Rezelj, V.V., Guo, J.Z., Swaney, D.L., et al. (2020). A SARS-CoV-2 protein interaction map reveals targets for drug repurposing. Nature 583, 459-468.
Edgar R.C. (2010). Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460-1.
Gurdeep Singh, R., Tanca, A., Palomba, A., Van der Jeugt, F., Verschaffelt, P., Uzzau, S., Martens, L., Dawyndt, P., and Mesuere, B. (2019). Unipept 4.0: functional analysis of metaproteome data. J. Proteome Res. 18, 606-615.
Chong, J., Wishart, D.S., and Xia, J. (2019). Using MetaboAnalyst 4.0 for comprehensive and integrative metabolomics data analysis. Curr Protoc Bioinformatics 68, e86.

AdditionalfileTableS1S7FigureS1.rar
Supplementary file legends Figure S1. Deleterious metabolite identification by comparison between experimental and the reference MS/MS spectra, related to Figure 5. A-D. The experimental (upper penal) and reference (lower penal) MS/MS spectra show high similarities. Structures of key fragments are annotated.
Table S1. Individual pathological and clinical details of the COVID-19 patients and controls, Related to Figure 1.
Table S2. Expression matrices of global metaproteome, Related to Figure 2. 2.1. All identified protein groups. 2.2. All identified peptides. 2.3. Altered human protein groups. 2.4. Altered microbial protein groups.
Table S3. Expression matrices of HILIC enriched metaproteome, Related to Figure 3. 3.1. All identified protein groups. 3.2. All identified peptides. 3.3. Altered human protein groups. 3.4. Altered Microbial protein groups.
Table S4. N-glycoproteomics and O-glycoproteomics, Related to Figure 3. 4.1. PSM of all identified human N-glycopeptides using pGlyco. 4.2. Frequency of site-specific protein N-glycosylation. 4.3. Frequency of site-specific and glycan-specific protein N-glycosylation. 4.4. LFQ abundance of all identified glycopeptides using PEAKS. 4.5. Altered human O-glycopeptides in COVID-19 using PEAKS.
Table S5. Relative abundance of taxonomic and functional groups of gut microbiome, Related to Figure 4. S5.1. Relative abundance of taxonomic groups of gut microbiome. S5.2. Relative abundance of functional groups of gut microbiome. S5.3. Associations between altered host proteins and gut microbiome.
Table S6. Altered metabolome, Related to Figure 5. S6.1. Metabolites differentiating COVID-19 patients from healthy controls. S6.2. Correlations between altered gut microbiome and metabolites.
Table S7. Expression matrices and differential analysis of lipidomics dataset, Related to Figure 6.

Download PDF

Version 1

posted

You are reading this latest preprint version

Molecular Phenomics of Gut Ecosystem in COVID-19 Patients

Status:

Version 1

Abstract

Figures

Introduction

Results

Disturbed host proteome

Glycosylation insight into mucosal immunological pathogenesis

Reduction of beneficial gut bacteria and potential host-bacteria interactions

Functional alteration of gut microbiome

Enrichment of bacterial related deleterious metabolites

Alerted microbial lipidome profiles

Increased lipid peroxidation and disturbed redox homeostasis in host and microbiome

Disscussion

Conclusions

Materials And Methods

Subject details

Metaproteomics sample preparation

Metaproteomics data acquisition

Glycopeptide enrichment

Glycoproteomics data acquisition

Metabolomics sample preparation

Metabolomics data acquisition

Lipidomics sample preparation

Lipidomics data acquisition

Metaproteomics data analysis

Taxonomy and functional analysis of gut microbiota

Glycoproteomics data analysis

Metabolomics data analysis

Lipidomics data analysis

Statistical analysis

Declarations

References

Supplementary Files

Status:

Version 1