Composition and dynamics of intestinal fungi during the postnatal 2 months of very low birth weight infants

TL;DR

This study explores how the gut fungi of very low birth weight infants change during their first two months of life.

Contribution

The study provides new insights into the composition and dynamics of intestinal fungi in very low birth weight infants using modern sequencing techniques.

Findings

The detectable rate of fungi in very low birth weight infants peaked on day 3 and showed fluctuations over time.

VLBWI had distinct fungal diversity patterns compared to controls, with Ascomycota and Basidiomycota being dominant.

The functional metabolic activity of fungi in VLBWI was lower than in the control group.

Abstract

It has been found that intestinal fungi play a role in the composition of the intestinal microecology and in the formation and development of the immunity during childhood. We investigated the gut fungi composition of preterm infants to analysis composition and dynamics of intestinal fungi during the postnatal 2 months of very low birth weight infants. We collected feces from 34 very low birth weight infants (VLBWI) and 28 preterm infants with birth weight >1500 g. We extracted total fungal DNA from feces and analyzed the composition of gut fungus through ITS sequencing. The fungal detectable rate in the experimental group peaked on day 3 (85.19%), then gradually decreased and started to show an increasing trend again by day 28. There were significant differences in the alpha diversity of intestinal fungus between VLBWI and controls, and the VLBWI had its own characteristics at different time points in richness and diversity. A total of 10 phylums and 342 genera were identified in all VLBWI samples. The dominant fungal phylum of the VLBWI group is Ascomycota (50.3%)and Basidiomycota (48.8%). The functional metabolic activity of the experimental group was lower than that of the control group.

Conclusion: The composition and abundance of VLBWI intestinal fungal showed several alterations during the first 2 months of life. The prediction of gut microbiota function suggests that intestinal metabolic function may be altered in VLBWI.
What is Known:• A limited number of studies has been found that symbiont fungi may be able to calibrate host immunological responses, promote development of peripheral lymphoid organs, promote T cell responses, and even may be associated with the development of certain diseases, such as inflammatory bowel disease (IBD), NEC, and allergic diseases. However, previous studies on intestinal microecology have mainly focused on adults while neglecting the role of fungi in the gut of children due to the much lower abundance of intestinal fungi than bacteria, limitations of techniques for detecting fungi, the difficulty of obtaining samples, and the absence of largescale reference databases.What is New:• In recent years, the discovery and development of fungal detection technologies such as 18s rDNA sequencing technology, Internal Transcribed Spacer(ITS), and DNA fingerprinting technology have further broadened the perspective on the impact of intestinal fungal exposure in early life.

What is Known:

• A limited number of studies has been found that symbiont fungi may be able to calibrate host immunological responses, promote development of peripheral lymphoid organs, promote T cell responses, and even may be associated with the development of certain diseases, such as inflammatory bowel disease (IBD), NEC, and allergic diseases. However, previous studies on intestinal microecology have mainly focused on adults while neglecting the role of fungi in the gut of children due to the much lower abundance of intestinal fungi than bacteria, limitations of techniques for detecting fungi, the difficulty of obtaining samples, and the absence of largescale reference databases.

What is New:

• In recent years, the discovery and development of fungal detection technologies such as 18s rDNA sequencing technology, Internal Transcribed Spacer(ITS), and DNA fingerprinting technology have further broadened the perspective on the impact of intestinal fungal exposure in early life.

Figures

Specaccum species accumulation curves: sample size is on the abscissa, the number of observed species (ASVs) is on the ordinate, and the shading reflects the confidence intervals of the curves. The results reflect the rate of addition of new species observed when the sample size is continuously expanded over the course of sampling the population of the sample

Venn diagram of species differences between the two groups

Folding line graph of fungal detectable rate in each group

Two-dimensional sorting chart of samples for PCoA analysis. Each point in the figure represents a sample and points of different colors indicate different samples (groups). The closer the projection of two points on the axes, the more similar the community composition of these two samples in the corresponding dimension. The elliptical dashed circle refers to the 95% confidence ellipse (i.e., 95 out of 100 samples in this sample group will fall in it)

Histogram of LDA effect values for marker species. The vertical coordinates are the categorical units with significant differences between groups, and the horizontal coordinates visualize the logarithmic score of the LDA analysis for each categorical unit in a bar chart. The classification units are ranked according to the size of their scores, which describes their specificity in the sample grouping. Longer lengths indicate more significant differences in the classification units

Predictive function analysis of KEGG pathways for the fecal fungus

Full text

Introduction

Human beings have a rich and unique microecology in the channels that connect them to the outside world, such as the skin, digestive tract, respiratory tract, oral cavity, and genitourinary tract [1]. Among them, the intestinal microbiota is the most representative microbial community [1]. Previous studies [2, 3] on intestinal microecology have mainly focused on intestinal bacteria or selected adults as subjects while neglecting the role of fungi in the gut of children due to the much lower abundance of intestinal fungi than bacteria, limitations of techniques for detecting fungi, the difficulty of obtaining samples, and the absence of large-scale reference databases. In particular, people still know little about the dynamic influence of fungi on the formation of intestinal microecological populations in neonates.

VLBWI, defined as newborns with birth weight <1500 g, are a special category of newborns characterized by low immunity and susceptibility to infectious diseases such as late-onset sepsis (LOS) and necrotizing enterocolitis (NEC) [4]. It has been reported that the development of complications such as feeding intolerance (FI), NEC, and LOS in VLBWI may be related to their unique gut microflora composition [5, 6].

A limited number of studies has been found that symbiont fungi may be able to calibrate host immunological responses [7, 8], promote development of peripheral lymphoid organs [9], promote T cell responses [10], and even may be associated with the development of certain diseases, such as inflammatory bowel disease (IBD), NEC, and allergic diseases [11–15].

There still are few studies on the intestinal fungi of very low birth weight infants. Based on the above-mentioned, first of all, we collected stool samples from 62 preterm infants at different time points, including 34 VLBWI and 28 preterm infants with birth weight >1500 g, then analyzed their intestinal fungi by ITS sequencing. At last, we explored the characteristics of intestinal fungi in VLBWI.

Method

Trial participants

After obtaining informed consent from parents, preterm infants who were admitted to the NICU directly after birth with a birth weight of less than 1500 g were eligible to participate in the experimental group, and preterm infants who had a birth weight of greater than 1500 g were eligible to participate in the control group. Exclusion criteria were preterm infants with severe asphyxia at birth or had congenital malformations, genetic metabolic disorders, or the presence of shock or multi-organ failure. The trial was in line with the Declaration of Helsinki and approved by the Ethics Committee of Ningbo Women’s and Children’s Hospital, Zhejiang Province, China.

Sampling

Medical personnel collected 2 g of stool samples on day 1, day 3, day 7, day 14, day 21, day 28, and day 60 (experimental group only) after the admission of study subjects who met the criteria with a special specimen box for stool cleaning, and stored them immediately in a −80 °C refrigerator. Collection was stopped for study subjects who were discharged or died in the middle of the study.

DNA extraction and sequencing of the ITS genes

Microbic DNA was extracted by a QIAamp DNA Stool Mini Kit (Qiagen) according to the manufacturer’s instructions. Then, DNA was quantified by Nanodrop, and the quality of DNA extraction was assessed by 1.2% agarose gel electrophoresis. The fungal-specific gene fragment sequence was amplified by the universal fungal primers. The purification of polymerase chain reaction PCR products was performed using Vazyme VAHTSTM DNA Clean Beads, and Illumina bridge PCR-compatible primers were introduced. The amplified products were detected by Quant-iT PicoGreen dsDNA Assay Kit fluorescent reagent and Microplate reader (BioTek, FLx800) quantification instrument. Sequencing libraries were prepared using Illumina’s TruSeq Nano DNA LT Library Prep Kit, and finally, the resulting samples were subjected to ITS high-throughput sequencing.

The raw data obtained from the Illumina platform turned into effective Amplicon Sequence Variants (ASVs) after removing the primer sequence, trimming, merging, and filtering. The alpha diversity refers to the fungal diversity within each sample, and it was calculated by using Simpson’s reciprocal index, which describes how many ASVs prevail in each sample. The beta diversity expresses the difference between the samples in terms of the number and abundance of ASVs within an age group, and it was calculated with the Bray-Curtis dissimilarity index, predicting the metabolic function of a sample’s fungus, identifying differential pathways, and obtaining the species composition of specific pathways. Functional prediction of ITS gene sequences in the MetaCyc database using PICRUSt2 was followed by functional unit PCoA analysis, i.e., using Bray-Curtis distance matrix combined with principal coordinates analysis to expand sample functional differences in low dimensions. After obtaining the abundance data of metabolic pathways, we used Student’s t-test to try to identify metabolic pathways with significant differences between groups.

Statistical analyses

Student’s t-test and analysis of variance were performed for the measurement data, and chi-square test analysis was performed for the count data by SPSS 23.0. Graphpad Prism v.8.0.2. was used for plotting.

Result

Participants’ characteristics

Sixty-two preterm infants were divided into two groups according to birth weight, including 34 VLBWI and 28 preterm infants with birth weight > 1500 g. There were no differences between the two groups in terms of sex (t = 0.01, P = 0.921) and delivery mode (t = 0.01, P = 0.921). However, there were significant differences in terms of gestational age, preterm rupture of membranes, amniotic fluid contamination, mechanical ventilation, and central vein, and the positive results are shown in Table 1. Table 1. Participants’ characteristicsExperimental groupControl groupt/c2PPremature rupture of membraneYes914.38^a^0.036No2527ApneaYes3436.927^a^ < 0.001No3124Amniotic fluid contaminationYes2246.901^a^ < 0.001No3226Mechanical ventilationYes26618.627^a^ < 0.001No822Central venous catheterYes34933.266^a^ < 0.001No019Gestational age (weeks)29.9 ± 2.2233 ± 1.32− 6.833^b^ < 0.001Birth weight(g)1239.41 ± 172.641896.43 ± 289.43− 11.07^b^ < 0.001^a^χ^2^ calculations^b^Student’s t‐test

Data quality control and ASV cluster

A total of 15,222,063 original sequences were obtained. After removing low-quality sequences and denoising, 87.45% (13,311,092/15,222,063) of valid sequences were retained. There were 9,431,961 sequences in the experimental group with 84,213.93 sequences on average for each sample, and 3,938,448 sequences in the control group with 75,739.38 sequences on average, and the number of denoised sequences in the two groups was significantly different by Student’s t-test (t = 2.08, P = 0.039). Specaccum species accumulation curves (Fig. 1) were plotted for the total number of ASVs whose end of the curve was flattened out, indicating that the sample size was sufficient to reflect the species composition of the community.Fig. 1. Specaccum species accumulation curves: sample size is on the abscissa, the number of observed species (ASVs) is on the ordinate, and the shading reflects the confidence intervals of the curves. The results reflect the rate of addition of new species observed when the sample size is continuously expanded over the course of sampling the population of the sample

All valid sequences were clustered into ASVs at 100% similarity, and the ASV abundance was counted for each sample, a total of 2704 ASVs were detected in the two groups of samples. The unique ASVs accounted for 62.03% (1439/2320) of all ASVs in the experimental group and 21.42% (497/2320) in the control group. 16.55% (384/2320) ASVs were common to the two groups (Fig. 2). This result suggested that the two groups have their own characteristics in the ASVs.Fig. 2. Venn diagram of species differences between the two groups

Detectable rate

Successful annotation to fungal taxa was obtained from 164 samples from 319 samples of participants. The positive rate of all samples is 51.41% (Table 2). The fungal detectable rate was 61.88% (112/181) for the experimental group, which was significantly higher than 37.68% (52/138) for the control group (t = 13.599, P < 0.001) (Table 3). Among the time points, there was no statistically significant difference in the detection rates between the experimental group and the control group at day 14 (t = 1.175, P = 0.278), day 21 (t = 0.383, P = 0.536), and day 28 (t = 0.830, P = 0.362), whereas the detection rates of fungi in the experimental group at days 1, 3, and 7 were significantly higher than those in the control group and the specific positive results are shown in Table 3. The fungal detectable rate in the experimental group peaked on day 3 (85.19%), then gradually decreased and started to show an increasing trend again by day 28, while the detection rate of fungi in the control group was relatively stable within the first week of life and then showed an overall upward trend. However, the overall detection rate maintained a stable trend over time (Fig. 3). Table 2. Fungal detectable rate of all samples in different timeTotal sampleDetectable sampleDetectable rate (%)Day 1613150.820Day 3523159.620Day 7573052.630Day 14532750.940Day 21482143.750Day 28331648.480Day 6015853.330Total31916451.410Table 3The counts of samples with detected and non-detected fungal DNA for two groupsExperimental groupControl groupt**PNegative (n)Positive (n)Positive rate (%)Negative (n)Positive (n)Positive rate (%)D1122264.71018933.3305.9260.015D342385.1901783215.251 < 0.001D772275.86020828.57012.779 < 0.001Summary6911261.880765237.68013.599 < 0.001Fig. 3Folding line graph of fungal detectable rate in each group

Analysis of species composition

Alpha diversity

The Alpha diversity indices used in our study include Chao1, Simpson, Shannon, Pielou’s evenness, Observed species, and Goods’ coverage. Richness was characterized by Chao1 and Observed species; diversity was characterized by Shannon and Simpson; evenness was characterized by Pielou’s evenness; depth was characterized by Good’s coverage. By Student’s t-test of the diversity indices of the two groups, Chao1 and Observed species were significantly higher in the experimental group than in the control group (Table 4), suggesting that the species richness of the experimental group was significantly higher than that of the control group (Fig. 4a–f). Analyzing the alpha diversity results by ANOVA in the VLBWI group found that there are significant differences in Chao1, Shannon, and Observed species at different time points (Table 5), suggesting that the VLBWI had its own characteristics at different time points in richness and diversity. Further comparisons revealed that there were no significant difference between days 1, 3, 7, 14, 21, and 28 within the VLBWI group, while the alpha diversity at day 60 was significantly higher than at other time points (Fig. 5a–c). Table 4. The alpha diversity in two groupsExperimental groupControl groupPChao126.580 ± 22.77724.857 ± 8.4180.029Simpson0.226 ± 0.1910.266 ± 0.1630.177Shannon0.824 ± 0.8770.813 ± 0.5390.252Pielou’s evenness0.148 ± 0.1320.151 ± 0.0960.416Observed species25.363 ± 21.39723.946 ± 8.1240.030Goods’ coverage0.600 ± 0.000050.600 ± 0.000030.645P* < 0.05Fig. 4a**–f The box plots comparing the abundance between the two groups. a Chao1; b Observed species; c Simpson; d Shannon; e Pielou’s evenness; f Goods’ coverageTable 5Distribution of alpha diversity at different points in time of VLBWIChao1SimpsonShannonPielou’evennessObserved’speciesGoods’coverageDay 126.924 ± 22.5270.239 ± 0.1970.872 ± 0.8900.155 ± 0.13625.45 ± 19.9170.59994 ± 0.000063Day 326.384 ± 17.4390.262 ± 0.1910.894 ± 0.7860.162 ± 0.12624.843 ± 14.9260.59994 ± 0.000066Day 728.088 ± 23.6850.259 ± 0.1710.935 ± 0.8240.172 ± 0.11926.818 ± 21.9400.59994 ± 0.00005Day 1421.212 ± 19.8130.141 ± 0.1860.507 ± 0.7900.096 ± 0.12620.621 ± 18.7970.59996 ± 0.00004Day 2118.654 ± 5.6160.135 ± 0.1750.444 ± 0.5590.087 ± 0.10718.1 ± 5.5520.59996 ± 0.00002Day 2820.361 ± 7.1250.179 ± 0.1510.557 ± 0.3910.112 ± 0.07819.527 ± 6.4700.59995 ± 0.00003Day 6052.848 ± 49.5080.349 ± 0.2401.676 ± 1.5820.253 ± 0.20750.463 ± 48.1780.59991 ± 0.00006F2.3991.9442.3722.0942.5011.142P0.0330.0800.0340.0600.027*0.343Fig. 5a–c The box plots comparing the abundance between different time points in VLBWI. a Chao1; b Simpson; c Observed species

Beta diversity

PcoA analysis (Fig. 6) was performed, suggesting a trend of clustering of intestinal fungi between the two groups. Significance assessment of two groups using PERMANOVA revealed a significant difference in the distribution of intestinal fungal composition between the two groups (t = 1.85256, P = 0.008).Fig. 6. Two-dimensional sorting chart of samples for PCoA analysis. Each point in the figure represents a sample and points of different colors indicate different samples (groups). The closer the projection of two points on the axes, the more similar the community composition of these two samples in the corresponding dimension. The elliptical dashed circle refers to the 95% confidence ellipse (i.e., 95 out of 100 samples in this sample group will fall in it)

Overview of the gut microbiota composition and detective rate in VLBWI

A total of 10 phylums and 342 genera were identified in all VLBWI samples. The ten fungi that were detected were Ascomycota, Basidiomycota, unclassified_Fungi, unidentified, Mucoromycota, Mortierellomycota, Chytridiomycota, Basidiobolomycota, Glomeromycota, and Blastocladiomycota. The top two fungi in phylums, Ascomycota (50.3%), Basidiomycota (48.8%), comprised 99.1% of the fungal population abundance, who showed a zigzag pattern of change in the first 2 weeks. Ascomycota abundance was higher on days 1 and 7, lower on days 3 and 14, and then entered a more stable period, with an overall decreasing trend after day 21. However, the trend of Basidiomycota was just opposite to that of Ascomycota (Fig. 7a).Fig. 7a, b Folding line graph of fungal abundance in VLBWI at different time points a in the level of phylum; b in the level of genus

At the genus level, the top 50 genera in abundance accounted for 97.38% of the total fungal genera detected. The top ten dominant genera in abundance were Malassezia (21.2%), Candida (10.8%), Leucosporidium (8.45%), Cystofilobasidium (8.11%), Pseudogymnoascus (6.17%), Mycosphaerella (3.92%), Alternaria (3.74%), Cutaneotrichosporon (3.60%), Cladosporium (3.54%), and Aspergillus (3.37%) (Fig. 7b). Their combined abundance accounted for 72.9% of the total genera. The abundance of Malassezia and Candida, as the two dominant fungi, showed similar trends, with an overall increasing trend within 2 weeks after birth, then peaked in 2 weeks and a leveling off from 2 weeks to 28 days before starting to decline. A similar trend was observed for Cladosporium but differed in that it started to show an increasing trend again at day 60. The abundance of Cystofilobasidium and Pseudogymnoascus increased to a high level in the early stage, but after 1 week of birth, it had a decreasing trend over time. Leucosporidium were present at relatively high abundance in samples retained on day 1 after birth but showed an overall decreasing trend over time. At day 28, most fungal genera began to show a decreasing trend or reached a trough, while Cutaneotrichosporon reached a peak. Alternaria, Aspergillus, and Mycosphaerella maintained low abundance levels with small fluctuations throughout the observation period (2 months after birth).

Fungal variance analysis

Species difference analysis was performed on the phyla and genera between the two groups, and a total of 4 phyla and 25 genera with significant differences were detected. Ranking according to the significance of the differences from least to most. At the phylum level, the abundance of Mortierellomycota, Sarocladium, and Holtermanniella was significantly higher and Leohumicola was significantly lower in the experimental group. At the genus level, the abundance of Mortierella, Humicola, Pseudopithomyces, Lophiotrema, Pestalotiopsis, Chaetomium, Filobasidium, Cladosporium, and Malassezia was significantly higher in the experimental group, and the significantly higher abundance of Umbelopsis, Coniosporium, Neosetophoma, Schizophyllum, Curvularia, Scytalidium, Debaryomyces, Cystofilobasidium, Mycocentrospora, Leucosporidium was significantly lower (Fig. 8).Fig. 8. Histogram of LDA effect values for marker species. The vertical coordinates are the categorical units with significant differences between groups, and the horizontal coordinates visualize the logarithmic score of the LDA analysis for each categorical unit in a bar chart. The classification units are ranked according to the size of their scores, which describes their specificity in the sample grouping. Longer lengths indicate more significant differences in the classification units

Fungal community functions differ according to gestational age

The ITS feature sequences were aligned with the reference sequences to construct the evolutionary tree. Using the Castor algorithm, the nearest sequence species of the feature sequences were inferred based on the gene family copy number corresponding to the reference sequence in the evolutionary tree, and thus their gene family copy numbers were obtained. The gene family copy number of each sample is calculated by combining the abundance of each sample. Finally, the gene families were “mapped” to the MinPath database to infer the presence of metabolic pathways, and the abundance of metabolic pathways in each sample was obtained, and the MetaCyc database was used for metabolic pathway prediction. We compared the two sets of results and found that there were 9 differences between two groups in functional abundances, which were in glycolysis III (P = 0.042), gluconeogenesis I (P = 0.013), palmitate biosynthesis I (P = 0.024), D-myo-inositol (1,4,5)-trisphosphate biosynthesis (P = 0.01), pyrimidine deoxyribonucleotides de novo biosynthesis I (P = 0.043), phospholipid remodeling (P = 0.035), stearate biosynthesis III (P = 0.044), phosphatidylglycerol biosynthesis I (P = 0.013), and phosphatidylglycerol biosynthesis II (P = 0.013) (Table 6). We found that the expression of all function pathways was reduced in the experimental group compared to the control group, except for gluconeogenesis I, D-myo-inositol (1,4,5)-trisphosphate biosynthesis (Fig. 9). This result suggests that the functional metabolic activity of the experimental group was lower than that of the control group. Table 6. Names of metabolic pathways and metabolitesCoding of pathwaysPathway metabolic functionANAGLYCOLYSIS-PWYGlycolysis III (from glucose)GLUCONEO-PWYGluconeogenesis IPWY-5994Palmitate biosynthesis I (animals and fungi)PWY-6351D-myo-inositol (1,4,5)-trisphosphate biosynthesisPWY-7184Pyrimidine deoxyribonucleotides de novo biosynthesis IPWY-7409Phospholipid remodeling (phosphatidylethanolamine, yeast)PWY3O-355Stearate biosynthesis III (fungi)PWY4FS-7Phosphatidylglycerol biosynthesis I (plastidic)PWY4FS-8Phosphatidylglycerol biosynthesis II (non-plastidic)Fig. 9. Predictive function analysis of KEGG pathways for the fecal fungus

Discussion

Previous studies have mainly investigated the composition and role of intestinal bacterial communities, but little is known about fungal communities, especially about the long-term effects of fungi on the formation of early intestinal microecological populations [16]. However, intestinal fungi are a non-negligible part of the intestinal microbiota [17]. It has been found that colonization of neonates by intestinal fungi may not only be associated with the development of invasive fungal disease, but may also be involved in the establishment and maturation of the human immune system [12, 18–20], which may even be associated with the development of certain future diseases.

In this study, 34 VLBWI and 28 preterm infants with a birth weight of more than 1500 g who were hospitalized in Ningbo Women’s and Children’s Hospital from January 2021 to April 2022 were enlisted. The two groups differed in terms of gestational age, birth weight, preterm rupture of membranes, asphyxia, amniotic fluid contamination, mechanical ventilation, and central veins. Disease severity tends to be negatively correlated with gestational age, so gestational age can be associated with many covariates and potential factors of disease associated with preterm birth [21], which lead to these differences in participant’s characteristics between the two groups. By collecting stool samples from these participants at 7 time points in the experimental group and 6 time points in the control group for sequencing, trends in their intestinal fungal composition and dynamics were described. The result of quality control indicated that the sample size of this study was sufficient for data analysis, and the sequencing depth met the requirements.

Willis et al. [2] were able to detect Microeukaryotes (mainly of fungal origin) in the first meconium sample after birth, suggesting that the presence of intestinal fungi may retrace to the fetal period. In the present study, the fungal detectable rate was 50.82% on the first day for all stool samples and up to 64.71% for the VLBWI, which corroborate the above findings.

Experimental results shown that the fungal detectable rate was significantly higher in the experimental group (χ^2^ = 13.59, P =  < 0.001). Microbial infestation in the uterine cavity may cause preterm labor [22], such as intrauterine inflammation. And as a fungal causative agent, Candida intrauterine inflammation caused by Candida is a rare but recognized cause of preterm birth [22]. Therefore, we hypothesized that intrauterine fungal infection was one of the reasons for the higher detection rate of fungal infection in VLBWI.

The total fungal detectable rate in the experimental group increased from the first day (64.71%) to the third day (85.19%), followed by a smoothly decreasing trend, while in the control group, in contrast, gradually increased from a lower postnatal detectable rate (33.33%) to 62.5% on day 28. Candida detection rate showed a fluctuating trend within 2 months after birth, which peaked on the third day and troughs on the 21st day. Since there are few known studies using high-throughput sequencing technology for the analysis of neonatal intestinal fungi, therefore the literature available for comparison is insufficient. However, in the study [23] by Nahid Kondori et al. on the fungal cultures of Swiss infants, it was found that the Candida detectable rate from 3 days of age maintains an increase to 18 months of age; this difference with our experimental results may be caused by the different detection method [3], resulting in bias in the final results. More studies of the same experimental method are needed to compare with this experiment.

A comparison of alpha diversity between the two groups revealed a greater abundance of fungus in the VLBWI group. We speculate that the fact that VLBWI face more exposure to the environment, healthcare professionals, and the occurrence of more invasive operations (mechanical ventilation, central venous line placement, etc.) [16, 24, 25] after birth may contribute to their having a more abundance of fungus. There were differences in alpha diversity among VLBWI groups at time point. At day 60, alpha diversity was significantly higher than ever before, suggesting a significant increase in fungal diversity in the stool at day 60. In the prospective study of parent-offspring fungal colonization by Kasper Schei et al. [26], it was found that fungal alpha diversity in infants reached a minimum in 10-day samples and thereafter start increasing steadily from birth to 2 years of age. It was also found that species diversity gradually increased with infant age and diet [25]. In addition, a study [27] found that fungi were also present on the surface of the NICU environment and in breast milk. All of the above suggest that the differences in alpha diversity at different times may be related to the environment, the age itself, and the dietary intake in which the VLBWI was born.

In VLBWI, the combined abundance of Ascomycota (50.3%) and Basidiomycota (48.8%) was 99.1%, suggesting that the majority of intestinal fungus in the samples belonged to the above two phyla. This is generally consistent with the results of previous studies [2, 25, 26, 28]. At the genus level, most of the fungi in the VLBWI group consisted of a few fungal genus, with the most dominant fungal genera being Malassezia (21.2%) and Candida (10.8%). In contrast to the second-highest abundance of Candida found in this study, previous studies found that Candida ranked first in abundance among the detected intestinal fungi [2, 23, 25, 26]. Malassezia and Candida were present at different samples at different points in time, so we consider they more likely to be intestinal colonizing fungus rather than transient population. A study [2] found that the main composition of the fungal population shifted from Malassezia to Candida could be observed at 1 to 5 months, but the abundance of Malassezia was consistently higher than Candida during the period observed (postnatal to 2 months of age) in our study. A longer follow-up may be needed to verify previous observations.

In our study, in addition to the two fungal genera mentioned above, other dominant fungal genera are Leucosporidium (8.45%), Cystofilobasidium (8.11%), Pseudogymnoascus (6.17%), Mycosphaerella (3.92%), Alternaria (3.74%), Cutaneotrichosporon (3.60%), Cladosporium (3.54%), and Aspergillus (3.37%). This is not exactly the same as the previous studies by James et al. [25]. Previous studies found that the composition and abundance of intestinal fungi varied greatly between individuals [13, 25, 29], and symbiotic fungal populations are more variable than bacterial populations [30]. Some of the fungal genera found in these studies on neonates were similar to adults [3, 31, 32]. In the Human Microbiome Project, ITS2 and 18S rRNA sequencing results indicated that Saccharomyces, Malassezia, Candida, Cyberlindnera, Penicillium, Cladosporium, Aspergillus, Debaryomyces, Pichia, Clavispora, and Galactomyces are the most common fungal genera in the human gut [29]. Based on the above observation, we speculate that the fungal community in the gut during the earliest stages of life has an influence on the composition of the gut fungi in adults.

The most basic condition for fungal colonization in humans is the ability to survive at a temperature of 37 °C. Some of the detected genera of fungi, such as Cladosporium, are not able to grow at human temperatures, and its ability to produce spores suggests that it may be caused by the inhalation of spores present in the environment by the host [25]. Many fungi are usually present in soil and air and bound to plants as pathogens or saprophytes and enter the digestive tract by ingestion (as foodborne contaminants) or inhalation [25]. In our study, some fungi such as Sonoraphlyctis, Anthracocystis, and Thecaphora occur only a few times in very low abundance, so we consider that they are of environmental origin and are not able to colonize the intestine.

The VLBWI group showed enhanced D-myo-inositol (1,4,5)-trisphosphate biosynthesis. This pathway mediates the biological response of a large number of hormones and neurotransmitters in target cells by regulating calcium release from intracellular stores and make roles in controlling calcium homeostasis, transferring calcium between intracellular stores, and regulating calcium entry across the plasma membrane [33], while the expression of all other metabolic functions was diminished. Rozlyn et al. [13] reported that at 1 year of age, fungal communities in the infant gut demonstrate an increased capacity for functions related to energy metabolism and a decrease in degradation pathways relative to communities at 3 months of age, but no more experiments were performed to verify this finding.

In conclusion, this study describes the composition and dynamics of intestinal fungi over time during the first 2 months of life in VLBWI and provides information for future studies on the intestinal fungus of children. However, there are some limitations in this study: firstly, we can not determine whether the detected fungus is intestinal colonization or temporary foreign contamination, which requires more study subjects and a longer observation period to verify. Since fungi are ubiquitous in the environment, and fungal DNA represents a relatively low percentage of fecal DNA content (especially in meconium), rigorous extraction, purification, and amplification techniques are required [2], and it is difficult to differentiate environmental fungi from colonizing fungi from the sample. A prominent reason why the study of fungi lagging behind that of the microbiome is the lack of standardization of fungal bioassay methods [29]. Secondly, we did not include full-term newborns as a separate control group, so we cannot compare VLBWI with healthy infant for further analysis. Finally, the observation period of this study was only 2 months after the birth of the study subjects which is too short to identify and summarize characteristics. In the future, we will select full-term infants as a control group and extend the observation period to explore the correlation and influence of the establishment, change, and maturation of intestinal fungal community with some clinical factors, such as gestational age, sex, birth weight, invasive manipulation, and food intake. We also tried to further analyze the relationship between intestinal fungi and disease occurrence in the context of clinical diseases.