The metagenomic and metabolomic profile of the gut microbes in Chinese full-term and late preterm infants treated with Clostridium butyricum

Ethics approval and consent to participate

The study was approved by the Human Research Ethics Committees of the Children’s Hospital of Soochow University (Reference 2020CS017). All specimens were collected according to the guidelines set by the Children’s Hospital of Soochow University. All authors confirm that all methods were performed in accordance with the relevant guidelines and regulations (Declaration of Helsinki). Written informed consents were obtained from the parents of all the enrolled infants before specimen collection.

Subject and selection criteria

The present study aimed to evaluate the composition, abundance, and diversity of gut microbes in Chinese full-term infants and late preterm infants. A total of 44 full-term infants and 50 late preterm infants admitted to the Neonatal Department between November 2020 and April 2021 were enrolled as participants in Suzhou, East China. All the enrolled infants were otherwise healthy without any symptoms of pediatric diseases. The infants’ parents were neither infected nor treated with any antimicrobials.

Detailed infants background and selection criteria could be found in Tables 1, 2, and Supplementary Table 1. The full-term group included 24 girls and 20 boys with gestational ages ranging from 37 to 40 + 5 weeks and birth dates from the first day of life to the 12th day. They are given a mixture of breast milk and formula milk and have an average birth weight of 3152 g.

The late-preterm group included 24 girls and 26 boys with gestational ages ranging from 33 to 36 + 6 weeks and a birth age ranging from 1 to 7 days. There were 21 vaginal births and 29 cesarean births. All were fed Nestle preterm baby milk and had an average birth weight of 2329.2 g. Stool samples were collected from the late-preterm group before probiotic treatment.

Then, the late-preterm newborns were provided probiotics and subsequently placed into two groups: day feeding and evening feeding. The day-feeding group included 12 female newborns and 13 male infants with gestational ages ranging from 33 + 4 to 36 + 2 weeks, with a minimum age of one day and a maximum age of seven days. There were 12 vaginal births and 13 cesarean deliveries. The average birth weight was 2338.4 g.

All the mothers had no history of infection during pregnancy or long-term use of antimicrobials, hormones, or blood products. All infants were free of complications, such as severe infections, cardiopulmonary diseases, congenital genetic diseases, or gastrointestinal malformation. While hospitalized, latamoxef was the only acceptable antimicrobial medication. The term infants are breastfed, while special formula milk for preterm neonates was given, and another nutritional supplementation was administered intravenously.

Group and subgroup

The full-term infants were clustered as Group A (n = 44) whiles late preterm infants were defined as Group B (n = 50). In addition, the 50 late preterm infants were further treated for C. butyricum therapy and stool samples were obtained as Group C. After that the infants from Group C were divided into sub-groups of day-fed C. butyricum (Group C_D, C1-C25, n = 25) and night-fed C. butyricum (Group C_N, C26-C50, n = 25). Figure 1 displays the detailed workflow. Latamoxef was given intravenously after being calculated using the norm of 40 mg/kg/per according to the weight of late-preterm newborns twice a day. The date of Latamoxef consumption started between November 2020 and June 2021, and the average duration of antibiotic treatment was 12 to 13 days. The probiotic of C. butyricum MIYARI 588 was lyophilized powder purchased from Miyarisan Pharmaceutical Co., Ltd. The CFU of C. butyricum is 1X10⁹ per gram. Each newborn was given 500 mL of milk with 40 mg of C. butyricum MIYARI 588 twice a day at 8 a.m. and 3 p.m. for the day group, and 8 p.m. and 3 a.m., for night group, with feces collected one week later.

Study workflow. A total of 144 genomes were enrolled for 16S rRNA metagenomic analyses and fecal metabolome analysis.

Sample

For Group A, the full-term infants were fed without antibiotics or probiotics, and the stool collection was the first stool in hospitalization within 24h, with a collection date between November 13, 2020, and May 21, 2021. For Group B, antibiotics were administered to premature babies hospitalized between October 27, 2020, and June 4, 2021. Feces were collected from preterm newborns between November 2, 2020, and June 7, 2021, without any probiotics consumed. For Group C, probiotics were first given to the preterm infants between November 2, 2020, and June 7, 2021. After probiotic feeding for one week, feces were collected from those infants between November 9, 2020 to June 14, 2021.Additionally, the stool samples were split into groups that received probiotics during the day and those who received them at night, depending on the feeding schedule. Then, each specimen was about 5 g in weight, gathered using a sterile disposable stool sampler and kept in a – 80 ℃ refrigerator for subsequent 16S rRNA sequencing to analyze the composition, abundance, and diversity of the gut microbes. Aseptic precautions were exercised during all procedures.

DNA extraction and 16S rRNA sequencing

After being removed from the sampler, the stool was diluted 5X with molecular grade water and homogenized with a vortex. The stool suspension (250 µl) was used for DNA extraction. The genomic DNA was purified and isolated with MagaBio Soil/Feces Genomic DNA Purification Kit (Hangzhou Bioer Technology, Hangzhou, Zhejiang, China). The extracted DNA was quantified at an A260/A280 nm ratio with the NanoDrop ND-1000 (Thermo Fisher Scientific, Cleveland, OH, USA). Then extraction from each sample was diluted to approximately 5 ng/µl and stored at – 20 °C for 16S rRNA sequencing.

The 12-bp barcoded primers synthesized by Invitrogen (Invitrogen, Carlsbad, CA, USA) were used to amplify the bacterial 16S rRNAV3-V4 fragments. The PCR mixture contained 25 μl reactions Taq (Takara Biotechnology, Dalian, Liaoning, China), 1 μl of each primer (10 mM), and 3 μl DNA (20 ng/μl) template (final volume: 50 μl). The PCR protocol was as follows: 94 °C for 30 s followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 30 s; followed by a final extension at 72 °C for 10 min. The PCR products were subjected to 1% agarose gel electrophoresis and then sequenced on an Illumina Miseq (PE 300) in the MAGIGENE Genomic Institute. The PCR products were mixed according to the Instructions of the GeneTools Analysis Software (Version 4.03.05.0, SynGene).

Bioinformatic analyses of genomic data

Pair-end reads for each sample were merged using the BBmerge module in the BBtools suite v38.94 with options of “loose = t mininsert = 120 mininsert0 = 100 qtrim2 = t qout = 33 entropy = t maxns = 0 trimq = 10”. Low-quality sequences in the merged DNA fragments were removed using the “fastq_filter” module in Vsearch v2.18.0 with a fastq_maxee_rate of 0.02. The “cluster” module in the Vsearch package was applied to the remaining set of high-quality sequences to generate clusters of nearly identical sequences (identities ≥ 0.997) and extract one centroid sequence for each cluster. The same package was applied to remove potential sequencing errors and chimeric sequences from the centroid sequences using the “cluster_unoise” module and the “uchime3_denovo” module, respectively¹⁷. The resulting set of non-repetitive, high-quality sequences was treated as ASVs (Amplicon sequence variants) and subjected to taxonomic assignment based on their comparison results with the SILVA database of 16S rRNA (release 138.1; www.arb-silva.de/) using BLASTn. We kept hits that shared ≥ 85% identities in ≥ 50% of the sequences of the reads. The final taxonomic information of each read was determined based on hybrid LCA criteria¹⁸.

The Quantitative Insights in Microbial Ecology (QIIME) software (Version 1.9.1) was enrolled to calculate Observed-otus, Chao1, Shannon, Simpson, ace, Goods-coverage, and PD_whole_tree index¹⁹. R software (Version 4.1.2) was applied to draw dilution curve, Rank abundance curve, species accumulation curve and analyze the difference between groups of Alpha diversity index. The used T-test and wilcox test were used to analyze the difference between groups of Alpha diversity index. Specifically, the Chao1 estimator and the ACE estimator were used to calculate community richness, while the Shannon index and the Simpson index were applied in community diversity analysis (scikit-bio.org/docs/latest/generated/skbio.diversity.alpha). R software (Version 4.1.2) was used to plot PCA, PCoA and NMDS plots. The stats package and ggplot2 package of R software was used for PCA and PCoA analysis, and the vegan package of R software was used for NMDS analysis. Use R software to analyze the differences between groups of Beta diversity index and perform parametric test and non-parametric test respectively. Stamp analysis was performed using Stamp software, and the filter value of Score was set to default. Metastats analysis uses Mothur software at each classification level (Phylum, Class, Order, Family, Genus, Species) to perform a permutation test between groups to obtain a p-value, and then use the Benjamini and Hochberg False Discovery Rate method to correct the p-value. Anosim, MRPP and Adonis analysis use the anosim function, mrpp function and adonis function of the R vegan package, respectively. AMOVA analysis was performed using the above function of the mothur software. Species with significant differences between groups were analyzed using R software for T-test between groups and plotted.

Fecal metabolomes evaluation

Background of the 40 late-preterm infants enrolled in metabolome analysis was showed in Table 4. The sample stored at − 80 °C refrigerator was thawed on ice. A 400 μL solution (Methanol:Water = 7:3, V/V) containing internal standard was added into 20 mg sample, and vortexed for 3 min. The sample was sonicated in an ice bath for 10 min and vortexed for 1 min, and then placed at − 20 °C for 30 min. The sample was then centrifuged at 12,000 rpm for 10 min (4 °C). And the sediment was removed, then centrifuged the supernatant was at 12,000 rpm for 3 min (4 °C). A 200 μL aliquots of supernatant were transferred for LC–MS analysis²⁰.

Detail LC–MS method was described as follows^{21, 22}. The analytical conditions of LC–MS were as follows, UPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 um, 2.1 mm*100 mm); column temperature, 40 °C; flow rate, 0.4 mL/min; injection volume, 2 μL; solvent system, water (0.1% formic acid): acetonitrile (0.1% formic acid); The column was eluted with 5% mobile phase B (0.1% formic acid in acetonitrile) at 0 min followed by a linear gradient to 90% mobile phase B (0.1% formic acid in acetonitrile) over 11 min, held for 1 min, and then come back to 5% mobile phase B within 0.1 min, held for 1.9 min. The data acquisition was operated using the information-dependent acquisition (IDA) mode using Analyst TF 1.7.1 Software (Sciex, Concord, ON, Canada). The source parameters were set as follows: ion source gas 1 (GAS1), 50 psi; ion source gas 2 (GAS2), 50 psi; curtain gas (CUR), 35 psi; temperature (TEM), 550 °C, or 450 °C; declustering potential (DP), 60 V, or-60 V in positive or negative modes, respectively; and ion spray voltagefloating (ISVF), 5500 V or-4500 V in positive or negative modes, respectively. The TOF MS scan parameters were set as follows: mass range, 50–1000 Da; accumulation time, 200 ms; and dynamic background subtract, on. The product ion scan parameters were set as follows: mass range, 25–1000 Da; accumulation time, 50 ms; collision energy, 30 or-30 V in positive or negative modes, respectively; collision energy spread, 15; resolution, UNIT; charge state, 1 to 1; intensity, 100 cps; exclude isotopes within 4 Da; mass tolerance, 50 mDa; maximum number of candidate ions to monitor per cycle.

All samples were acquired by the LC–MS system machine orders in Metware company, China. The original data file acquisition by LC–MS was converted into mzML format by ProteoWizard software. Peak extraction, peak alignment and retention time correction were respectively performed by XCMS program. The “SVR” method was used to correct the peak area. The peaks with detection rate lower than 50% in each group of samples were discarded. After that, metabolic identification information was obtained by searching the laboratory’s self-built database, integrated public database, AI database and metDNA. For two-group analysis, differential metabolites were determined by VIP (VIP > 1) and P-value (P-value < 0.05, paired t-test). VIP values were extracted from OPLS-DA result, which also contains score plots and permutation plots, and was generated using R package MetaboAnalystR v1.0.1. The data was log transform (log2) and mean centering before OPLS-DA. Identified metabolites were annotated using KEGG Compound database (www.kegg.jp/kegg/compound/), and annotated metabolites were then mapped to KEGG Pathway database (www.kegg.jp/kegg/pathway.html). Significantly enriched pathways are identified with a hypergeometric test’s P-value for a given list of metabolites. The compound name of internal standard used in LC/MS analysis were described in Supplementary Table 2.

Statistical analysis

The analysis of diversity, variance and similarity with repeated measures was used to determine the statistical significance for each experimental group. Statistical analyses were carried out in R (4.0.3). Vegan (v 2.5–6) was used for all alpha-diversity calculations: index of Shannon, Simpson, Chao1, and ACE diversity (alpha diversity measurement of evenness and richness); evenness (homogeneous the distribution of taxa counts) and richness (number of taxa in a community). Pairwise Bray–Curtis dissimilarity was used to assess beta-diversity, or the overall variation between each sample²³. The Bray–Curtis dissimilarity metric compares two communities based on the number or relative abundance of each taxon present in at least one of the communities. The differential metabolites were determined by VIP (VIP > 1) and P-value (P-value < 0.05, paired t-test). Analyses were conducted using GraphPad Prism (v7) with statistical significance accepted as *P < 0.05, **P < 0.01, ***P < 0.005, ****P < 0.001. The images of Barplot, PCA, PcoA, and Stamp were created on Tutools platform (www.cloudtutu.com). The ANOSIM analysis, Violin and Heatmap were plotted by www.bioinformatics.com.cn, a free online platform for data analysis and visualization. Statistical analysis of differential metabolites and MSEA enrichment annotated by KEGG pathways were displayed in Supplementary Tables 3 and 4.