A total of 30 children with GERD and 30 healthy controls (HCs) who were admitted to Beijing Children’s Hospital Affiliated with Capital Medical University from September 2022 to March 2023 were included in this study. The subjects ranged from 3 to 14 years old. The diagnosis of GERD was made using the relevant criteria in the “Pediatric Gastroesophageal Reflux Disease Diagnosis and Treatment Plan (Trial)” (Wang and Jiang, 2006), and the diagnosis of GERD was confirmed by endoscopy or esophageal pH examination over 24 h. HC children were recruited from child health centers. The exclusion criteria were as follows: (1) use of proton pump inhibitors and gastrointestinal mototropic drugs during the week prior to the study; (2) abnormal anatomical structure of the digestive tract; (3) eosinophilic esophagitis; (4) Helicobacter pylori infection; (5) metabolic diseases; (6) coagulation dysfunction; (7) heart, liver, kidney and other organ functional diseases; and (8) blood and immune system diseases. General information (age, sex, body mass index (BMI)), gastrointestinal symptoms, and routine blood test results were collected from the study subjects, and gastroscopy results and 24-hour esophageal pH test results were also collected from patients with GERD. None of the participants had taken antibiotics, probiotics, or prebiotics in the 2 months prior to stool collection. Thirty children with GERD agreed to provide stool samples collected for 16S rRNA sequencing analysis, 28 children agreed to provide serum samples collected for metabolomics detection analysis, 20 children in the control group provided stool samples collected for 16S rRNA sequencing analysis, and 28 children agreed to provide serum samples used for metabolomics detection analysis (see recruitment Flow Chart 1 for details). This study complied with the principles set out in the Declaration of Helsinki. The study protocol was approved by the Ethics Committee of Beijing Children’s Hospital, and all subjects signed the informed consent form (Approval no. [2022]-E-175-Y).

According to the sample collection instructions, feces were collected in a hospital or at home, transported to Beijing Children’s Hospital affiliated with Capital Medical University under low temperature storage, and then stored in a -80°C refrigerator.

Total genomic DNA from samples was extracted using the CTAB method (Yang et al., 2022), and the DNA purity and concentration were measured by agarose gel electrophoresis. Appropriate amounts of sample DNA were diluted to 1 ng/μl in a centrifuge tube, and the diluted genomic DNA was used as a template. Using specific primers with barcodes 341 F (5’- TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG - 3’) and 806 r (5’- GTCTCGTGGGCTCGGAGATGTATATAAGAG ACAGGGACTACHVGGGTWTCTAAT-3’) for 16S bacterial amplicon sequencing (V3-V4 region), the TruSeq^® DNA PCR-Free Sample Preparation Kit was used to construct the libraries. After qualified quantitative analysis by Qubit and Q-PCR, NovaSeq 6000 was used for on-machine sequencing. The reads of each sample were spliced by FLASH, and the resulting splicing sequence was the original tag data. The tag sequences were compared with the species annotation database to detect the chimera sequences, and the final valid data were obtained. The Uparse algorithm was used to cluster all valid data obtained from all samples, and sequences with ≥97% similarity were classified as operational taxonomic units (OTUs). The OTU sequence was annotated by species in the SSUrRNA database.

After the DNA samples were assessed for quality, libraries were constructed using the NEBNext^® Ultra DNA Library Prep Kit for Illumina (NEB, USA). After the DNA samples were assessed for quality, Illumina PE150 sequencing was performed (Treiber et al., 2020). Readfq was used to process the raw data to obtain clean data for subsequent analysis, and the data were assembled and analyzed by MEGAHIT software (v1.0.4-beta). MetaGeneMark was used for gene prediction, CD-HIT software was used to delete redundant genes, DIAMOND software was used for sequence comparison, and MEGAN software’s LCA algorithm was used to determine the species annotation information of the sequences.

Nontargeted metabolomics analysis was performed by LC/MS technology (Kui et al., 2022). Serum samples with volumes of 100 μL were placed in an eppendorf (EP) tube, 400 μL of 80% methanol solution was added, vortex shocking was carried out, the samples were placed in an ice bath for 5 min, and centrifugation was performed at 15,000 × g and 4°C for 20 min. The supernatant was diluted to achieve a 53% methanol content, and the supernatant was collected for subsequent detection after centrifugation. Ultrahigh-performance liquid chromatography–tandem mass spectrometry (UHPLC–MS/MS) analysis was performed using the Vanquish UHPLC system (Thermo Fisher, Germany) in conjunction with the Orbitrap Q Exactive™ HF Mass spectrometer (Thermo Fisher, Germany). The chromatographic separation conditions were as follows: chromatography was performed on a Hypesil Gold column (C18) at a column temperature of 40°C and a flow rate of 0.2 mL/min. In positive mode, mobile phase A was 0.1% formic acid, and mobile phase B was methanol; in negative mode, mobile phase A was 5 mM ammonium acetate, and mobile phase B was methanol. The detection parameters were as follows: spray voltage, 3.5 kV; sheath gas flow rate, 35 psi; auxiliary gas flow rate, 10 L/min; ion transfer tube temperature, 320°C; ion import RF level, 60; auxiliary gas heater temperature, 350°C; polarity, positive ion and negative ion mode; and MS/MS secondary scanning, data-dependent scanning. The original data were imported into the CD 3.1 library search software for processing, with a mass deviation of 5 ppm, signal strength deviation of 30%, signal-to-noise ratio of 3, and minimum signal strength; ion and other information were used for peak extraction, and quantification of the peak area and integration of the target ion were performed. Then, the molecular formula was predicted using the molecular ion peaks and fragment ions, which were compared with the mzCloud, mzVault and Masslist databases. Finally, the original quantitative results were standardized to obtain the identities and relative levels of metabolites.

SPSS 23.0 was used for statistical analyses of all demographic and clinical data. The quantitative demographic and clinical data with a normal distribution are expressed as the mean ± standard deviation, and the differences between groups were tested by t tests. Quantitative demographic and clinical characteristics with nonnormal distributions are expressed as the medians (P25, P75), and differences between groups were tested using the Wilcoxon rank sum test. Qualitative demographic and clinical data are expressed as percentages. The chi-square test was used to test the differences between groups. The results with p<0.05 were considered statistically significant.

Quantitative Insights Into Microbial Ecology (QIIME) software (Version 1.9.1) was used to calculate the Shannon index, Simpson index and Wilcox test for intergroup difference analysis. Based on the unweighted UniFrac β diversity analysis, principal coordinates analysis (PCoA) visualizations were analyzed using the ggplot2 software package R (Version 2.15.3), and the differences at the phylum, genus and species levels and the functional module abundance between the two groups were analyzed using the Wilcoxon rank sum test. The p value was corrected by the Benjamini–Hochberg method to obtain the q value to control the error rate of multiple tests. Linear discriminant analysis effect size (LEfSe) software was used for LEfSe analysis. DIAMOND software was used to compare the nonredundant gene set with the Kyoto Encyclopedia of Genes and Genomes (KEGG), eggNOG and CAZy databases to obtain the corresponding functional information of gene sequences.

Principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed after the metabolite data were transformed by MetaX software. The analysis of intergroup metabolite differences was performed using a t test, and the fold change (FC) between groups was calculated. The default criteria for differential metabolite screening were variable importance of projection (VIP)>1, P <0.05, FC≥2 or FC ≤ 0.5. The ggplot2 package of R software was used to draw the volcano map, and the Pheatmap package was used to draw the clustering heatmap using the z score to normalize the metabolite data. The identified metabolites were annotated by the KEGG, HMDB and LIPID Maps databases. The KEGG metabolic pathways enriched by differential metabolites were identified, and results with P <0.05 were considered to indicate significant enrichment. Spearman correlation coefficients were calculated using the levels of the differentially enriched bacteria and metabolites in R software.

There were no statistically significant differences in age, sex, or BMI among all participants. Details of the demographic and clinical characteristics of the children with GERD and HCs are shown in Supplementary Table S1 . To determine whether the gut microbiota composition of children with GERD was changed as compared to HC, 16S rDNA sequencing was used to detect fecal microbes in children with GERD and HC children. The species Venn diagram showed a total of 834 OTUs in the two groups, with 2250 OTUs specific to the GERD group and 595 OTUs specific to the control group ( Figure 1A ). The α diversity results showed that the Shannon and Simpson indices in children with GERD were significantly lower than those in the HC group, indicating a decrease in gut microbiota diversity in children with GERD ( Figures 1B, C ). Analysis of similarities (ANOSIM) based on unweighted UniFrac distance showed significant differences in the distribution of the microbial community structure between the two groups ( Figure 1D ). PCoA and nonmetric multidimensional scaling (NMDS) showed a difference in the gut microbiota composition between the two groups ( Figures 1E, F ). At the phylum level, the relative abundances of Proteobacteria and Bacteroidetes in the GERD group were significantly higher than those in the HC group, while the relative abundances of Firmicutes and Actinobacteria were significantly lower ( Figure 1G ). At the family level, the relative abundances of Enterobacteriaceae, Bacteroidaceae and Prevotellaceae in children with GERD were high, and those of Bifidobacteriaceae, Ruminococcaceae, and Lachnospiraceae were low ( Figure 1H ). At the genus level, the relative abundances of Bacteroides, Prevotella-9, Escherichia-Shigella, Klebsiella and Haemophilus in children with GERD were significantly higher than those in the HC group. The relative abundances of Bifidobacterium, Faecalibacterium and Streptococcus were lower in the GERD group than in the HC group ( Figure 1I ). LEfSe analysis was used to process the samples to identify species with significant differences between the groups. The results showed that Prevotella-9 and Bacteroides had the highest abundances in the GERD group, while Bifidobacterium and Blautia had the highest abundances in the HC group ( Figures 1J, K ).

To determine the changes in the gut microbiota at the species level in children with GERD, metagenomic sequencing was performed using stool samples from 10 children with GERD and 10 age-matched HC children after sample selection through microPITA. The Venn diagram results showed that the two groups had a total of 858,667 genes, with 278,228 specific genes in the GERD group and 178,342 specific genes in the HC group ( Figure 2A ). The PCoA results based on the Bray–Curtis distance showed significant separation of the microbial composition at the species level between the GERD and HC groups ( Figure 2B ), and ANOSIM confirmed the significant differences in the gut microbiota composition at the species level between the two groups ( Figure 2C ). The clustering heatmap of bacterial species enriched by differences between the GERD and HC groups is shown in Figure 2D . Figures 2E, F show the top 20 species of bacteria that were significantly enriched in the gut microbiota of the GERD group and HC group, respectively. The abundances of Bacteroides stercoris, Bacteroides vulgatus, Alistipes putredinis, Bacteroides dorei, and Bacteroides fragilis increased significantly, while those of Fusicatenibacter saccharivorans, Blautia wexlerae, Blautia obeum, Anaerostipes hadrus and Eubacterium hallii decreased significantly.

To identify the functional characteristics of the gut microbiota of children with GERD, this study further annotated the gut microbiota metagenomic data using the KEGG database, and ANOSIM indicated that the gut microbiota differed significantly between the GERD group and the HC group at the KEGG Orthological (KO) level ( Figure 3A ). The PCoA results confirmed that there were significant differences in the distribution of microbial functions based on the KEGG module between the GERD and HC groups ( Figure 3B ). The activities of metabolic pathways such as energy, amino acids, vitamins, carbohydrates and lipids were lower in the GERD group than in the HC group, while the activities of glycan synthesis and metabolic pathways were higher ( Figures 3C–E ).

The eggNOG orthologous group (og) in the GERD group was also significantly separated from that in the HC group ( Supplementary Figure S1A ). The PCoA results showed significant differences in the distribution of functional genes between the two groups ( Supplementary Figure S1B ). The functions of cell cycle control, cell division, cytoskeleton, extracellular structure, nucleotide transport and metabolism, energy production and conversion, etc., were more impaired in the GERD group than in the HC group, while cell motility, cell wall/membrane/envelope biogenesis, inorganic ion transport and metabolic pathway enrichment were affected ( Supplementary Figures S1C–E ).

The above metagenome-based functional annotation analysis revealed functional changes in the gut microbiota in children with GERD. In this study, the serum samples of the two groups of children were used for metabolomic analysis by LC/MS technology. Principal component analysis (PCA) and PLS-DA of serum metabolites were performed. The PLS-DA results showed significant differences in the metabolic characteristics of serum samples between the GERD group and the HC group ( Figures 4A–D ). The PLS-DA model verified that there was no overfitting phenomenon, and the sample characteristics could be used for subsequent analysis ( Supplementary Figure S2 ). Furthermore, the differential metabolites between the two groups were analyzed by volcano plot, and the threshold values were set as VIP > 1.0, FC > 1.5 or FC < 0.667 and P value< 0.05. The results showed that 182 differential metabolites were identified in the positive ion mode between the GERD and HC groups, of which 130 metabolites had increased levels in the GERD group. The levels of 52 metabolites were lower in the GERD group. In negative ion mode, 106 differential metabolites were identified between the groups, of which 49 had higher levels and 57 had lower levels in the GERD group ( Figures 4E, F ).

The results of the clustering heatmap showed that samples obtained from the same tissue were clustered into a group, and the relative content of metabolites in different sample groups significantly differed ( Figure 5A ). According to the differential multiple value of the metabolites in the two groups, logarithmic conversion was performed with the base of 2, and the top 20 metabolites in this dataset were selected and displayed in the matchstick chart. As shown in Figure 5B , the levels of metabolites such as 4-methoxycinnamaldehyde, 12-epi leukotriene B4, propionyl-L-carnitine, prostaglandin A1, and prostaglandin G2 in the GERD group were significantly increased. The levels of metabolites such as 4-benzoic acid, gentisic acid, hydroquinone, tetradecanediacid, 2-phenylglycine and theophylline decreased significantly. Further KEGG annotation of the differential metabolites showed that the top three enrichment pathways were global and overview maps, amino acid metabolism, and lipid metabolism ( Figure 5C ). The KEGG bubble map shows the top 20 pathways with the most significant enrichment. The metabolic pathways affected by the differences in metabolites between children with GERD and HCs mainly included arachidonic acid (AA) metabolism; tyrosine metabolism; sulfur metabolism; glycine, serine, and threonine metabolism; glutathione metabolism; and caffeine metabolism ( Figure 5D ).

Pearson correlation analysis was used to analyze the associations between the differential gut microbes in children with GERD and serum metabolites, and the results showed that the changes in serum metabolites in children with GERD were significantly correlated with the changes in the levels of various bacteria in the gut microbiota. The levels of Chondromyces, Labilithrix and Stigmatella were significantly positively correlated with the levels of a variety of metabolites ( Figure 6 ), suggesting significant interactions between the gut microbiota and metabolites that may affect the occurrence of GERD.