The study was approved by the Ethics Committee of the promoting center (Comitato Etico Indipendente di Area Vasta Emilia Centro, n. 283/2020/Oss/AOUBo).

We prospectively enrolled consecutive adult (≥18 years) patients diagnosed of SARS-CoV-2 virus infection, hospitalized from April 21 to May 9, 2020, in two teaching hospitals (St. Orsola and Bellaria Hospital) and one tertiary non-teaching hospital (Maggiore Hospital) from the metropolitan area of Bologna (Italy). Subjects were excluded if they had not a laboratory confirmed diagnosis of SARS-CoV-2 infection, and/or clinical data were unavailable. Underlying conditions were recorded according to Charlson comorbidity index (Charlson et al., 1987). As for SARS-CoV-2 infection, we collected: i) date of symptoms onset; ii) date and symptoms at hospitalization; iii) vital signs, laboratory tests and imaging findings at hospitalization; iv) clinical severity at hospitalization classified according to sequential organ failure assessment; and v) all administered antiviral and immunomodulatory treatments. Admission to ICU, management of respiratory failure, renal replacement therapy and inotropic support, as well as don’t resuscitate order established by attending physician were recorded. Occurrence of bacterial superinfection was assessed according to Centers for Disease Control and Prevention (CDC) criteria (Horan et al., 2008). A fecal sample was collected at the time of infectious disease consultation in a sterile plastic container and kept at -80°C until further processing.

Microbiological diagnosis of SARS-CoV-2 infection was performed by detection of SARS-CoV-2 RNA in respiratory samples (oropharyngeal-nasopharyngeal swab, bronchoalveolar lavage or broncoaspirate), as previously described (Bartoletti et al., 2020). Briefly, total DNA/RNA was extracted from samples by QiaSymphony (QIAGEN) and detection of SARS-CoV-2 was performed by real-time RT-PCR targeting regions in the N gene following the US CDC protocol.

All fecal samples were processed in a biosafety level 3 laboratory [Centro di riferimento regionale per le emergenze microbiologiche (CRREM), St. Orsola Hospital, Bologna, Italy]. Microbial DNA was extracted using the repeated bead-beating plus column method (Yu and Morrison, 2004), with slight modifications as previously described (Turroni et al., 2017). For library preparation, the V3-V4 hypervariable region of the 16S rRNA gene was amplified by using the 341F and 785R primers (D’Amico et al., 2019). The final libraries, indexed and purified, were sequenced on an Illumina MiSeq platform, with a 2 × 250 bp paired-end protocol according to the manufacturer’s instructions (Illumina). Sequencing reads were deposited in the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA; BioProject ID PRJNA700830). Further details are available in the Supplementary Methods .

Raw sequences were processed using a combined pipeline of PANDASeq (Masella et al., 2012) and QIIME 2 (Bolyen et al., 2019). Length and quality-filtered reads were binned into amplicon sequence variants (ASVs) using DADA2 (Callahan et al., 2016). Taxonomic assignment was performed using VSEARCH (Rognes et al., 2016) against the Greengenes database (May 2013 release). Chimeras were discarded during analysis. Publicly available 16S rRNA gene sequences of 69 healthy subjects matched by age [median, IQR (days): 59, 48-95], sex (39 females, 30 males) and geography (across Italy) were downloaded from databases and processed as above [33 subjects: NCBI SRA, Bioproject ID SRP042234 (De Filippis et al., 2016); 36 subjects: MG-RAST ID mgp17761 (Biagi et al., 2016)]. GM sequences of 16 patients admitted to ICU at St. Orsola Hospital after undergoing liver transplantation in October 2019-February 2020, were also used for comparative purposes (unpublished data). We selected non-COVID-19 patients who were best matched for known microbiota-associated confounding factors (i.e., age [median, IQR (years): 57, 47.8-63] and sex ratio [50%]) as well as for other therapy and hospitalization-related confounders (i.e., exposure to antibiotics in the 2 weeks prior to fecal sampling [69%], time interval between hospitalization and sampling [median, IQR (days): 12, 8-15.5], and development of bacterial infection after sampling [12.5%]). All fecal samples had been collected by the authors and processed in the same laboratory, then subjected to the same wet and in silico analysis steps. Alpha diversity was calculated using the number of observed ASVs and the inverse Simpson index, while beta diversity was estimated by computing Bray-Curtis distances between the genus-level profiles, which were then used as input for Principal Coordinates Analysis (PCoA).

All statistical analysis was performed in R (https://www.r-project.org/). PCoA plots were generated using the “vegan” (http://www.cran.r-project.org/package-vegan/) and “Made4” (Culhane et al., 2005) packages, and data separation was tested by a permutation test with pseudo-F ratio (function “Adonis” in “vegan”). For each PCoA plot, ellipses including 95% confidence area based on the standard error of the weighted average of sample coordinates were overlaid. The bacterial genera most contributing to the ordination space were identified using the function “envfit” of “vegan”. Linear discriminant analysis (LDA) effect size (LEfSe) algorithm with LDA score threshold of 3 (on a log10 scale) was applied on genus-level tables to identify discriminating taxa (Segata et al., 2011). Group differences in alpha diversity and taxon relative abundance were assessed by Wilcoxon test or Kruskal-Wallis test followed by post-hoc comparisons. P values were corrected for multiple comparisons using the Benjamini–Hochberg method. A false discovery rate (FDR) ≤ 0.05 was considered statistically significant. Sequences assigned to Enterococcus were aligned to the NCBI 16S rRNA database (September 2019 release) using BLASTn (Altschul et al., 1990); only hits with ≥80% identity were considered.

As for clinical data, categorical variables are presented as absolute numbers and relative frequencies; continuous variables are presented as mean and standard deviation if normally distributed or as median and IQR if non-normally distributed. The incidence rate of BSI due to Enterococcus spp. (i.e., E-BSI) was calculated as the number of E-BSIs divided by the total number of 10,000 person-days at risk in ICU at St. Orsola Hospital in Bologna. The incidence rate was calculated for the first 4 months of 2020 (during the COVID-19 pandemic) and for the last three years over the same period (January-April). These values were compared by Poisson regression.

Of the 76 enrolled patients, 5 were excluded due to unavailable clinical data and 2 because of unconfirmed SARS-CoV-2 diagnosis, thus 69 were analyzed. Most of them were hospitalized at Bellaria hospital (N=51), 16 at Sant’Orsola and 2 at Maggiore hospital. The general characteristics of the study population are shown in Table 1 . The median time from symptom onset to hospitalization was 3 days, 77% presented with a moderate/severe pneumonia according to respiratory parameters and imaging study, median of sequential organ failure assessment (SOFA) at hospitalization was 2. During hospitalization, 33% of patients developed severe respiratory failure, 23% were admitted to ICU and 14% were mechanically ventilated. Hydroxychloroquine and LMWH (low-molecular-weight heparin)-based treatments were administered in most patients (88.4%), while tocilizumab was used in 36% of enrolled patients. In-hospital death occurred in 9 patients (13%).

The GM of COVID-19 patients was profiled by 16S rRNA gene-based next-generation sequencing and compared with that of age- and sex-matched healthy Italians, whose sequences are publicly available (Biagi et al., 2016; De Filippis et al., 2016) (see Materials and Methods). A total of 2,792,602 reads (mean ± SD, 40,472 ± 11,243) were obtained and clustered into 12,889 Amplicon Sequence Variants (ASVs).

A significant reduction in alpha diversity was observed in COVID-19 patients compared to healthy controls (p value = 0.0008, Wilcoxon test) ( Figure 1A ). The PCoA of inter-individual variation, based on Bray-Curtis dissimilarity between the genus-level profiles, showed a significant separation between the groups (p value = 0.001, permutation test with pseudo-F ratio) ( Figure 1B ). As regards the composition, the GM of COVID-19 patients showed extensive destructuring even at high phylogenetic levels (see Figure S1 ). It was in fact characterized by alterations affecting the two dominant phyla, Firmicutes and Bacteroidetes, as well as subdominant components, mainly Actinobacteria and Synergistetes (p value ≤ 0.02, Wilcoxon test). At the family level, COVID-19 patients showed an enrichment in generally subdominant taxa, such as Enterococcaceae, Coriobacteriaceae, Lactobacillaceae, Veillonellaceae, Porphyromonadaceae and Staphylococcaceae (p value ≤ 0.03). On the other hand, their GM was characterized by a reduction of the dominant families Bacteroidaceae, Lachnospiraceae and Ruminococcaceae, as well as Prevotellaceae and Clostridiaceae (p value ≤ 0.001) ( Figures 2A and S1 ). At the genus level ( Figure 2B ), as expected, taxa belonging to Bacteroidaceae (i.e., Prevotella and Bacteroides), Lachnospiraceae (i.e., Coprococcus, Blautia, Roseburia and Lachnospira) and Ruminococcaceae (i.e., Faecalibacterium, Ruminococcus, Oscillospira and Anaerofilum) were the main discriminants of the GM of healthy subjects (p value ≤ 0.008). On the other hand, the GM of COVID-19 patients showed a distinctive pattern, with the enrichment of known or potential opportunistic pathogens, such as Enterococcus, Staphylococcus, Serratia and Collinsella, as well as Lactobacillus, Parabacteroides, Lactococcus, Phascolarctobacterium, Odoribacter, Actinomyces, Methanobrevibacter and Akkermansia (p value ≤ 0.02) (see Figure S2 ). According to a species-level analysis, the sequences assigned to Enterococcus, i.e., the dominant genus in the GM of COVID-19 patients (mean relative abundance in COVID-19 patients vs healthy controls, 18.5% vs 0.05%), were mostly attributable to E. faecium (8.4%) along with E. hirae (5.5%), E. faecalis (1.8%) and E. villorum (1.1%) ( Figure 3 ).

PCoA analysis based on Bray-Curtis dissimilarity between GM profiles of COVID-19 patients showed no segregation by age, sex, antibiotic intake in the 2 weeks prior to fecal sampling, length of hospital stay, time interval between fecal sampling and hospital admission, and outcome (death/discharge) (p value > 0.05, permutation test with pseudo-F ratio) (see Figure S3 ). On the other hand, the GM profiles were found to stratify by intensive care unit (ICU) admission and occurrence of bloodstream infection (BSI) (p value ≤ 0.05) ( Figure 4 ). When looking for the bacterial genera driving the clustering patterns, we found that the same taxa were consistently depleted or enriched in relation to the two covariates (p value ≤ 0.001, “envfit” function). Specifically, ICU patients and those developing BSI were characterized by the over-representation (sometimes mono-dominance) of Enterococcus compared to the respective counterparts (p value ≤ 0.001). On the other hand, Streptococcus, Oscillospira, Blautia and other Ruminococcaceae, Lachnospiraceae and Clostridiales taxa were associated with COVID-19 patients who had not entered ICU and those who had not developed BSI (p value ≤ 0.001) ( Figure 4 ). The severe destructuring of ICU and BSI-related GM profile was also evidenced by the further loss of alpha diversity (p value ≤ 0.004) ( Figure 4 ). It should be noted that of the 16 COVID-19 patients who had entered ICU, 14 (87.5%) had developed BSI.

In an attempt to further explore the GM dysbiosis of COVID-19 patients, the sequences in this study were compared with those of critically ill non-COVID-19 patients who had been admitted to the ICU at St. Orsola hospital in October 2019-February 2020 following liver transplantation. Groups were best matched for known microbiota-associated confounding factors (i.e., age and sex ratio) as well as for other therapy and hospitalization-related confounders (i.e., antibiotic intake, time interval between hospitalization and fecal sampling, and development of bacterial infection) (see also Materials and Methods). Regardless of underlying disease, patients admitted to ICU showed comparable alpha diversity levels (p value > 0.05, Wilcoxon test), lower than COVID-19 patients who had not entered ICU (p value ≤ 0.004) ( Figure 5A ). In contrast, the GM structures of the three patient groups (i.e., COVID-19 patients admitted or not to ICU and other ICU patients) clearly segregated in the Bray-Curtis-based PCoA space (p value = 0.001, permutation test with pseudo-F ratio). Several genera were found to drive this separation (p value ≤ 0.001, “envfit” function). Notably, Enterococcus was far overrepresented in both groups of COVID-19 patients compared to other ICU patients, but closely associated with COVID-19 patients admitted to ICU (p value = 0.0001, Wilcoxon test). Similarly, Ruminococcus was enriched in COVID-19 patients but associated with those who had not entered ICU (p value = 0.0003), as well as Oscillospira, Dorea and Coprococcus (p value ≤ 0.01). Critically ill non-COVID-19 patients were mainly discriminated by Enterobacteriaceae genera, particularly Klebsiella (p value ≤ 0.03) ( Figure 5B ).

In order to evaluate the dynamic change in BSI due to Enterococcus spp. (E-BSI) during the COVID-19 pandemic, we compared the E-BSI incidence in critically ill patients over the same period (January-April) since 2017 to 2020 ( Figure 6 ). During the first 4 months of 2020, 1,317 patients were admitted to the ICU with a mean length of stay of 5.3 days and a total patient days of 6,924. The incidence rate of ICU-acquired E-BSI reached 27.4 (95% CI, 1.75-4.29) per 10,000 patient-days during the COVID-19 outbreak. In detail, we observed 19 episodes of E-BSI in the first 4 months of 2020 among patients recovered in ICU, mainly due to E. faecium (57.9%). A significant increase in the incidence rate of E-BSI was observed between 2017 and 2020 (p-value = 0.01, Poisson regression). This incidence rate was 14.8 (95% CI, 0.74–2.96) in 2019 and 15.2 (95% CI, 0.79–2.92) in 2018. The relative risk of ICU-acquired E-BSI during the first 4 months of 2020 was 1.84/3.14-fold higher than in previous years.