Associations between infant fungal and bacterial dysbiosis and childhood atopic wheeze in a nonindustrialized setting

Background Asthma is the most prevalent chronic disease of childhood. Recently, we identified a critical window early in the life of both mice and Canadian infants during which gut microbial changes (dysbiosis) affect asthma development. Given geographic differences in human gut microbiota worldwide, we studied the effects of gut microbial dysbiosis on atopic wheeze in a population living in a distinct developing world environment. Objective We sought to determine whether microbial alterations in early infancy are associated with the development of atopic wheeze in a nonindustrialized setting. Methods We conducted a case-control study nested within a birth cohort from rural Ecuador in which we identified 27 children with atopic wheeze and 70 healthy control subjects at 5 years of age. We analyzed bacterial and eukaryotic gut microbiota in stool samples collected at 3 months of age using 16S and 18S sequencing. Bacterial metagenomes were predicted from 16S rRNA data by using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States and categorized by function with Kyoto Encyclopedia of Genes and Genomes ontology. Concentrations of fecal short-chain fatty acids were determined by using gas chromatography. Results As previously observed in Canadian infants, microbial dysbiosis at 3 months of age was associated with later development of atopic wheeze. However, the dysbiosis in Ecuadorian babies involved different bacterial taxa, was more pronounced, and also involved several fungal taxa. Predicted metagenomic analysis emphasized significant dysbiosis-associated differences in genes involved in carbohydrate and taurine metabolism. Levels of the fecal short-chain fatty acids acetate and caproate were reduced and increased, respectively, in the 3-month stool samples of children who went on to have atopic wheeze. Conclusions Our findings support the importance of fungal and bacterial microbiota during the first 100 days of life on the development of atopic wheeze and provide additional support for considering modulation of the gut microbiome as a primary asthma prevention strategy.

Background: Asthma is the most prevalent chronic disease of childhood. Recently, we identified a critical window early in the life of both mice and Canadian infants during which gut microbial changes (dysbiosis) affect asthma development. Given geographic differences in human gut microbiota worldwide, we studied the effects of gut microbial dysbiosis on atopic wheeze in a population living in a distinct developing world environment. Objective: We sought to determine whether microbial alterations in early infancy are associated with the development of atopic wheeze in a nonindustrialized setting. Methods: We conducted a case-control study nested within a birth cohort from rural Ecuador in which we identified 27 children with atopic wheeze and 70 healthy control subjects at 5 years of age. We analyzed bacterial and eukaryotic gut microbiota in stool samples collected at 3 months of age using 16S and 18S sequencing. Bacterial metagenomes were predicted from 16S rRNA data by using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States and categorized by function with Kyoto Encyclopedia of Genes and Genomes ontology. Concentrations of fecal short-chain fatty acids were determined by using gas chromatography. Results: As previously observed in Canadian infants, microbial dysbiosis at 3 months of age was associated with later development of atopic wheeze. However, the dysbiosis in Ecuadorian babies involved different bacterial taxa, was more pronounced, and also involved several fungal taxa. Predicted metagenomic analysis emphasized significant dysbiosisassociated differences in genes involved in carbohydrate and taurine metabolism. Levels of the fecal short-chain fatty acids acetate and caproate were reduced and increased, respectively, in the 3-month stool samples of children who went on to have atopic wheeze. Conclusions: Our findings support the importance of fungal and bacterial microbiota during the first 100 days of life on the development of atopic wheeze and provide additional support for considering modulation of the gut microbiome as a primary asthma prevention strategy. (J Allergy Clin Immunol increased dramatically in incidence in industrialized countries since the 1960s, making it the most common chronic disease of childhood. 1,2 The increase in prevalence has been attributed primarily to environmental alterations rather than changes in genetic susceptibility. 2 Epidemiologic studies have revealed several potentially protective environmental factors, such as growing up on a farm, vaginal birth, breast-feeding, the presence of household pets, birth order, and number of children, as well an increased risk of asthma being associated with antibiotic use during late pregnancy and the first year of life. 3 The nature of these exposures has led many to hypothesize that microbial alterations (dysbiosis) driven by these factors and leading to modulation of the infant immune response are causally related to the increase in asthma and other atopic diseases.
Previously, we provided strong support for this hypothesis in mice [4][5][6] and human subjects, 7,8 in which we found a critical window early in life during which changes in the gut microbiome are influential in immune dysregulation and can lead to atopy and asthma. 9 Importantly, we identified specific bacterial taxa associated with reduced asthma risk in a prospective longitudinal cohort study of Canadian babies (Canadian Healthy Infant Longitudinal Development [CHILD] study). 10 These microbial taxa were not only associated with disease but also implicated causally in ameliorating experimental asthma in a murine model. 7 A more recent microbiome study in a US pediatric cohort also showed bacterial dysbiosis in 1-month-old infants at high risk of asthma. 11 Interestingly, this study also measured the fecal fungal community and found several fungal taxa also involved in the gut dysbiosis associated with asthma. 11 Traditionally not included in microbiome studies, fungi are important members of the gut microbiome, especially in the first months of life, when fungi are present at a much higher diversity than later in life. 11 Because asthma rates continue to climb in affluent countries, the less developed regions of the world are now also experiencing increases in asthma prevalence, particularly urban and urbanizing populations. The International Study of Asthma and Allergies in Childhood 12,13 reported asthma prevalence worldwide and showed that some of the highest rates were in Latin American countries, such as Peru, Brazil, and Costa Rica. Although this has been taken as evidence against the hygiene hypothesis, rural populations within these countries have lower rates of asthma than urban populations, and rural-to-urban migration and modernization can contribute to increases in asthma prevalence. 14 Thus changes in microbial exposure driven by urbanization might help explain regional differences in asthma within Latin America.
The intestinal microbiome differs by geographic location, which is likely explained by differences in diet, lifestyle, and environmental exposures. 15 Not surprisingly, the limited data available on the microbiome of Latin Americans show notable differences between the gut microbiome from this part of the world and North American or European populations. 15,16 Given what we have reported previously in Canadian infants, we explored whether similar patterns can be observed in a geographically distinct population with similar reported rates of asthma prevalence to Canada (approximately 10% of the population). 17 To do this, we completed a case-control study nested within a birth cohort from a rural district in a tropical coastal region of Ecuador, in which we selected children who had atopic wheeze at 5 years of age and healthy control subjects. We compared the bacterial and fungal microbiomes and levels of short-chain fatty acids (SCFAs) as a measure of bacterial carbohydrate metabolism in fecal samples collected at 3 months and contrasted our findings with those obtained from a Canadian cohort. Here we show that bacterial and fungal dysbiosis can be detected in early life in babies who go onto have atopic wheeze in later childhood in a nonindustrialized setting.

Study design
We conducted a case-control study nested within the ECUAVIDA birth cohort study designed to study prospectively the effects of early-life parasite infections on the development of atopy and allergic diseases. The cohort of 2404 newborns was recruited between 2006 and 2009 in the public hospital serving the rural district of Quinind e, Esmeraldas Province, a tropical region of coastal Ecuador. Detailed information for the cohort is provided elsewhere. 18 Cases and control subjects were selected from children with a stool sample collected at 3 months of age. As outlined in more detail below, cases were defined as children with maternally reported wheeze in the previous 12 months at 5 years of age and with evidence of atopy based on a positive skin prick test response. Healthy control subjects were a random sample of children with no previous history of wheeze and no evidence of atopy at 5 years of age selected by using the sample option in Stata 11 software (StataCorp, College Station, Tex). The study protocols were approved by the Ethics Committees of the Hospital Pedro Vicente Maldonado and Universidad San Francisco, Quito, Ecuador. The study is registered as an observational study (ISRCTN41239086). Informed written consent was obtained from the child's mother.

Data collection and study procedures
Data on wheezing were collected by using maternal questionnaires administered by a trained physician (M.E.C. and M.V.). Wheeze was defined as any episode of reported wheeze during the previous 12 months. Other data were collected by using periodic questionnaires administered during the first 2 weeks of life and then at 7 and 13 months. Blood was collected at 7 and 13 months and analyzed for eosinophil counts by using standard methods. Stool samples were collected from the child at 3 months and analyzed by using microscopic methods, as previously described. 19 An aliquot was stored at 2808C until analysis.
Atopy was measured based on skin prick test reactivity at 5 years of age using 9 common allergen extracts (Greer Laboratories, Lenoir, NC): Dermatophagoides pteronyssinus/farinae mix, American cockroach (Periplaneta americana), cat, dog, grass pollen (9 southern grass mix), fungi, egg, milk, and peanut, with positive histamine and negative saline controls. A positive reaction was defined as a mean wheal diameter of at least 3 mm greater than that elicited by the saline control, and a positive skin prick test response was defined as a positive reaction to any of the allergens tested. No children were taking antihistamines at the time of skin testing.

Microbial community analysis
Extraction and amplicon preparation. All samples were amplified by means of PCR in triplicate by using barcoded primer pairs flanking the V4 region of the 16S gene with the forward primer 59-GTGCCAGCMGCCGCGGTAA-39 and the reverse primer 5-GGACTACHVGGGTWTCTAAT-39, as previously described. 20 16S. Each 50 mL of PCR contained 22 mL of water, 25 mL of TopTaq Master Mix, 0.5 mL of each forward and reverse barcoded primer, and 2 mL of template DNA. The PCR program consisted of an initial DNA denaturation step at 958C for 5 minutes, 25 cycles of DNA denaturation at 958C for 1 minute, an annealing step at 508C for 1 minute, an elongation step at 728C for 1 minute, and a final elongation step at 728C for 7 minutes. Controls without template DNA were included to ensure that no contamination occurred. Amplicons were run on a 2% agarose gel to ensure adequate amplification. Amplicons displaying bands at approximately 250 bp were purified with the Illustra GFX PCR DNA Purification Kit (GE Healthcare, Fairfield, Conn). Purified samples were diluted 1:50 and quantified with PicoGreen (Invitrogen, Carlsbad, Calif) in the TECAN M200 (excitation at 480 nm and emission at 520 nm).
18S. DNA was sent to Integrated Microbiome Resource at Dalhousie University (Halifax, Nova Scotia, Canada) for amplification and sequencing. The 18S gene was amplified with the primers E572F (59-CYGCGGTAATTC CAGCTC-39) and E1009R (59-AYGGTATCTRATCRTCTTYG-39), and the reaction included a PNA blocking primer to reduce amplification of mammalian sequences (59-TCTTAATCATGGCCTCAGTT-39). Amplification was done in duplicate, with one reaction using undiluted DNA and one reaction using DNA diluted 1:10 in PCR water. Amplification was done according to previously described protocols. 21 Sequencing. 16S. Each PCR pool was analyzed on the Agilent Bioanalyzer (Agilent Technologies, Santa Clara, Calif) by using the High Sensitivity DS DNA assay to determine approximate library fragment size and to verify library integrity. Pooled library concentrations are determined by using the KAPA Library Quantification Kit (Illumina, San Diego, Calif). Library pools are diluted to 4 nmol/L and denatured into single strands by using fresh 0.2 N NaOH. The final library loading concentration is 8 pmol/L, with an additional PhiX spike-in of 20%. The amplicon library is sequenced on the MiSeq by using the 8MiSeq 500 Cycle V2 Reagent Kit (250 3 2).
18S. PCR products were visualized on E-Gels (Thermo Fisher Scientific, Waltham, Mass). The PCR product was quantified with the Invitrogen Qubit with PicoGreen and pooled at equal concentrations, according to the method of Comeau et al. 21 PhiX was spiked in at 5%, and the resulting library was sequenced at Dalhousie University on the Illumina MiSeq with the 8MiSeq 500 Cycle V2 Reagent Kit (250 3 2).
Bioinformatics. Sequences were preprocessed, denoised, and quality filtered by size by using Mothur MiSeq SOP (16S) 22 or QIIME (18S). 23 Demultiplexed reads were trimmed to a uniform length of 250 bp by using the FastX Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) and clustered into operational taxonomic units (OTUs) by using the Minimum Entropy Decomposition (MED) method, 24 as implemented in the Oligotyping microbial analysis software package. 25 MED is used to perform de novo taxonomic clustering with Shannon entropy to separate biologically meaningful patterns of nucleotide diversity from sequencing noise; processed data are partitioned into phylogenetically homogeneous units (MED nodes) for downstream bacterial diversity analyses. This analysis was carried out with the minimum substantive abundance parameter (-M) set at 250 reads. All other parameters were run with default settings; the maximum variation allowed per node (-V) was automatically set at 3 nucleotides.
Representative sequences were classified by clustering against the Greengenes database at 97% similarity (16S) 26 or SILVA123 at 99% similarity (18S). 27 The 16S data set was filtered to remove mitochondria and chloroplast sequences. The 18S data set was filtered to remove mammalian sequences, plants, and all OTUs present in fewer than 3 samples. After filtering, a cutoff of 1000 reads per sample was applied. All 16S samples passed the cutoff. In the 18S data set many samples had fewer than 1000 sequences per sample after filtering. Thus 10 samples from atopic patients with wheeze (AWs) and 22 control samples were retained for statistical analysis of the 18S data set. It is important to apply this sequencing cutoff because a low number of sequences might be a consequence of amplification or sequencing artifacts that might not reflect a real microbial community and will not provide enough statistical power for analysis. Singleton OTUs and any OTUs present less than 3 times among all samples were removed from both the 16S and 18S analyses. We constructed a phylogenetic tree for Pichia species by first extracting sequences from the aligned 99% SILVA 123 data set that corresponded to the genus Pichia or to the family Pichiaceae with RAxML. 28 The short Illumina reads were then aligned to the SILVA 123 data set by using PyNAST 29 and placed in the reference Pichia species tree by using the maximum likelihood RAxML evolutionary placement algorithm. 30

Quantitative PCR
Quantitative PCR (qPCR) was performed on genomic DNA isolated from stool samples by using the QIAamp Fast DNA Stool Mini Kit (Qiagen) to validate the high-throughput sequencing results. After extraction, all DNA concentrations were measured with the Quant-iT PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific), and template DNA used in the qPCR reactions was obtained from stocks normalized to 1 or 0.2 ng/mL (samples 2, 7, 11, 12, 14, 28, 54, 56, 59, 61, 66, 67, 71, 73, 79, and 93) and diluted as necessary.
Pichia kudriavzevii-specific primers targeting the 18S rRNA gene of this organism (forward: 59-CAT TCC GGG GGG CAT GCC T-39 and reverse: 59-GGC CAG CTT CGC TCC CTT TCA-39) and fungi-specific 18S primers (FR1: 59-AIC CAT TCA ATC GGT AIT-39 and FF390: 59-CGA TAA CGA ACG AGA CCT-39) 31 were used to measure the absolute abundance of P kudriavzevii and fungal 18S gene copies, respectively, by using the absolute quantification method. Primers for P kudriavzevii were designed by using CLUSTALW and National Center for Biotechnology Information Primer BLAST software, and their specificity for P kudriavzevii was verified by sequencing the amplicons obtained from a PCR reaction done with genomic DNA combined from 5 chosen subjects shown to be enriched in P kudriavzevii based on the sequencing data.
qPCR standard curve generation. 18S standard curves (see Fig E1 in this article's Online Repository at www.jacionline.org) were generated by using 1:5 dilutions of a 1 ng/mL stock of 18S amplicons obtained from PCR reactions done with the fungi-specific 18S primers (FR1 and FF390) and purified P kudriavzevii template DNA. Template DNA was extracted by using the YeaStar Genomic DNA Kit (Zymo Research, Irvine, Calif) from pure liquid cultures (ATCC 6258) grown for 24 hours at 308C in YPD media. A P kudriavzevii standard curve was generated by performing 1:5 dilutions of a 1 ng/mL stock of the same purified P kudriavzevii DNA (see Fig E1). qPCR reactions. All qPCR reactions were performed in duplicate on all stool samples for which there was sufficient DNA (n 5 95; samples 3 and 19 were not included) by using the 7500 Fast Real-Time System (Applied Biosystems, Foster City, Calif). Each 10-mL qPCR reaction done with the P kudriavzevii-specific primers contained 5 mL of iQ SYBR Green Supermix (Bio-Rad Laboratories, Hercules, Calif), 0.5 mL of each of the forward and reverse primers, and 4 mL of 0.2 ng/mL template DNA. The qPCR program consisted of 2 initial steps of 20 seconds at 558C and 10 minutes at 958C, followed by 40 cycles of 30 seconds at 958C and 60 seconds at 608C and a final cycle of 958C at 15 seconds, 608C at 1 minute, 958C at 15 seconds, and 608C at 15 seconds.
Each 20-mL reaction used for 18S rDNA quantitation contained 10 mL of iQ SYBR Green Supermix, 2.5 mL of each of the FR1 and FF390 primers, 3 mL of nuclease-free water, and 2 mL of 1 ng/mL genomic DNA. The qPCR program consisted of an initial step at 958C for 10 minutes, followed by 40 cycles of 15 seconds at 958C, 30 seconds at 508C, and 60 seconds at 708C, with a final cycle of 958C at 15 seconds, 608C at 1 minute, 958C at 15 seconds, and 608C at 15 seconds.

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States
We used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) 32 to generate a profile of putative functions (through metagenomic prediction) from the 16S rRNA OTU data. Predicted metagenomes from all samples were categorized by function at Kyoto Encyclopedia of Genes and Genome orthology level 3. To test for significant differences in functional category abundances between the AW and control samples, we used the Welch t test implementation of STAMP. 33 We also tested for differentially abundant predicted genes with DESeq2 34 under default settings. DESeq2 estimates variance and mean dependence in count data from high-throughput data sets and tests for a differential expression based on a model using negative binomial distribution. The test statistics' P values were adjusted for multiple testing by using the procedure of Benjamini and Hochberg (false discovery rate threshold, 5%).

SCFA analysis
Infants' fecal samples were combined with 25% phosphoric acid, vortexed, and centrifuged until a clear supernatant was obtained. Supernatants were submitted for gas chromatography to the Department of Agricultural, Food and Nutritional Science of the University of Alberta (Edmonton, Alberta, Canada). Samples were analyzed, as previously described, 35 with modifications. Briefly, samples were combined with 4-methyl-valeric acid as an internal standard, and 0.2 mL was injected into the Bruker SCION 456 gas chromatograph by using a Stabilwax-DA 30-m 3 0.53-mm 3 0.5-mm column (Restek, Bellefonte, Pa). A standard solution containing acetic acid, proprionic acid, isobutyric acid, butyric acid, isovaleric acid, valeric acid, and caproic acid, combined with internal standard, was injected in every run.
Programmable temperature vaporization injector and flame ionization detector temperatures were held at 2508C for the entire run. The oven was started at 808C and immediately ramped to 2108C at 458C/min, where it was held for 5.11 minutes. Total run time was 8.00 minutes. Helium was used at a constant flow of 20.00 mL/min. Sample concentrations were normalized to the wet weight of feces.

Statistical analysis
Differences in frequencies in categorical and continuous variables between cases and control subjects were evaluated by using x 2 and Mann-Whitney tests, respectively. The Fisher exact test was used for the reduced sample available for fungal sequencing analysis. We assessed fecal microbial diversity and relative abundance of bacterial and fungal taxa by using phyloseq, 36 along with additional R-based computational tools [37][38][39][40][41][42] in Rstudio software. Principal coordinates analysis was conducted by using phyloseq (Bray-Curtis as distance metric) and statistically confirmed with permutational multivariate ANOVA. The Chao1 index was calculated by using phyloseq and statistically confirmed with the Mann-Whitney test (GraphPad Prism software, version 5c; GraphPad Software, La Jolla, Calif). The R packages DESeq2 34 and multivariate association with linear models (MaAsLin) 43 were used to calculate differentially abundant OTUs between the control and AW samples. From the available metadata, the following covariates were forced into the MaAsLin model based on significant associations with the 5-year AW phenotype and on previously reported associations with asthma or with microbiome shifts: antibiotic use during pregnancy or the first year of life, duration of antibiotic use during pregnancy or the first year of life, type of delivery, household potable water, number of respiratory tract infections during the first year of life, eosinophilia at 7 months, and number of diarrheal episodes during the first year of life. All ecological measurements, except for the Chao1 index, were performed on the 100 most abundant OTUs after normalizing OTUs for relative abundance, which represent 98.14% and 99.83% of total OTUs in the 16S and 18S data sets, respectively.

Study subject characteristics: AWs and control subjects
From 1066 stool samples collected from infants in the ECUAVIDA cohort at 3 months of age with follow-up data at 5 years of age (2090), we identified 27 AWs and a random sample of 70 healthy control subjects with no history of wheeze and no atopy at 5 years. The required sample size for microbiome analysis was determined by using MetSizeR, a software package designed to estimate sample size for high-throughput experiments. Assuming at least 100 parameters measured and that 10% would differ significantly between cases and control subjects, a sample size of 24 was required to maintain a false discovery rate of less than 0.05. Characteristics of AWs and control subjects are shown in Table I. Most variables were distributed similarly between the 2 groups with the exceptions of mode of delivery (caesarean delivery: 37% for AWs vs 16% for control subjects, P 5 .022) and number of documented acute respiratory tract infections, which were more frequent in AWs compared with control subjects (P 5 .011). There were no significant differences with respect to baseline characteristics between AWs included in this analysis and those not included for lack of a 3-month stool sample (see Table E1 in this article's Online Repository at www.jacionline.org).

Bacterial and fungal alterations associated with AW phenotype
As expected, the gut microbiome between 3-month-olds in the Canadian CHILD Study and ECUAVIDA study differed substantially (see Fig E2 in this article's Online Repository at www.jacionline.org). In line with what we reported previously in Canadian infants, the AW phenotype did not explain any significant changes in a or b bacterial diversity (Mann-Whitney and permutational multivariate ANOVA, respectively; see Fig E3, A and B, in this article's Online Repository at www.jacionline.org). In addition to the bacterial profile, we also sequenced the V4 region of the 18S rRNA gene to obtain a profile of the intestinal eukaryotic microbiome in these infants. After subtracting plant and human sequences, our analysis found only sequences of fungal origin in the fecal microbiome of these 3-month-old babies. In this data set many samples had fewer than 1000 sequences per sample after applying filtering steps, which prompted their removal from the data set. Thus only 10 AW and 22 control samples were retained for statistical analysis of the 18S data set. However, AWs and control subjects used for fungal sequence analysis were reasonably representative of the original case and control samples (see Table E2 in this article's Online Repository at www.jacionline.org).
Similar to the 16S analysis, no significant differences were found in fungal a or b diversity between infants with atopic wheeze and control subjects (see Fig E3, C and D). Despite the lack of broad microbial community shifts between AWs and control subjects, several OTUs of bacteria (Fig 1) and especially fungi (Fig 2) were found to be differentially abundant by using DESeq2, a method of differential analysis of high-throughput sequencing data based on a model of negative binomial distribution. Among the observed bacterial alterations were 2 OTUs identified as Veillonella species, 3 OTUs identified a Streptococcus species, 1 OTU identified as Bacteroides species, 1 OTU identified as Ruminococcus gnavus, and 1 OTU identified as Bifidobacterium species (DESeq2; P < .05, Wald test and false discovery rate). Interestingly, 5 of these OTUs were among the top 10 most abundant bacterial OTUs in the data set. Notably, more than twice number of fungal OTUs (19) were significantly different between AWs and control subjects compared with bacterial OTUs, with 7 of them identified as the yeast P kudriavzevii (Fig 2, B). The larger differences in the fungal versus bacterial microbiome were also evident in relative abundance plots of the top 100 most abundant bacterial and fungal taxa, organized by genus (Figs 1, A, and 2, A). Additionally, we found several correlations between bacterial and fungal OTUs, suggesting potential ecologic synergistic and antagonistic associations occurring in the infant gut (see Fig E4 in this article's Online Repository at www.jacionline.org).
Given extensive metadata collected in this cohort, we captured the effect of the gut microbiome on the atopic wheeze phenotype (used as a surrogate for asthma risk) while controlling for potential confounding using MaAsLin. 43 MaAsLin is a multivariate linear modeling tool with boosting that tests for associations between specific microbial taxa and continuous and/or Boolean metadata. In microbiome studies certain variables, such as diet, age, and sample origin, can have large effects on microbiome composition. Thus it is important to subtract their influence from the potential effect caused by other metadata in the study, such as the phenotype. With MaAsLin, we found a significant association between most of the differential taxa obtained through DESeq2 analysis and atopic wheeze, including the bacterial genera Bifidobacterium, Streptococcus, and Veillonella, as well as the fungi P kudriavzevii, Saccharomyces, and Saccharomycetales (see Tables E3 and E4 in this article's Online Repository at www.jacionline.org), supporting the strong association between these microbial changes and the phenotype. By using both statistical approaches, it was evident that fungal alterations were more marked than bacterial dysbiosis.
Because of the strength of the association between the yeast P kudriavzevii and atopic wheeze risk, we validated the sequencing result using maximum likelihood to place individual sequences assigned as P kudriavzevii into a phylogenetic tree of all sequences assigned to Pichia species yeast in the reference database and found that all of the OTUs fell within the P kudriavzevii clade (see Fig E5 in this article's Online Repository at www.jacionline.org). We further validated this result by designing and optimizing specific primers for this yeast. Using qPCR on all AW and control samples, we confirmed the increased prevalence of Pichia species in the stool of 3-month-old babies who went on to have atopic wheeze in this Ecuadorian cohort (Fig 3, A). Although there is overlap in values between groups, almost a third of child AWs (8/27) had a higher abundance of P kudriavzevii than the highest abundance among control subjects (P 5 .001).
In addition, we found that AW samples contained a higher proportion of fungal reads over all eukaryotic reads (Fig 3, B). To confirm this bioinformatic observation, we quantified total fungal 18S rRNA gene copies using qPCR in all samples. Notably, we found a small yet significant increase in 18S rRNA gene copies in the AW samples (Fig 3, C), suggesting that fungi are overrepresented in the gut microbiome of babies who go on to have atopic wheeze.

Functional microbiome alterations in babies who had atopic wheeze
To determine the functional potential of the bacterial dysbiosis detected by using 16S sequencing, we predicted the microbiome Fungus-specific qPCR and Illumina sequence data. A, qPCR quantification (standard curve method) of P kudriavzevii in feces of 3-month-old infants from the ECUAVIDA study (control subjects, 68; AWs, 27). B, Relative abundance of fungal reads over the total number of eukaryotic reads from 18S data set. Only samples with more than 1000 sequences after filtering were included in the 18S high-throughput sequencing analysis (control subjects, 23; AWs, 10). C, qPCR quantification (standard curve method) of fungus-specific 18S copies in feces of 3-month-old infants from the ECUAVIDA study (control subjects, 68; AWs, 27). In Fig 3, A and C, graphs are based on a logarithmic scale. All graphs show Mann-Whitney tests: *P < .05 and ***P < .001. metagenome using PICRUSt. 32 PICRUSt is a computational method that predicts functional composition based on the 16S rRNA gene content and copy number and reference bacterial genomes. By comparing the predicted genetic composition between the 2 groups, we found taurine and hypotaurine metabolism, polyketide sugar unit biosynthesis, and carbohydrate digestion and absorption pathways were decreased in AWs (Welch t test, see Fig  E6 in this article's Online Repository at www.jacionline.org).
We investigated further the predicted change in bacterial carbohydrate metabolism by measuring SCFA concentrations in feces. Similar to what we found in babies from the CHILD study, there was a significant decrease in acetate concentrations in the stool of babies in the ECUAVIDA study who went on to have atopic wheeze (Fig 4, A), suggesting similar gut bacterial fermentation patterns between the 2 populations. In addition, we detected an increase in caproate levels associated with the AW phenotype in this cohort (Fig 4, F).

DISCUSSION
In the present study we explored the association between gut bacterial and fungal dysbiosis in infancy and the development of atopic wheeze in later childhood in a representative sample of an Ecuadorian birth cohort (ECUAVIDA), which has studied a population living in a rural district in a tropical nonindustrialized country. This represents a distinct environment from all previous studies of the effects of the microbiome on asthma development done in developed settings. In agreement with previous analyses from Canadian 7,8 and US 11 birth cohorts, our findings provide further confirmation of the presence of a microbial signature in the first 100 days of life associated with pediatric risk. A direct comparison between these new findings and the dysbiosis present in the Canadian CHILD study is not entirely possible because of differences in study design. In the CHILD study atopic wheeze was defined at 1 year of age, whereas this study measured atopic wheezing at 5 years as a surrogate for asthma risk. However, the observed differences in gut dysbiosis associated with the risk of atopic wheezing in both cohorts could also originate from differences in gut microbial composition, given that the microbiome of babies in ECUAVIDA differed substantially from that in age-matched Canadian babies in the CHILD study. As has been shown previously, 44 modernization processes (or westernization) strongly affect the human gut microbiome. Although this study was not designed to study the factors that lead to such differences, the humid warm weather in this region and less hygienic conditions are likely to increase exposures to microbes, including those that cause childhood infectious diseases, which might affect health and modify the infant microbiome. 14,45 The differences in the gut microbiome detected in this cohort yielded important information. First, the microbial signature associated with atopic wheeze in these infants was more pronounced than in the Canadian cohort. In addition, except for Veillonella species, the bacterial dysbiosis associated with atopic wheezing in this cohort involved different taxa to those detected in the Canadian study. Despite major compositional differences associated with atopic wheezing between the Canadian and Ecuadorian cohorts, we found important metabolic similarities between the 2 cohorts. Notably, a decrease in fecal acetate in atopic infants with wheeze was a common finding in the 2 populations, suggesting that alterations in fermentative patterns that led to a decrease in acetate production might be a common feature associated with atopic wheezing. A reduction in carbohydrate metabolism predicted by PICRUSt also supports the association between SCFA production and atopic wheeze. Acetate, 46 butyrate, 47 and propionate 46 have been shown to stimulate regulatory T-cell proliferation and attenuate allergy-associated T H 2 immune responses. 48 Given our data in 2 very different cohorts, fecal acetate levels reductions should be studied further in larger longitudinal cohorts as a possible biomarker of asthma risk in young babies.
Prompted by the results from a recent study in US infants, 11 we profiled the eukaryotic microbiome of the Ecuadorian babies using the 18S rRNA marker gene and a custom-annotated taxonomic database. Interestingly, after removing mammalian and plant-derived sequences, only fungal sequences remained in the data set. Given that we have successfully detected other eukaryotic commensals using the same methodological approach in adult human cohorts (data not shown), this finding indicates that the eukaryotic microbiome in early life consists mainly of fungi, at least in this population. Filamentous fungi and yeasts are common members of the normal gut microbiome, although in a much smaller proportion than bacteria. 49 However, it has been reported that during the first 3 months of life, they are present at a significantly higher diversity than later in life. 11 Furthermore, their role in immune regulation and asthma development has been documented previously in relevant models of murine lung inflammation. In particular, modulation of the mycobiome has been shown to increase the severity of murine lung inflammation. 50 Moreover, antibiotic-induced blooms of Candida species were associated with increased experimental murine lung inflammation. 51,52 This effect has been shown to be mediated by Candida species-produced prostaglandin E 2 , 53-55 an arachidonic acid derivative known play an important role in orchestrating immune interactions relevant in asthmatic patients.
Among the 19 fungal OTUs associated with the AW phenotype (twice as many as the differential bacterial OTUs detected), 7 of them corresponded to the yeast P kudriavzevii. The robustness of the statistical association between this yeast and the AW phenotype prompted us to further validate and confirm its taxonomic identity bioinformatically and through qPCR, especially given the reduced sample size used in the 18S sequencing analysis. We note that there is sufficient variation in this 18S region to identify Pichia species and distinguish between relatives. 56 P kudriavzevii has been described as a commensal member of the gut mycobiome, 57 and to our knowledge, this is the first report of this species implicated in gut microbial dysbiosis associated with atopic disease. Notably, qPCR analysis in all the samples detected not only an increase in P kudriavzevii but also an increase in fungi-specific 18S gene copies, suggesting, also for the first time, that fungal overrepresentation in the infant gut microbiome is associated positively with the development of atopic wheeze. However, given the size of the sample studied and the small but significant differences detected in fungal gene copies and reads, caution should be taken when extrapolating these findings to other populations. These findings should be replicated in other human cohorts before considering gut fungal overrepresentation, as associated with asthma risk.
Evidently, fungal dysbiosis was more pronounced than bacterial dysbiosis in this study. Thus the role of P kudriavzevii and other yeasts should be explored in mechanistic studies using animal models. In addition to being important members of the gut microbiome, fungal species can experience an ecological advantage over bacteria after antibiotic treatment, a very strong risk factor associated with asthma development in most human cohorts. 3 Collectively, these data support our previously proposed concept of a critical window in early life during which microbial dysbiosis can be detected in strong association with a subsequent allergic phenotype predictive of asthma. 9 The nature of this dysbiosis appears to reflect the microbiome composition in a given population, prompting additional analysis in prospective studies around the world and inclusion of fungal analyses in all microbial surveys. Optimized biomarker studies of the microbial taxa and the metabolites involved in regional asthma-associated dysbiosis could help identify infants at risk of asthma before symptoms occur, as well as provide a scientific rationale for future therapeutic strategies aimed at restoring an altered infant gut microbiome.
We thank the ECUAVIDA study for their dedicated work and the cohort mothers and children for their continuing participation. We acknowledge also the directors and staff of the Hospital Padre Alberto Buffoni in Quininde, Esmeraldas Province. Finally, we thank Carlos Sandoval, Matthew Lemay, Evan Morien, Heidi Britton, Miriam Nieves, and Lisa Nikolai for technical assistance.
Clinical implications: Our findings support the importance of early-life infant gut microbial dysbiosis in the development of atopic wheeze and could be used for the development of novel asthma primary prevention strategies through modification of the microbiome in early life.  Fig E3, A and B) and fungal (Fig E3, C and D) diversity in fecal samples from 3-month-old infants who presented with atopic wheeze at 5 years (AW) versus healthy control subjects (CTRL). Sequencing results from all 16S samples and from 33 18S samples (control subjects, 23; AWs, 10) were retained after applying a cutoff of 1000 sequences per sample. The atopic wheeze phenotype did not explain significant changes in these measurements (P 5 .323 and .676 [Mann-Whitney test and permutational multivariate ANOVA] for Fig E3, A and C, respectively). FIG E5. Phylogenetic placement of P kudriavzevii OTUs within Pichia species sequences from the SILVA 123 reference database. 18S ribosomal DNA sequences labeled as Pichia species or Pichiaceae in the 99% SILVA reference database were aligned and then used to construct a reference tree with RAxML. OTUs from this study assigned to Pichia species were placed in this tree with the RAxML evolutionary placement algorithm. All OTUs cluster within a clade of P kudriavzevii.     Direction of the correlation coefficient is in relation to the control group. P values are corrected for false discovery rate. Boldface represents the same OTUs identified as significant by using DESeq2 analysis (see Table E3).