Urinary metabolomic profiles uncover metabolic pathways in children with asthma

Abstract

Objective

The prevalence of asthma has gradually increased worldwide in recent years, which has made asthma a global public health problem. However, due to its complexity and heterogeneity, there are a few academic debates on the pathogenic mechanism of asthma. The study of the pathogenesis of asthma through metabolomics has become a new research direction. We aim to uncover the metabolic pathway of children with asthma.

Methods

Liquid chromatography (LC)–mass spectrometry (MS)–based metabolomic analysis was conducted to compare urine metabolic profiles between asthmatic children (n = 30) and healthy controls (n = 10).

Results

Orthogonal projections to latent structures-discrimination analysis (OPLS-DA) showed that there were significant differences in metabolism between the asthma group and the control group with three different metabolites screened out, including traumatic acid, dodecanedioic acid, and glucobrassicin, and the levels of traumatic acid and dodecanedioic acid in the urine samples of asthmatic children were lower than those of healthy controls therein. Pathway enrichment analysis of differentially abundant metabolites suggested that α-linolenic acid metabolism was an asthma-related pathway.

Conclusions

This study suggests that there are significant metabolic differences in the urine of asthmatic children and healthy controls, and α-linolenic acid metabolic pathways may be involved in the pathogenesis of asthma.

Keywords:

• Asthma
• metabolomics
• α-linolenic acid metabolism
• children
• urine

Introduction

Asthma is a complex heterogeneous disease and is the most common chronic disease in childhood (1–3). With the progress of urbanization and industrialization over the years, the prevalence of asthma has increased year by year, affecting approximately 14% of children worldwide and 8.6% of children in the United States, which has now made asthma a global public health problem (4,5). At the same time, according to the results of an epidemiological survey about pediatric asthma in China, the 2-year prevalence rate of asthma in children aged 0–14 years increased from 0.91% in 1990 to 1.54% in 2000 to 2.32% in 2010, and the total prevalence rate increased from 1.09% in 1990 to 1.97% in 2000 to 3.02% in 2010 (6,7), which means that the number of children with asthma in China is constantly increasing. The rising prevalence aggravates the burden on individuals, families, and society, especially in low- and middle-income countries, making the prevention and treatment of asthma an urgent public health issue (8).

Asthma is a heterogeneous disease with a wide range of variability at the genetic level and is subject to gene–environment interactions that allow different patients to differ in disease severity, natural disease course, and response to available treatments (9). Because of the difficulty of early diagnosis and definitive treatment in asthma, we can elucidate the pathological mechanism of asthma at the molecular level by identifying potential diagnostic biomarkers/therapeutic targets to achieve early diagnosis and early treatment (10). Metabolomics relies on advanced analytical chemistry technology and complex statistical methods to systematically analyze the dynamic changes in biological specimens (such as biological fluids, cells, and tissues), comprehensively evaluate the impact of small molecule metabolites due to disease progression or environmental changes, and then explore the feasibility of its application in pathogenesis research and clinical diagnosis (11,12). In recent years, metabolomics which aims at identification and quantitation of small molecules involved in metabolic reactions, provides a comprehensive assessment of biological processes in the development of childhood asthma. A multicenter prospective cohort study showed that the metabotype of infants with severe bronchiolitis characterized by high abundance of inflammatory amino acids and low abundance of polyunsaturated fatty acids had the highest risk for developing asthma (13). One study showed that potential protective prenatal nutrients (vitamin D, n-3 polyunsaturated fatty acids) may alter risk of asthma and allergic disease into adolescence (14). Based on the existing research results, we hypothesize that fatty acid metabolism plays an important role in the development of childhood asthma. Because of its noninvasive and convenient sampling method, urine is a widely used body fluid in metabolomic studies (15). For example, Saude et al.’s study provided evidence that urine metabolomics could be used to distinguish between healthy controls, stable asthma patients, and patients with asthma exacerbations (16). Taking into consideration all these issues, our study further systematically studied the metabolic mechanism of children with asthma by liquid chromatography (LC)–mass spectrometry (MS)-based untargeted metabolomics analysis, which is popular in metabolomic studies due to its high throughput, soft ionization, and good coverage of metabolites, on urine samples of children with mild to moderate asthma in China and healthy controls. The findings of our study on metabolic pathways in children with asthma may provide new therapeutic target for childhood asthma.

Methods

Study population

Thirty children with mild to moderate asthma and 10 healthy controls aged 5–10 years were enrolled from Deqing People’s Hospital between August 2022 and January 2023. The clinical manifestations of children with asthma met the following diagnostic criteria: (1) four or more cumulative episodes of recurrent coughing, wheezing, shortness of breath, and chest tightness; (2) wheezing sounds or prolonged exhalation in the lungs could be examined during the attacks; (3) previous lung function indicated a decrease in forced expiratory volume in the first second (FEV₁), and the diastolic test was positive; (4) the symptoms of wheezing improved after controlled treatment. All asthmatic children were in stable condition by using inhaled corticosteroids as a maintenance medication. The control group consisted of healthy children without wheezing, and those with a history of related respiratory diseases were excluded. All participants were free of respiratory infections and had no special dietary arrangements before sampling. At recruitment, children underwent routine blood tests, lung function assessments, and total immunoglobulin E (IgE) measurements, and their demographic and atopic index data were collected. The Ethics Committee of Deqing People’s Hospital approved this study (LL2021-K21). All parents signed written informed consent before sample and data collection. All methods were carried out in accordance with the principles of the Helsinki Declaration.

Urine sample collection and preparation

One milliliter of early morning midstream urine specimens was dispensed into centrifuge tubes, 10 µL of sodium azide working solution (0.5 mg/L) was added to inhibit bacteria, and the samples were stored at −80 °C until analysis. After being thawed at room temperature, 100 µL of sample was transferred to an Eppendorf tube. After the addition of 400 µL of extract solution (methanol:acetonitrile = 1:1, containing isotopically labeled internal standard mixture), the samples were vortexed for 30 s, sonicated for 10 min in an ice-water bath, and incubated for 1 h at −40 °C to precipitate proteins. Then, the sample was centrifuged at 12,000 rpm (relative centrifugal force [RCF] = 13,800 (×g), R = 8.6 cm) for 15 min at 4 °C. The resulting supernatant was transferred to a fresh glass vial for analysis. The quality control (QC) sample was prepared by mixing an equal aliquot of the supernatants from all of the samples and was used to monitor the repeatability, stability, and reliability of the analytical methods.

LC–MS analysis

LC–MS/MS analyses were performed using an ultra-high performance liquid chromatography (UHPLC) system (Vanquish, Thermo Fisher Scientific) with a UPLC BEH Amide column (2.1 mm × 100 mm, 1.7 µm) coupled to a Q Exactive HF-X mass spectrometer (Orbitrap MS, Thermo Fisher Scientific) (17). The mobile phase consisted of 25 mmol/L ammonium acetate and 25 mmol/L ammonia hydroxide in water (pH = 9.75) (A) and acetonitrile (B). The autosampler temperature was 4 °C, and the injection volume was 2 µL. The Q Exactive HF-X mass spectrometer was used for its ability to acquire MS/MS spectra in information-dependent acquisition mode in the control of the acquisition software (Xcalibur, Thermo Fisher Scientific). In this mode, the acquisition software continuously evaluates the full scan MS spectrum. The electron spray ionization source conditions were set as follows: sheath gas flow rate of 30 Arb, aux gas flow rate of 25 Arb, capillary temperature of 350 °C, full MS resolution of 120,000, MS/MS resolution of 7500, collision energy of 10/30/60 in NCE mode, and spray voltage of 3.6 kV (positive) or −3.2 kV (negative).

Data processing

The raw data were converted to the mzXML format using ProteoWizard and processed with an in-house program, which was developed using R and based on XCMS, for peak detection, extraction, alignment, and integration (18). Then, an in-house MS2 database (BiotreeDB) was applied for metabolite annotation. The cutoff for annotation was set at 0.3.

After the relative standard deviation of the detected peaks was denoised, the missing values were filled up by half of the minimum value, and the internal standard normalization method was employed in this data analysis. The final dataset containing the information of peak number, sample name, and normalized peak area was imported to the SIMCA16.0.2 software package (Sartorius Stedim Data Analytics AB, Umea, Sweden) for multivariate analysis. Data were scaled and logarithmic transformed to minimize the impact of both noise and high variance of the variables. After these transformations, principal component analysis (PCA), an unsupervised analysis that reduces the dimension of the data, was carried out to visualize the distribution and grouping of the samples. 95% confidence interval in the PCA score plot was used as the threshold to identify potential outliers in the dataset. To visualize group separation and find significantly changed metabolites, supervised orthogonal projections to latent structures-discrimination analysis (OPLS-DA) were applied. Then, a 7-fold cross-validation was performed to calculate the value of R² and Q². R² indicates how well the variation of a variable is explained and Q² means how well a variable could be predicted. To check the robustness and predictive ability of the OPLS-DA model, a 200-times permutation was further conducted. Afterward, the R² and Q² intercept values were obtained. Here, the intercept value of Q² represents the robustness of the model, the risk of overfitting, and the reliability of the model, which will be the smaller the better. Furthermore, the value of variable importance in the projection (VIP) of the first principal component in OPLS-DA analysis was obtained (19). It summarizes the contribution of each variable to the model (20). The metabolites with VIP > 1 and p < 0.05 (Student’s t-test) were considered significantly changed metabolites. The results of the differentially abundant metabolites are visualized in the form of a volcano plot and were analyzed by correlation analysis to further understand the mutual regulatory relationship between metabolites in the process of biological state change. A receiver operating characteristic (ROC) curve was plotted for each well-defined differentially abundant metabolite, and the area under the curve (AUC) was calculated to judge the diagnostic efficiency of the metabolite.

In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway database (http://www.genome.jp/kegg/) was used for pathway enrichment analysis (21,22), and the corresponding differential abundance scores were determined.

Statistical analysis

The baseline characteristics of the study subjects were statistically analyzed using SPSS 25.0, with a difference of p < 0.05 being statistically significant. Binary categorical variables are represented as numbers (percentages) and were analyzed with chi-square tests. Continuous variables normally distributed are represented as the mean ± standard deviation and were analyzed with Student’s t-tests. Continuous variables that did not obey the normal distribution are represented as median (lower quartile, upper quartile) and were analyzed with Wilcoxon Mann–Whitney U tests.

Results

Baseline characteristics

A total of 40 children participated in this study: 30 asthmatic patients and 10 healthy controls. A comparison of the baseline characteristics of asthmatic and healthy subjects is presented in Table 1. There were no differences in age, sex, height, weight, BMI, white blood cell counts, percentage and count of neutrophils and eosinophils, IgE levels, or allergens (cat hair, dog hair) between children with asthma and healthy controls. There were significant differences in allergen of dust mites and lung function indicators such as FEV₁, forced vital capacity (FVC), FEV₁/FVC, peak expiratory flow (PEF), and forced expiratory flow (FEF)_25–75 between the two groups.

Table 1. Baseline characteristics of the 40 children.

	Asthma (n = 30)	Controls (n = 10)	95% CI	p Value
Age (years)	6.40 (5.98, 8.20)	6.90 (5.75, 8.13)	(−1.00, 1.00)	0.814
Sex
Male (n, %)	15 (50.0)	6 (60.0)	(0.16, 2.85)	0.855
Female (n, %)	15 (50.0)	4 (40.0)
Height (cm)	120.30 (116.65, 130.85)	123.90 (115.25, 131.08)	(−7.00, 7.30)	0.975
Weight (kg)	21.85 (19.98, 24.40)	23.00 (19.40, 26.58)	(−3.40, 2.20)	0.574
BMI (kg/m^2)	14.85 ± 1.40	15.25 ± 1.24	(−1.4, 0.61)	0.432
Blood routine examinations
WBC count (×10^9/L)	6.85 (6.08, 8.20)	6.80 (5.25, 7.73)	(−0.70, 1.50)	0.522
N (%)	61.16 ± 4.46	61.59 ± 5.20	(−3.86, 3.01)	0.803
NC	4.44 (3.67, 4.76)	3.89 (3.25, 4.80)	(−0.42, 1.06)	0.333
E (%)	4.40 (3.60, 5.80)	4.45 (3.50, 5.28)	(−0.80, 1.60)	0.731
EC	0.31 (0.26, 0.40)	0.25 (0.20, 0.38)	(−0.03, 0.14)	0.142
IgE (IU/mL)	173.50 (83.00, 279.75)	87.25 (60.95, 175.75)	(−21.0, 153.00)	0.174
Allergen
Dust mites (n, %)	29 (96.7)	1 (10.0)	(14.79, 4607.52)	<0.001
Cat hair (n, %)	1 (3.3)	0 (0.00)	–	1.000
Dog hair (n, %)	2 (6.7)	0 (0.00)	–	1.000
Lung-function tests
FEV1 (%)	74.00 (68.00, 78.00)	87.50 (79.00, 98.50)	(−24.00, −9.00)	<0.001
FVC (%)	91.43 ± 9.32	102.40 ± 12.27	(−18.43, −3.50)	0.005
FEV1/FVC (%)	0.79 (0.76, 0.83)	0.87 (0.85, 0.89)	(−0.11, −0.04)	<0.001
PEF (%)	83.27 ± 8.45	101.80 ± 11.88	(−25.47, −11.60)	<0.001
FEF25 (%)	81.33 ± 7.25	100.60 ± 11.30	(−25.46, −13.07)	<0.001
FEF50 (%)	75.50 (64.00, 80.00)	85.50 (72.00, 93.25)	(−22.00, −4.00)	0.006
FEF75 (%)	60.27 ± 6.51	78.30 ± 11.90	(−26.73, −9.33)	0.001

Data are presented as numbers (percentages), mean ± standard deviation or median (lower quartile, upper quartile) as required.

CI: confidence interval; BMI: body mass index; WBC: white blood cell; N: neutrophil; NC: neutrophil count; E: eosinophil; EC: eosinophil count; IgE: immunoglobulin E; FVC: forced vital capacity; FEV₁: forced expiratory volume in 1 s; PEF: peak expiratory flow; FEF: forced expiratory flow.

Metabolic profiling

From five QC samples and 40 experimental samples in this study, 9389 peaks were detected (Dataset S1), and 472 metabolites were retained after filtering of deviation and missing values (Dataset S2), filling of missing values, and data normalization. The dimensionality of the data was reduced by PCA, which effectively highlighted the overall distribution trend of the metabolomics data and the degree of differences between the samples. It can be seen from the dispersion point plot obtained by PCA (Figure 1 and Figure S1) that the samples are basically in the 95% confidence interval. However, due to the high dimensionality and small sample size of the metabolomic data, the differences between the groups were affected by confounding factors other than the grouping variables, and the overall metabolic level difference between the asthma group and the healthy control group, as shown in the PCA scatter plot, was not significant. Subsequently, through OPLS-DA, the irrelevant variables in metabolites were filtered out to obtain more reliable information about the intergroup differences in metabolites and the degree of correlation of the experimental group. The scatter plot of the OPLS-DA model (Figure 2a) shows that the difference between the two groups of samples is large and that the repeatability is good. Through the replacement test, the corresponding OPLS-DA model was established 200 times to obtain the R² and Q² values of the random model, which were higher than their original values (R²_x = 0.416; R²_y = 0.694; Q² = 0.272, p < 0.05; Figure 2b), indicating that the OPLS-DA model has good predictive ability.

Figure 1. Score scatter plot of PCA model for group asthma vs. control. The x- (PC1) and y- (PC2) axes indicate the first and second principal components, respectively. The circle represents the 95% confidence interval.

Figure 2. Score scatter plot and permutation plot test of OPLS-DA model for group asthma vs. control. (a) Score scatter plot of OPLS-DA model for group asthma vs. control. The x- (t[1]P) and y- (t[1]O) axes indicate predictive and orthogonal directions, respectively. The circle represents the 95% confidence interval. (b) Permutation plot test of OPLS-DA model for group asthma vs. control. R² and Q² indicate interpretability and predictability, respectively.

Potential metabolite-biomarker identification

According to VIPs > 1 and p < 0.05 of the OPLS-DA models, three metabolites associated with asthma could be selected, including traumatic acid (fatty acids), dodecanedioic acid (fatty acids), and glucobrassicin (sugars) (Table 2). As shown in Figure 3, the levels of traumatic acid and dodecanedioic acid in the urine samples of asthma patients are lower than those of healthy controls, and the expression of corresponding enzyme-encoding genes may be strongly inhibited, while glucobrassicin is the opposite. To further understand whether the changes between metabolites are coordinated, differentially abundant metabolite correlation analysis was performed. According to heatmaps, there was a positive correlation between traumatic acid and dodecanedioic acid (Figure 4). ROC curves were constructed, and the AUC was calculated to further evaluate the diagnostic efficacy of the above metabolites. All three metabolites had a certain accuracy (AUC 0.7–0.8) (Figure 5).

Figure 3. Volcano plot for group asthma vs. control. The x-axis indicates the change in multiples of each substance in the group (taking the logarithm base 2), and the y-axis represents the p value of the Student’s t-test (taking the negative logarithm base 10). Significantly up-regulated metabolites are shown in red, significantly down-regulated metabolites are shown in blue, and non-significantly differentiated metabolites are shown in gray.

Figure 4. Heatmap of correlation analysis for group asthma vs. control. The x- and y-axes represent the differential metabolites in this group of contrasts, with red indicating a positive correlation and blue indicating a negative correlation.

Figure 5. ROC analysis for distinguishing asthma patients from healthy controls. (a) ROC analysis of dodecanedioic acid (AUC = 0.75). (b) ROC analysis of glucobrassicin (AUC = 0.74). (c) ROC analysis of traumatic acid (AUC = 0.74).

Table 2. Identification of differentially abundant metabolites in urine between asthmatic and healthy subjects.

ID	Peak	RT	MZ	Mean asthma	Mean control	VIP	p Value	q Value	Fold change	Log_fold change
Upward
641	Glucobrassicin	318.132	227.1288629	0.009093848	0.022281495	1.13	0.012	0.396	6.52	2.70
Downward
432	Dodecanedioic acid	286.505	229.1445467	0.005731273	0.016753875	1.14	0.032	0.396	0.34	−1.55
307	Traumatic acid	132.739	447.0566545	0.001480022	0.000227166	1.34	0.029	0.394	0.41	−1.29

RT: retention time; MZ: mass-to-charge ratio; VIP: variable importance in the projection.

Construction of metabolic pathways based on differential metabolites

Through the analysis of metabolic and regulatory pathways, a more comprehensive and systematic understanding of the changes in biological processes caused by changes in experimental conditions can be obtained. The differentially abundant metabolites were introduced into the KEGG pathway database for analysis, obtaining a metabolic pathway – α-linolenic acid metabolism-Homo sapiens (human) (hsa00592) – in which the level of traumatic acid decreased (Figure 6). According to the enrichment results of differentially abundant metabolites in the KEGG metabolic pathway, this metabolic pathway belongs to lipid metabolism (p < 0.05, Rich Factor = 0.023) (Figure 7a). Differential abundance scores showed that the overall expression of the pathway tended to be downregulated (Figure 7b).

Figure 6. Metabolic pathways. Blue represents down-regulated significantly differentiated metabolites.

Figure 7. KEGG enrichment and differential abundance score for group asthma vs. control. (a) KEGG Enrichment for group asthma vs. control. The x-axis represents the Rich Factor, and the y-axis represents the KEGG metabolic pathway name. The color indicates the size of the p value, and the smaller the p value, the redder the color. (b) Differential Abundance Score for group asthma vs. control. The x-axis represents the DA Score and the y-axis represents the KEGG metabolic pathway name. The DA Score reflects the overall change of all metabolites in the metabolic pathway, and the score −1 indicates a downward trend in the expression of all annotated differential metabolites in the pathway.

Discussion

Asthma has an early onset compared to other chronic diseases and has a long-term negative impact on patients’ quality of life (23). However, the current diagnosis of asthma is mainly based on patients’ clinical symptoms, pulmonary function tests, and drug treatment response, which cannot determine whether patients have the disease in an early stage (8). Therefore, the prevention, diagnosis, and treatment of early asthma can be successfully achieved through noninvasive detection of airway dysfunction and inflammation and early diagnosis by identifying biomarkers (24). Metabolomics can reveal the underlying mechanism of disease through the quantitative analysis of small molecules in organisms. Therefore, in this study, we analyzed and compared the urine metabolic profiles of asthmatic children and healthy controls based on the LC–MS metabolomics method to deepen the understanding of the pathophysiological pathways of childhood asthma. In this study, we observed differences in three metabolic substances and one metabolic pathway between children with asthma and healthy controls, which means that we may be able to use these differences to prevent and manage asthma disease in the future.

Some previous studies have identified some metabolites in the metabolic profile of asthma as possible biomarkers by using biological fluids from clinical samples combined with metabolomic approaches. Among them, several studies have reached similar conclusions. For example, differences in the gut microbiota result in lower expression levels of α-linolenic acid metabolism in patients with allergic rhinitis than in nonallergic rhinitis individuals (25). Rohmann et al. reported that analysis of serum markers in allergic patients based on Fourier transform ion cyclotron resonance mass spectrometry (FT-ICR-MS) metabolomics revealed enrichment of the α-linolenic acid metabolic pathway and the found metabolic marker OPC6-CoA, which overlapped with the metabolic pathway of traumatic acid in our study (26). Asthma, as a heterogeneous disease of multiple gene–environment interactions, is often associated with type 2 inflammation, and many allergen components can induce the onset of asthma (9,27). Therefore, research on allergies can support the close relationship between α-linolenic acid metabolism and asthma to a certain extent. In addition, Zhu et al. analyzed the rhinitis metabolomic analysis of infants hospitalized with bronchiolitis and concluded that infants with metabolic type B (high abundance of amino acids and low abundance of polyunsaturated fatty acids) have a significantly higher risk of asthma with an enriched amino acid and α-linolenic acid/linoleic acid metabolic pathway (13). Taken together, these observations confirm that dysregulation of α-linolenic acid metabolism plays a key role in the development of childhood asthma and may serve as a new therapeutic target for asthma.

α-Linolenic acid is an indispensable polyunsaturated fatty acid in the human body, and its role as an immunomodulator has been shown to be associated with asthma (28). α-Linolenic acid can be converted into eicosapentaenoic acid in the human body, thereby inhibiting arachidonic acid metabolism and blocking leukotriene production in vitro and inflammatory processes in vivo (29,30). Observational studies and randomized controlled trials have also shown that polyunsaturated fatty acids have a protective effect on the development of asthma: maternal intake of high levels of omega-3 polyunsaturated fatty acids during pregnancy reduces the risk of asthma and wheezing in offspring (14). Meanwhile, a higher plasma ratio of omega-3 polyunsaturated fatty acids in childhood and adolescence was associated with a 22% lower risk of asthma in early adulthood (31). According to our experimental results, the downregulated traumatic acid that is the product of α-linolenic acid metabolism, along with dodecanedioic acid, belongs to the category of fatty acids and shows a positive correlation (Figure 4). Therefore, it can be inferred that fatty acids, including α-linolenic acid, will be affected to some extent in children with asthma, which may provide new therapeutic target for children with asthma in the future.

Our study has several potential limitations. First, the sample size of this study is relatively small, and the results cannot be completely extrapolated to all patients with asthma, which needs to be verified by larger studies. Second, this study is a cross-sectional study that lacks time and cannot confirm the causal relationship between urine metabolic changes and asthma. A longitudinal study to assess the temporal dynamics of the identified metabolic changes in children with asthma will be needed. Third, the molecular/cellular biological mechanisms of asthma and differentially abundant metabolites need to be further elucidated. Fourth, there are other possible confounding factors in this study, such as race, lifestyle, and diet.

Conclusions

In conclusion, this study suggests that metabolic differences can be observed in the urine samples of asthmatic patients and healthy controls, and α-linolenic acid metabolic pathways may be involved in the pathogenesis of asthma, which will help further study metabolically targeted strategies for asthma prevention and treatment.

Author contributions

Data curation: Jianlong Gao and Genli Pan. Formal analysis: Jianlong Gao, Sainan Jiang, Yiting Zhou, and Siyi Jin. Investigation: Jianlong Gao and Genli Pan. Methodology: Yiting Zhou, Siyi Jin and Yuan Jiang. Supervision: Genli Pan and Yuan Jiang. Validation: Genli Pan and Yuan Jiang. Writing – original draft: Sainan Jiang and Yiting Zhou. Writing – review and editing: Sainan Jiang and Yuan Jiang.

Ethical approval

The study was conducted in accordance with the Declaration of Helsinki and approved by the Ethics Committee of Deqing People’s Hospital (protocol code LL2021-K21 on 21 June 2021).

Consent form

Informed consent was obtained from all subjects and guardians involved in the study. Written informed consent has been obtained from the patients’ guardians to publish this paper.

Acknowledgements

We thank all children enrolled in this study and all guardians for their support. Also, we thank all doctors and nurses who have participated in this study. We would like to express our gratitude to BioTree Biotechnology Co., Ltd., Shanghai, China for the suggestions regarding the metabonomic analysis.

Disclosure statement

The authors have no conflicts of interest to declare.

Data availability statement

Data available on request due to restrictions, e.g., privacy or ethical.

Additional information

Funding

This research was funded by the National Natural Science Foundation of China [Grant No. 81801552] and the Huzhou Science and Technology Plan Project [Grant No. 2021GY54]. The funder had no role in the design, data collection, data analysis, and reporting of this study.