Informed consent of all participants was obtained for this study. The collection of stool samples and data analysis were approved by the institutional ethics board of the Sixth Affiliated Hospital, Sun Yat-sen University (NO.2023ZSLYEC-231). All interventions performed on patients were based on clinical judgment, according to the needs of the patients and not for the purposes of this study. Data were analyzed anonymously. We confirm that all methods are implemented in accordance with relevant guidelines and regulations. The study was also conducted in accordance with the Declaration of Helsinki (revised in Fortaleza, Brazil, October 2013).

In January and February 2023, COVID-19 patients were categorized into two groups, namely the EN group and the PN group, based on their clinical treatment requirements. Throughout their hospitalization, all participants received personalized nutritional support treatment as per the standard of care. Prior investigations have indicated a substantial alteration in the intestinal flora of patients after a week of nutritional support therapy (20, 21). Consequently, in our study, all patients were administered enteral nutrition or parenteral nutrition support for a duration of one week, following which stool samples were obtained from COVID-19 patients. Eligibility criteria for patients included being 60 years of age or older, according to the Chinese definition of the elderly individuals (22), and tested positive for SARS-CoV-2 by nasopharyngeal swabs using quantitative RT-PCR according to the recommendations of the National reference laboratory and the National Health Bureau. This study employed a prospective non-interventional design to investigate the patient inclusion process, as depicted in Figure 1 . Within 48 hours of admission, each group received distinct forms of nutritional support. The EN group had access to oral nutrition powder, nasogastric tube, or nasojejunal access, while the PN group received parenteral nutrition via a central venous catheter (subclavian or jugular vein). The decision was made following a clinical evaluation of bowel function, with no rigid criteria imposed, thereby allowing the individual clinician’s discretion. Consequently, additional factors were taken into account, including the patients’ underlying condition and the anticipated rate of postoperative or intervention-related recuperation. The recommended daily intake of enteral nutrition was determined as 20-25 kcal/kg, while protein requirements were estimated at 1.5 g/kg (23, 24). The caloric intake of parenteral nutrition was ascertained to be 20-25 kcal/kg per day, while the calorie/nitrogen ratio was computed to be 120-150:1 (23). Glucose accounts for 50% to 70% of the total energy requirements, with lipid utilization contingent upon serum triglyceride levels. Furthermore, the intravenous solution was supplemented with sufficient quantities of vitamins, electrolytes, trace elements, and insulin. The medical records contained comprehensive information regarding COVID-19 patients, encompassing routine demographic data, past medical history, epidemiological exposure history, pre-existing comorbidities (such as diabetes, hypertension, chronic respiratory diseases, immunosuppression, hematological oncological disease, and previous chronic therapy), symptoms, signs, laboratory results, chest x-rays, administered antiviral therapy, corticosteroids, rescue invasive treatment, intensive care units (ICU) care, supportive measures, and in-hospital complications. Furthermore, clinical outcomes such as the duration of mechanical ventilation, ICU length of stay, and 60-day mortality were also documented.

Stool samples were collected from elderly COVID-19 patients as soon as possible after admission and immediately refrigerated and thawed within 24 hours and frozen at -80 °C until nucleic acid extraction. Stool samples, whether transported short or long distances, must be kept in containers filled with liquid nitrogen and supervised by humans. All samples were collected during standard care rounds using all necessary precautions. The Clinical records and laboratory information system (LIS) were used to retrieve patient data, including laboratory test results and clinical presentation.

Stool DNA was tested for concentration, integrity and purity before library construction, and the qualified samples enter the process of library construction. The gut microbiota composition was determined by 16S rRNA gene sequencing in fecal samples collected from the participants. Firstly, the sample was interrupted, and 1μg metagenomic DNA was taken, which was interrupted with a Covaris (Covaris S2 System, Massachusetts, USA) ultrasonic fragmenter. The interrupted sample magnetic beads were subjected to fragment selection so that the sample bands were concentrated around 200 to 400 bp. Then the end is repaired, “A” is added, and the connector is connected. The reaction system was prepared, the reaction was suitable for A certain period of time, the double stranded cDNA end was repaired, and the A base was added to the 3’end, the linker was prepared, the reaction system was suitable for a certain period of time, the linker was connected with the DNA. The PCR reaction system was prepared, and the reaction program was set to amplify the ligation product. The amplified products were purified and recovered by magnetic beads. After denaturation of the PCR product to a single strand, the cyclized reaction system was prepared, and the reaction was thoroughly mixed at the appropriate temperature for a certain time to obtain the single stranded circular product. After digestion of the uncyclized linear DNA molecules, the final library was obtained. A single-stranded circular DNA molecule replicates through a rolling circle to form a DNA nanosphere (DNB) containing multiple copies. The obtained DNBs were added to the mesh holes on the chip by high-density DNA nanochip technology, and sequenced by combined probe anchor polymerization (cPAS). The metagenome was sequenced using the MGISEQ 2000 platform PE150. SOAPnuke (v1.5.0) (25) software was used to filter the raw data for quality control, and Bowtie2 (2.2.5) (26) was used to align the host sequence and remove the sequence in the ratio to generate Clean Data. MEGAHIT (v1.2.9) (27) was used for de novo assembly of Clean Data that passed quality control. Assembled sequences less than 200 bp in length were filtered out. Species annotation was performed using Kraken2’s default parameters, while species-level abundance was estimated for metagenomic samples using a Bayesian algorithm and Kraken classification results via Bracken. UHGG database of human gastrointestinal genome was used for database selection (28). The sequence number of species contained in the sample was calculated by comparison between Kraken2 and self-built database (screened UHGG database), and then Bracken2 was used to estimate the actual abundance of species in the sample and complete the species annotation. The R package was used to calculate the Alpha diversity of species, including chao1 index, shannon index and simpson index. It also measures the difference between samples or groups by calculating Bray-Curtis distance and Jsen-Shannon divergence, called Beta diversity, to reflect whether there are significant microbial community differences between sample groups. Besides, linear discriminant analysis effect size (LEfSe) method was used to identify differentially rich taxa, and the LDA threshold was set to be >4 (29). Methods Nonparametric test and linear discriminant analysis were combined to find biomarkers of each group. LEfSe searched for the biomarker function of each group (LDA> threshold function, with higher abundance in the corresponding group and lower abundance in other groups). That is, functions that are significantly more abundant in this group than in the other groups.

Descriptive statistics are presented as numbers and percentages for categorical variables, as means and standard deviations (SD) for continuous variables, and as medians of interquartile intervals (IQRs) if continuous variables are not normally distributed. Parametric tests (Student’s t-test and one-way ANOVA) and nonparametric tests (Mann-Whitney and Kruskal-Wallis tests) were used, as appropriate, taking into account the assumption of normality and the number of groups compared. Bonferroni tests were performed to correct for multiple testing. The Mann-Whitney test was used to measure the significance of metagenomic data. The Kolmogorov-Smirnov test was used to test the assumption of normality of the distribution of variables. Significant predictors of 60-day mortality (P<0.1) were used in a multiple regression model integrating odds ratio (OR) with 95% confidence interval (CI). The forward stepwise method was used, and the goodness of fit was tested by the Hosmer-Lemeshow test. For categorical variables, the chi-square test and Fisher’s exact test were used as appropriate. The Alpha diversity and beta diversity was measured with the principal coordinate analysis (PCoA) using Bray-Curtis and unweighted unifrac dissimilarity matrices. The dysbiosi index for the EN vs PN group was calculated using the selbal algorithm (in R-version 4.2.2), which takes into account the constitutive nature of the microbiota data (30, 31). The moderating effect of the association between the ecological imbalance index (as a continuous variable) and mortality was further assessed in a Cox regression, adjusting for covariates of differences in trend levels. Spearman correlation analysis was performed to determine the relationship between fecal microbiome abundance and clinical characteristics. The corrplot package is a graphical display of a correlation matrix (R-4.2.2 for Windows). All statistical analyses were conducted by using SPSS software (version 26.0), and statistical significance was defined as a P value<0.05.

The non-redundant genes were functional annotated using the blastp function of Diamond (0.8.24) (32), and then the intestinal microbiota functional analysis was performed. Functional metabolic pathway enrichment analysis was performed using Kyoto Encyclopedia of Genes and Genomes (KEGG) and the cluster of orthologous groups of proteins (COG) database based on the predicted metagenome. In addition, Statistical Analysis of Metagenomic Profiles (STAMP) was employed to analyze the number of marker genes assigned to a functional profile indicating the number of sequences assigned to different biological subsystems or pathways (33). The predicted KEGG and COG pathways were compared between groups using STAMP, and the P value of significance (by factorial Kruskal-Wallis test) was 0.05.

Table 1 presents the key attributes, clinical variables, and outcome measures associated with elderly individuals affected by COVID-19. The study encompassed a cohort of 71 older adults who were confirmed positive for SARS-CoV-2 through laboratory testing. The research was conducted between January 2023 and February 2023, and the participants were assigned to receive either EN or PN as a means of nutritional support. Within the EN group, there were 19 females and 29 males, with an average age of 77.67 ± 9.12 years. In contrast, the PN group comprised of 4 females and 19 males, with a mean age of 80.43 ± 8.29 years. There were no statistically significant differences in terms of age and gender between the two groups, as shown in Table 1 . Regarding comorbidities, the EN group consisted of 17 patients (35.4%) with diabetes, 23 patients (47.9%) with hypertension, 21 patients (6.3%) with chronic respiratory disease, 6 patients (12.5%) with chronic cardiac disease, and 5 patients (10.4%) with a history of cancer, as indicated in Table 1 . Additionally, the PN group demonstrated statistically significant reductions in PaO2/FiO2 ratio, respiratory index, and resting arterial oxygen saturation (P<0.05). Conversely, white blood cells (P=0.032), neutrophils (P=0.020), and procalcitonin (P<0.001) were significantly elevated in the PN group compared to the EN group. Furthermore, there was no significant difference in antibiotic exposure between the two groups, as 33 patients (68.8%) in the EN group received antibiotics during hospitalization, which was not significantly different from the PN group (P=0.099). A significant difference was not observed between the groups when it came to Covid-19-specific therapies such as steroids and intravenous immunoglobulins ( Table 2 ). The patients in the PN group were significantly more likely to undergo ICU admission (P<0.001), require invasive ventilation (P<0.001), suffer respiratory failure, shock, complications, organ failure, and die within 60 days (P=0.001). In the present study, a cohort of 11 patients was admitted to the ICU exclusively from the PN group for therapeutic purposes. The duration of their stay in the ICU was recorded as 12.00 days [IQR, 8.00-32.00 days]. The comprehensive breakdown of complications in the elderly patients diagnosed with COVID-19 pneumonia can be found in Table 3 .

The 60-day mortality rate was found to be 2.1% (1/48) in the EN group and 30.4% (7/23) in the PN group (P=0.001). To determine the independent predictors of 60-day mortality, a stepwise logistic regression analysis was conducted, incorporating age (as a continuous variable), sex (male or female), presence of concomitant diseases (yes or no), and body mass index (also continuous) as variables. Additionally, univariate logistic regression analysis revealed that nine other variables significantly influenced the 60-day mortality rate. The variables included in the analysis were nutrition type (enteral, parenteral), systolic blood pressure, diastolic blood pressure, respiratory index, invasive ventilation (no, yes), hemoglobin count (continuous), white blood cell count (continuous), neutrophil count (continuous), and creatinine count (continuous). Consequently, these parameters were also incorporated into the multiple regression analysis. To assess the adequacy of the logistic regression model, the Hosmer-Lemeshow statistic was employed. Multivariate regression analysis showed that high respiratory index [odds ratio (OR): 1.835, 95% confidence interval (CI): 1.181-2.851], invasive respiration [OR: 0.012, 95% CI: 0.001-0.121], low hemoglobin [OR: 0.947, 95%CI: 0.911-0.984], and PN support [OR: 0.049, 95%CI: 0.006-0.426] was a predictor of mortality. These models are both illustrated in Table 4 . To rule out the possibility that the choice of nutritional pattern is based on disease severity, we have established a multiple regression model that takes into account all confounding factors that may affect disease severity. Our results showed that the choice of nutritional approach was not significantly associated with the severity of illness at admission, but was associated with organ failure and complications when using regression analysis with baseline information at admission ( Table 5 ).

The presence of SARS-CoV-2 viral RNA in the gastrointestinal tract leads to infection of intestinal cells and the duodenum, indicating its impact on the gut microbiome’s composition. To evaluate this, 16S ribosomal RNA (rRNA) gene analysis was conducted on 71 fecal samples to determine the microbial composition. Changes in the relative abundance (median) of each bacterial genus were assessed to examine alterations in the gut microbiota composition among COVID-19 patients in both the PN and EN groups. Notably, no significant difference in bacterial gene alpha diversity was observed between the PN and EN groups, as depicted in Figure 2A . Next, we proceeded to analyze the Alpha diversity of the species in both the EN group and PN group. The Chao1, Shannon, and Simpson indices were employed as commonly utilized measures of Alpha diversity. Notably, the analysis of the Shannon diversity index revealed no statistically significant difference between the two groups (P>0.05, Wilcoxon rank sum test) ( Figure 2B ).

To ascertain potential dissimilarities in the overall gut microbial composition between the two groups, we investigated various markers of beta diversity. The differences in fecal microbial communities between the two groups (EN and PN) were visualized using Intestinal microbiota PCoA and Non-metric multi-dimensional (NMDS) based on Bray-Curtis difference scaling ( Figures 2C, D ). Additionally, Analysis of similarities (Anosim) was employed to examine the dissimilarities between the two groups. The findings indicated that while there were resemblances between the two groups, no statistically significant differences were observed (P=0.167) ( Figure 3 ). Subsequently, an examination was conducted to assess the alterations in the relative abundance of intestinal microflora among COVID-19 patients receiving varying forms of nutritional support at the phylum level. Specifically, the two patient cohorts primarily consisted of Bacteroidetes, Firmicutes, and Proteobacteria at this taxonomic level. Notably, these distinctions did not reach statistical significance ( Figure 4 ). However, employing STAMP analysis, it was observed that the EN group exhibited a higher relative abundance of Chloroflexi member in comparison to the PN group among elderly COVID-19 individuals (mean 0.58% vs 0.31%, P=0.043). In COVID-19 individuals, the abundance of Nucleocytoviricota members in the PN group was significantly higher compared to the EN group (mean 0.34% vs. 0.54%, P=0.033) ( Figure 5 ). Regarding bacterial composition, Enterococcus_faecium, Escherichia_coli, Phocaeicola_vulgatus, Bacteroides_thetaiotaomicron, Bacteroides_fragilis, and Bacteroides_Fragilis were identified as the most prevalent bacterial genera in both groups ( Figure 6 ).

To further examine alterations in bacterial abundance, the researchers employed linear discriminant analysis effect size (LEfSe) to compare the variations in gut microbiota between the two groups of COVID-19 patients. The LEfSe analysis successfully identified several taxa at various levels, which exhibited dissimilarities among elderly COVID-19 patients receiving distinct nutritional support therapies. LEfSe feature selection using LDA scores greater than 3.0 successfully identified a number of bacterial genera that exhibited discriminatory characteristics between subjects in the EN group and PN group. Notably, the elderly patients with COVID-19 in the PN group demonstrated a significant increase in Enterococcus faecium, Pseudomonas aeruginosa, and Corynebacterium striatum compared to patients in the EN group. Conversely, the EN group exhibited significant increases in Ruminococcus, Negativicutes, Bacteroides, and Blautia ( Figure 7 ).

For further analysis, we generated the dysbiosis index for EN vs PN group based on the R package selbal, which was used to select bacterial taxa that differentiated patients in the EN group from those in the PN group. This index is defined as the difference between arithmetic means of log-transformed abundances in these two groups of taxa (30, 31). As shown in Figure 8 , this index was associated with nutritional support therapy (P<0.001). We further performed Cox regression analysis of the dysbiosis index with the covariates that showed significant differences between the EN group and PN group in Table 2 . The results show that the clinical outcomes of patients in the EN group and the PN group, and separately adjusting for markers with different trend levels between the two groups did not change this association ( Figure 9 ). To further analyze the relationship between gut microbiota and related confounding factors in elderly COVID-19 patients, including the use of antibiotics or other drugs, we correlated the dysbiosis index of EN vs PN group with baseline data. We further grouped elderly COVID-19 patients with or without antibiotic or dexamethasone treatment and calculated the relationship with the ecological disequilibrium index. Our results showed that patients receiving antibiotic and immunoglobulin treatment had higher ecological disequilibrium indexes, while patients receiving dexamethasone treatment had lower ecological disequilibrium indexes, but without statistical difference (P>0.05) ( Table 6 ). In addition, there was no significant correlation between the ecological imbalance index and the severity of disease in patients ( Table 7 ).

The LEfSe analysis was employed to assess the disparity in gut microbiota between the two cohorts of COVID-19 patients. In comparison to the EN group, the PN group exhibited a noteworthy elevation in the relative abundance of Enterococcus_faecium and Pseudomonas_aeruginosa, whereas the relative abundance of Ruminococcus was significantly diminished. To ascertain the potential association between the distinct microbial compositions in the two cohorts and the clinical indicators, we conducted a Spearman correlation analysis. This analysis aimed to investigate the interplay between the gut microbiota and the clinical indicators in both groups of COVID-19 patients. The findings revealed a positive correlation between Enterococcus_faecium and Organ failure (r=0.285, P=0.016), as well as White blood cell parameter (r=0.428, P<0.001). The Neutrophil cell parameter (r=0.455, P<0.001), Alanine aminotransferase parameter (r=0.244, P=0.042), and Creatinine parameter (r=0.292, P=0.013) exhibited positive correlations. Additionally, the genus Pseudomonas_aeruginosa demonstrated a positive correlation with complications during hospitalization (r=0.281, P=0.018). Furthermore, the presence of Ruminococcus exhibited a negative correlation with the incidence of Shock (r=-0.236, P=0.048), Organ failure (r=-0.259, P=0.029), and ICU transition (r=-0.239, P=0.044). Additionally, it displayed a negative correlation with the Invasive ventilation parameter (r=-0.242, P=0.042) and complications during hospitalization (r=-0.259, P=0.029) as indicated in Table 8 . To visually represent the relationship between the relative abundance of three distinct flora and each clinical indicator, a heat map was employed ( Figure 10 ).

The perturbation of the metabolic activity of the gut microbiota serves as an indicator of the severity of disease in individuals affected by COVID-19. Additionally, alterations in the relative prevalence of significant anticipated genomic functions were observed between the enteral and PN groups within the gut microbiota. Notably, the metabolic signatures identified within the same KEGG pathway exhibited variations across samples from patients receiving distinct nutritional support modalities, implying a potential correlation with early colonization and the progression of the disease. Changes in the metabolic profile and energy generation of immune cells have the potential to impact their activation, consequently worsening the imbalanced immune response observed in individuals with COVID-19. Additionally, we conducted STAMP analysis to identify functional disparities between the studied groups. Notably, genes responsible for ribosome synthesis, starch and sucrose metabolism, phosphotransferase system, ascorbate and uronate metabolism, as well as microbial metabolism in diverse environments, exhibited higher expression levels in the intestinal microbiome of patients in the PN group when compared to those in the EN group (P<0.05, Figure 11 ).