- Research article
- Open Access
Human milk fungi: environmental determinants and inter-kingdom associations with milk bacteria in the CHILD Cohort Study
BMC Microbiology volume 20, Article number: 146 (2020)
Fungi constitute an important yet frequently neglected component of the human microbiota with a possible role in health and disease. Fungi and bacteria colonise the infant gastrointestinal tract in parallel, yet most infant microbiome studies have ignored fungi. Milk is a source of diverse and viable bacteria, but few studies have assessed the diversity of fungi in human milk.
Here we profiled mycobiota in milk from 271 mothers in the CHILD birth cohort and detected fungi in 58 (21.4%). Samples containing detectable fungi were dominated by Candida, Alternaria, and Rhodotorula, and had lower concentrations of two human milk oligosaccharides (disialyllacto-N-tetraose and lacto-N-hexaose). The presence of milk fungi was associated with multiple outdoor environmental features (city, population density, and season), maternal atopy, and early-life antibiotic exposure. In addition, despite a strong positive correlation between bacterial and fungal richness, there was a co-exclusion pattern between the most abundant fungus (Candida) and most of the core bacterial genera.
We profiled human milk mycobiota in a well-characterised cohort of mother-infant dyads and provide evidence of possible host-environment interactions in fungal inoculation. Further research is required to establish the role of breastfeeding in delivering fungi to the developing infant, and to assess the health impact of the milk microbiota in its entirety, including both bacterial and fungal components.
Fungi constitute an important yet frequently neglected component of the human microbiome [1, 2] with a possible role in health and disease [3, 4]. Fungal and bacterial colonisation occur in parallel during early life , yet most infant microbiome studies have overlooked fungi. Milk is a source of diverse and viable bacteria , but only a few studies have assessed fungi in human milk [7,8,9,10,11,12]. Additionally, although geographical differences were observed in a recent study of 80 women from 4 countries , the maternal, infant, and environmental determinants of milk fungi (mycobiota) are still mostly unknown.
Bovine and human studies have confirmed the presence of potentially viable fungi in milk using microscopy, culture-dependent, and culture–independent methods [7,8,9,10, 13]. One metagenomics study of human milk suggests that fungi likely constitute 0.5–2% of the milk microbial community . Another study of 76 mother-infant dyads found that maternal age, blood type, antibiotics, vaginal delivery, and infant sex were associated with Candida colonisation of the infant, and in a subset, maternal vaginal and rectal samples were identified as potential origins of this taxon . Bottle feeding and frequent pacifier use was associated with a higher rate of oral Candida carriage in asymptomatic infants [15, 16]. Host genetic variation in major histocompatibility genes were also associated with variations in bovine milk mycobiota . Home characteristics and season are associated with indoor fungi [17, 18] and thus could plausibly influence milk mycobiota; however, mothers’ milk and the home environment (another potential source of milk fungi) have not been examined.
Bacteria-fungi interactions are ubiquitous in microbial communities including the human gut microbiota. These interactions occur through direct physical interactions or several types of molecular communications affecting the composition and function of each respective assemblage . Both bacteria and fungi compositions are altered in different disease conditions including inflammatory bowel disease [20, 21]. Inter-kingdom bacteria-fungi interactions have rarely been investigated in milk [9, 13]. To address these open questions, we analyzed the presence, composition, and determinants of human milk mycobiota at 3–4 months postpartum in a subset of Canadian mother-infant dyads from the CHILD Cohort Study .
In our study population (Table S1) the majority of dyads were urban residents (96%) although population densities ranged widely from 0 to 36,170 persons/km2 (median 3447). 67% of mothers were atopic and 25% delivered by Caesarean section. 18% of infants developed possible or probable asthma by the age of 3 years. Mean ± SD bacterial richness and diversity of milk were 145 ± 46 and 15.2 ± 9.2, respectively (Table S2).
Presence of milk fungi was significantly and independently associated with environmental characteristics, human milk oligosaccharides, and milk bacterial composition
Using a minimum threshold of 1000 reads/sample informed by positive PCR results (Fig. 1a), 58/271 (21.4%) of mothers’ milk samples contained fungi. To address whether milk fungi were more likely acquired externally (e.g. from the environment, maternal skin and/or the infant oral cavity) or internally from the maternal gastrointestinal tract, we assessed associations of fungi presence with environmental vs. maternal or other factors. Both environmental (e.g. study city, population density, and season) and maternal characteristics (e.g. antibiotics and atopy) were associated with fungi presence in univariate analyses (Figs. 1b & 2a). Milk from mothers in Vancouver had the highest presence of fungi (32% vs. 18–21% in Edmonton, Toronto, and Manitoba; unadjusted OR = 2.06, 95%CI 1.11–3.82, for Vancouver vs. other cities; p = 0.021; Figs. 1b & 2a). Vancouver is a large coastal Canadian city with more humid climate, higher annual precipitation, and milder winters compared with the other study sites (Table S3), and thus likely has higher load of environmental fungi . In line with this, season and population density were also associated with fungal presence, with the lowest presence detected in spring and at low population density (Figs. 1b & 2a).
Among the maternal, infant, and early life factors assessed, maternal intrapartum antibiotics and infant antibiotics at the time of sample collection were associated with a higher likelihood of fungi presence (Fig. 2a) indicating that bacteria originating from the mother or infant might be associated with the presence of milk fungi. Indeed, milk bacterial taxonomic clusters, defined previously on the basis of hierarchical clustering of core taxa , were significantly associated with presence of fungi (p < 0.001; Figs. 2a & 4a). Cluster 1 (enriched in Moraxellaceae, Enterobacteriaceae, and Pseudomonadaceae) had the highest presence of fungi (42%) followed by Cluster 2 (enriched in Streptococcaceae, Staphylococcaceae, and Oxalobacteraceae) (30%) compared to Cluster 3 (enriched in Oxalobacteriaceae and Comamonadaceae) (18%) and Cluster 4 (Streptococcaceae and Comamonadaceae) (12%).
We also assessed the association of milk fungi with secretor status and the 19 most abundant human milk oligosaccharides (HMOs), which were previously measured in the same samples . Secretor status is genetically determined by polymorphisms in the fucosyl transferase 2 (FUT2) gene that influences the synthesis of fucosylated HMOs . The impact of maternal secretor status  and HMOs on infant gastrointestinal bacterial composition is well established , but little is known about their impact on milk mycobiota. There is some evidence that certain HMOs can inhibit fungi in vitro , but it is also plausible that other HMOs could be metabolised by fungi and support their growth. Here, we found no association between secretor status and milk fungi presence; however, two FUT2-independent HMOs (disialyllacto-N-tetraose (DSLNT) and lacto-N-hexaose (LNH)) were less abundant in milk containing fungi (Figs. 2a, 3a and b), suggesting that these HMOs might inhibit or be metabolised by milk fungi. Alternatively, these HMOs might indirectly influence the milk mycobiota via modulating the milk microbiota or fatty acids composition .
To identify the determinants of milk fungi while controlling for the above-mentioned factors, we performed unsupervised variable selection by least absolute shrinkage and selection operator (LASSO) and used the LASSO-selected variables to predict the presence of milk fungi. In agreement with the logistic regression results above, LASSO identified season, population density, study city, residential vegetation, bacterial clusters, bacterial outliers (previously identified based on variance in the first principal component score ), DSLNT, LNH, intrapartum antibiotics, and child antibiotics at 3 months as predictors of fungi presence. In addition, LASSO identified maternal asthma and atopy, suggesting that atopic mothers may have distinct skin fungal load and/or composition, which has been observed in other studies . Combined, these selected variables had a relatively good predictive accuracy for presence of fungi (AUC = 0.77, 95% CI 0.70–0.85) (Fig. 3c). In a multivariable logistic regression model adjusted for all of the LASSO-selected factors, DSLNT, LNH, bacterial taxonomic clusters, and child antibiotics at 3 months were the strongest independent predictors of fungi presence (Fig. 2b). As it is plausible that not all fungal genera respond to these variables similarly, we also examined the association of these factors with the presence of specific fungal taxa (see below).
Among other factors we evaluated, fungi presence was not associated with visible mould or moisture levels in the home, older siblings, pet ownership, maternal BMI, or history of infant oral thrush (a fungal infection caused by Candida spp.) in the first 3 months of life (Fig. 2a). We lacked information on breast thrush and mastitis, which are also known to be associated with Candida spp. in milk . Additionally, although others have found that milk collection using a pump and bottle feeding were associated with increased presence of Candida spp. in milk , we did not observe an association between fungi presence and mode of breastfeeding (nursing at the breast vs. pumped milk feeding at least once in the preceding 2 weeks) (Fig. 2a). Breastfeeding exclusivity (i.e. formula supplementation) at the time of sample collection was also not associated with fungal presence (Fig. 2a).
Milk microbiota composition differs in the presence vs. absence of breastmilk fungi
Since we observed that bacterial taxonomic clusters were associated with fungi presence (Figs. 2a & 4a), we further explored the relationship between milk bacteria and mycobiota. In our previous analysis of milk bacteria, we identified “outlier” samples based on variance in the first principal component score  (Figure S1). Outlier status was not associated with any technical or biological variables previously assessed, but here we found that outliers had significantly higher presence of fungi compared to the rest of the samples (43% vs. 19%, OR = 3.17 95%CI 1.49–6.65, p = 0.002) (Figs. 2a & 4a). Additionally, Proteobacteria richness and diversity were lower in samples containing fungi (Fig. 4b), and a similar trend was observed for total milk bacteria. Overall, the relative abundance of Proteobacteria was lower while Firmicutes and Bacteroidetes were higher in the presence of fungi (Fig. 4c). Using linear discriminant analysis , members of Actinobacteria, Bacilli, and γ-Proteobacteria were enriched in the presence of fungi while α-Proteobacteria and β-Proteobacteria were depleted (Fig. 4d).
Milk fungi profile was dominated by Basidiomycota and Ascomycota and tended to differ by study site, maternal antibiotic use, and maternal atopy
Next, among the 58 samples with detectable fungi, we assessed fungal composition and diversity. Ascomycota and Basidiomycota were the dominant fungal phyla in the milk (Fig. 5a) in agreement with previous reports [7, 9]. Composition was quite heterogeneous at the genus level (Fig. 5b). Overall, 12 genera had a minimum mean relative abundance of > 1% (Table S4). The most prevalent fungi were Candida (present in 60% of samples containing fungi, mean 28.6% ± 38.2%% relative abundance), Alternaria (50%, mean 6.9% ± 20.8%), and Rhodotorula (43%, mean 10.3% ± 24.7%). Some samples were dominated by one genus (relative abundance > 50%) - most frequently belonging to Candida (n = 17, 29% of samples), Alternaria (n = 4, 7%), or Rhodotorula (n = 6, 10%). Less prevalent yet dominant genera within some samples included Clavispora (n = 2, 3%), Exophiala (n = 2, 3%), and Penicillium (n = 2, 3%). Other samples contained multiple less prevalent genera (Fig. 5b). Acknowledging the low power due to small sample size, we detected borderline differences in the presence of Alternaria according to intrapartum antibiotics exposure, and Rhodotorula between study cities (Fig. 5c).
Mean fungal richness and diversity at amplicon sequencing variant (ASV) level were 17.9 ± 13.7 and 2.9 ± 2.6, respectively. Fungal diversity was higher in atopic vs. non-atopic mothers (3.6 vs. 1.9, p = 0.044) (Fig. 5d). We did not find any other significant associations between fungal taxonomic structure or diversity with other maternal, infant, early life, environmental, and milk factors (not shown).
Milk fungi exhibit co-exclusion associations with milk bacteria
Finally, we assessed the inter-kingdom associations of abundant milk bacteria and fungi (taxa with > 1% mean relative abundance). Interestingly, bacterial richness was positively correlated with fungal richness (both assessed at the ASV level, r = 0.48, p < 0.001; Fig. 5e) though bacterial and fungal diversities were not significantly correlated (r = 0.21, p = 0.11). In addition, bacterial taxonomic clusters and compositional outliers were significantly associated with fungi α (Figure S2) and β diversity (Fig. 5f). While several bacterial genera were positively associated with each other, indicative of co-presence, we did not detect any associations between different fungal genera. Candida was the only fungal genus associated with bacterial genera; these inter-kingdom associations were negative, suggesting a mutual exclusion relationship (Fig. 5g).
We found that approximately 20% of milk samples contained detectable fungi and provide evidence that fungal presence was associated with environmental characteristics and milk bacteria. Milk fungal taxonomy was consistent with previous reports [7, 9]. Ascomycota and Basidiomycota were the dominant fungal phyla in the milk and three genera of Candida, Alternaria, and Rhodotorula were the most prevalent taxa. Some of the previously reported fungal genera including Saccharomyces and Aspergillus [7, 9] were not among the abundant genera identified in our study, perhaps due to geographical variations and/or other differences in the participants characteristics. To our knowledge, this is one of the largest studies of human milk mycobiota performed to date, and the only milk mycobiota study to examine home environmental characteristics. Our results expand considerably upon existing knowledge about milk fungi, providing evidence for bacterial-fungal interactions within the human milk microbiome.
Are fungi universally present in milk?
It is intriguing that fungi were not present in the majority of samples in our study. Previous culture-dependent and -independent studies in bovine animals [13, 34] and humans [7,8,9,10, 16] have confirmed the presence of fungi in milk. However, while the majority of bovine milk samples contained detectable levels of fungi , the presence of fungi in human milk has varied from 40% (by culture), to 35–80% (qPCR), and 70–100% (sequencing) in different studies [7, 9, 11]. The presence of fungi in our study was lower than previously reported sequencing estimates; potentially due to the higher and more conservative threshold we applied to define fungal presence, other methodological differences (e.g. milk sample collection method, the volume of milk used for DNA extraction, sequencing depth, or target region: ITS2 in our study and 28S rRNA and ITS1 in previous reports), true geographical variations among countries (Canada in our study vs. Spain, South Africa, Finland, and China in previous reports), and/or differences between study cohorts, sampling time, and mode of breastfeeding. Our results are consistent with reports of non-universal fungal presence in the infant gastrointestinal tract over the first year of life . Given that fungi constitute the rare biosphere of the human microbiota , it is conceivable that a greater sampling effort may improve their identification in low biomass samples such as milk.
Origins and maternal determinants of milk fungi
The origins of milk fungi are unclear and it is not known whether milk is evolved to transfer maternal fungi to the infant. The external (retrograde inoculation) and internal (oro-entero-mammary) pathways hypothesized for the origins of milk bacteria [35, 36] could potentially be extended to milk fungi. However, the supporting evidence in general is lacking. Additionally, it is not clear whether fungi exist in the intramammary milk or only following expression potentially originating from maternal skin. Dominance of milk mycobiota by Candida (a common oral fungi ) in a subset of samples in our study and others [7, 10, 11] suggest the infant oral cavity as a potential source of milk fungi even in the absence of symptomatic oral thrush , while dominance of environmental fungi such as Penicillium, Exophiala, and Rhodotorula in other samples suggest contribution from maternal skin or other environmental surfaces.
Although previous studies have found that milk collection using a pump, bottle feeding, and frequent pacifier use were associated with increased presence of Candida spp. in human milk and the infant oral cavity [15, 16], we did not observe any differences in mycobiota presence, diversity, or composition according to mode of breastfeeding (nursing at the breast vs. pumped and bottle feeding at least once in the preceding 2 weeks). It remains possible that pumping has a transient effect on the milk mycobiota that was not detectable in our study.
Our findings suggest the diversity of milk fungi is influenced by maternal characteristics, such as atopy. The direction of association with maternal atopy warrants further investigation, as it is possible that host immunity influences both atopic sensitization and fungal colonization of the skin and/or milk, or alternatively, that environmental fungi influence both maternal atopy and milk mycobiota . Additionally, intrapartum antibiotic exposure was associated with lower presence of Alternaria in milk. Although antibiotics do not directly target fungi, their impact on the overall bacterial composition has been shown to eliminate constraints on fungal (mainly Candida) colonisation and growth . It is therefore plausible that intrapartum antibiotics indirectly modulate the fungal community in the maternal vaginal and/or gastrointestinal tract as well as the infant oral cavity – all of which are potential sources of milk fungi.
To our knowledge, ours is the first study to report associations between specific HMOs and milk fungi, which warrants further exploration. There is preliminary evidence that milk protein content might be associated with milk fungi , but the extent to which other milk components including HMOs and immunomodulatory factors influence milk fungi remains to be elucidated.
Regional and seasonal differences in the presence of fungi
The environmental pool of available species is an important factor influencing the formation and composition of microbial communities . Similar to Boix-Amorós et al., we observed geographical variations in the presence of milk fungi with the highest prevalence of fungi observed in milk from mothers residing in Vancouver (a large coastal Canadian city with humid climate, high annual precipitation, and mild winters) compared with smaller Canadian cities in the Prairies with lower average temperature precipitation. This difference was partly explained by differences in population density and season, but was not related to other characteristics of the outdoor environment (e.g. greenness) or home environment (e.g. mould, dust). Presence of fungi was lower in milk samples collected in spring consistent with seasonal variations reported in goat milk fungi  which could be partly attributed to fluctuations in outdoor  and indoor  environmental fungi levels.
Bacteria-fungi interaction in the milk: direct or indirect interaction?
We observed associations between bacterial composition and presence and relative abundance of fungi, suggesting that some milk bacterial communities were more permissive to fungal presence, or vice-versa. Correlation of bacterial and fungal richness but not diversity suggests that certain factors may impact microbial dispersal in general, and thus enhance both bacterial and fungal richness, whereas the diversity of each kingdom is separately determined by the niche relevant biotic and abiotic factors. To our knowledge, only one previous study has examined the potential interaction of fungi and bacteria in milk, finding a positive correlation between the skin inhabitant genus Malassezia and bacterial load . Our results provide preliminary evidence for antagonistic bacteria-fungi interactions, consistent with evidence that some milk bacteria demonstrate antifungal properties [43, 44], as observed in other environments . Overall, these results could reflect one or more of the following biological relationships: a) some milk bacterial communities are more permissive to fungal presence and proliferation, b) milk bacterial dynamics are influenced by fungi when they are present, and/or c) milk bacterial composition is influenced by the milk environment, which also independently facilitates fungal inoculation or colonisation . More research is needed to determine whether this observation is due to active fungi-bacteria ecological interactions within the milk environment, or perhaps the infant oral cavity.
Strengths and limitations
This study is among the first to assess milk mycobiota in association with home environment characteristics and milk bacterial composition. A key limitation is that we defined the threshold of fungi presence in relation to PCR amplification results. More accurate methods such as quantitative PCR or culture are required to confirm fungi presence. The availability of outdoor and indoor environmental features in addition to host variables allowed us to comprehensively assess determinants of milk mycobiota, although we lacked information on breast thrush (or use of anti-fungal medications) and mastitis, which are known to be associated with Candida spp. in milk . Although the infants were generally healthy at the time of sample collection, this study was slightly enriched for infants who developed asthma later in life and this could limit the generalizability to healthy infants. Finally, as with all culture-independent studies, we do not have information on the viability of milk fungi detected in our study.
In conclusion, we profiled human milk mycobiota in a well-characterised cohort of mother-infant dyads and provide evidence of possible host-environment interactions in fungal inoculation (Fig. 6). This research on the composition and potential determinants of milk fungi gives rise to several new avenues of exploration. For example, the origins of milk fungi are unclear and it is not known whether milk has evolved to transfer fungi to the infant nasopharyngeal and/or gastrointestinal microbiota. Our data suggest that environmental fungal load is a determinant of milk fungi presence, and that sources of milk fungi likely include the maternal skin and/or infant oral cavity and other environmental surfaces which are more directly influenced by environmental fungi. Further studies are needed to confirm these potential sources and determinants of milk mycobiota. The biological nature of the associations we have observed between bacteria and fungi in milk is also unclear, but could involve competition for nutrient sources or biofilm formation. It is possible that bacteria and fungi form complexes in the infant mouth or maternal skin and are translocated together. The role of milk fungi in infant health and disease should be explored since early life exposure to fungi and gastrointestinal mycobiota have been associated with altered risk of asthma and allergy in children [3, 4, 47, 48]. Further research is required to establish the role of breastfeeding in delivering fungi to the developing infant, and to assess the health impact of the milk microbiota in its entirety, including both bacterial and fungal components. Additionally, further investigation is required to assess the role of milk fungi in infant health outcomes.
Women with singleton pregnancies were enrolled between 2008 and 2012 into the general population CHILD birth cohort (n = 3455), from which a subset of 271 mother-infant dyads was selected for this study . The subset was selected from samples with available 16S rRNA microbiota data  and was enriched for children diagnosed with possible or probable asthma. Asthma was diagnosed by an expert study physician at the clinical assessment at age 3 years and was classified for this analysis as “possible or probable asthma” or “no asthma”. Samples with high concentration of DNA based on bacterial 16S rRNA V4 PCR amplification (≥15 ng/μL) were prioritised. Written informed consent was obtained from the participants. The study was approved by the Human Research Ethics Boards at McMaster University, the Hospital for Sick Children, and the Universities of Manitoba, Alberta, and British Columbia.
Milk sample collection
Each mother provided one sample of milk at 3–4 months postpartum [mean (SD) 17 (5) weeks postpartum] in a sterile milk container provided by the CHILD study. To control for differences in the milk composition of fore- and hindmilk  as well as the diurnal variation , a mix of foremilk and hindmilk from multiple feeds during a 24-h period was collected. Hand expression was recommended, but pumping was also acceptable. The sample was collected in a real life situation with no recommended hand or environment cleaning procedures. Samples were refrigerated at 4 °C at home for up to 24 h before being collected and processed by study staff . Samples were stored at − 80 °C until analysis.
DNA extraction and fungal sequencing
DNA was extracted from 1 ml of milk sample as previously described . Milk mycobiota was assessed by sequencing the internal transcribed spacer (ITS) 2 region of the eukaryotic ribosomal gene with modified ITS3/ITS4 primers  on a MiSeq platform (Illumina, San Diego, CA, USA) as previously described . For each sample, the PCR reaction was performed in duplicate, consisting of an initial denaturing step at 95 °C for 5 min followed by 32 amplification cycles at 95 °C for 40 s, 52 °C for 45 s, and 72 °C for 40 s, with a final extension step at 72 °C for 5 min in an Eppendorf Mastercycler pro (Eppendorf, Hamburg, Germany). Sterile DNA-free water was used as negative controls in sequencing library preparation. Fecal samples and mock community containing 8 bacteria and 2 fungi species (Saccharomyces cerevisiae and Cryptococcus neoformans) (Zymo Research, CA, USA) were used as positive controls. PCR amplification was assessed on all samples by agarose gel electrophoresis and visible bands were reported by two independent observations (SM, KF). Duplicate PCR products were pooled and the concentration of double-stranded DNA was quantified by Quanti-iT PicoGreen dsDNA Assay kit (Invitrogen, CA, USA). The concentrations of milk double-stranded DNA were used to normalise the samples 250 ng/μL of DNA and samples were pooled accordingly. For samples with very low DNA concentration, the maximum 10 μL of sample was pooled. Based on previous bovine milk mycobiome work, we anticipated that a subset of samples would be negative for fungi . Therefore, we pooled samples into three pools based on DNA concentration and PCR amplification (high, intermediate, and low DNA concentration). Pooled DNA was then cleaned using DNA Clean & Concentrator (Zymo Research, USA) and the concentration was measured using Qubit dsDNA HS Assay Kit (Invitrogen, CA, USA). To minimise the potential impact of well-to-well contamination for low DNA samples , the final pooled DNA consisted of high, intermediate, and low DNA with proportions of 80:10:10.
Fungal sequencing processing
Overlapping paired-end reads were processed with DADA2 pipeline  using the QIIME 2 v.2018.6 . Unique amplicon sequence variants (ASVs) were assigned a taxonomy and aligned to the 2017 release of the UNITE v.7 reference database at 99% sequence similarity . Demultiplexed sequencing data was deposited into the Sequence Read Archive of NCBI (accession number PRJNA536254).
Initial preprocessing of ASVs was conducted using Phyloseq v. 1.26.1 . Overall, 1421 unique ASVs were detected. ASVs belonging to kingdoms other than Fungi (e.g. plants, n = 334 ASVs) were removed. We obtained a mean (SD) of 16,520 (78,489) and median (IQR) of 52 (13–709) high quality fungal sequencing reads per milk sample, compared with 129,100 (69,044) reads from the positive controls (stool sample and mock community), and 9 (17) reads in negative controls. The mock community was dominated by Cryptococcus neoformans (Basidiomycota phylum) suggesting sequencing bias against Ascomycota. We did not identify any fungal contaminant ASVs using decontam . The majority of samples had very low sequencing reads with 184/271 having less than 300 reads per sample. The threshold applied to define presence of fungi was informed by the presence of visible PCR bands on gel electrophoresis. There was a clear trend between the presence of a PCR band and the number of sequencing reads, with 9% of samples containing 300–500 reads having a band vs. 71% of samples with 1000–5000 and 91% of samples with more than 5000 reads (Fig. 1a). To optimize the retention of samples that were PCR-positive by visual inspection, as well as the elimination of samples with very low sequencing depth, samples with at least 1000 sequencing reads were considered positive for fungi. These samples were rarefied to the minimum sequencing depth of 1000 sequences, resulting in 58 fungi-positive samples containing a total of 625 fungal ASVs. The number of reads for each ASV was relativized to the total sum. Fungal richness (observed taxa) and diversity (inverse Simpson index) were estimated at ASV and genus levels.
Participants were recruited from four Canadian cities: Vancouver (49.2827° N, 123.1207° W, British Columbia), Edmonton (53.5444° N, 113.4909° W, Alberta), Manitoba (49.8951° N, 97.1384° W, Winnipeg and two rural towns), and Toronto (43.6532° N, 79.3832° W, Ontario). Average climate information from 1981 to 2010 was obtained from Environment and Climate Change Canada (http://climate.weather.gc.ca; access date 13 March 2019; used in descriptive but not statistical analyses). Infant sex, birth weight, gestational age, and birth mode were documented from hospital records. Infant feeding, history of oral thrush, and season of milk sample collection were reported by standardized questionnaire. At the time of sample collection (3–4 months), breastfeeding status was classified as “exclusive” (human milk only) or “partial” (human milk supplemented with infant formula or supplementary food). The mode of breastmilk feeding was also reported and classified as “all direct breastfeeding” (nursing at the breast only, with no feeding of pumped milk), or “some pumped milk” (at least one serving of pumped milk in the past 2 weeks) . Population density was estimated for CHILD participants’ homes using 2006 census data on population counts within the smallest available census geographical boundary (250 m) and was categorized as high or low (above or below the median). Census data also provided the rural vs. urban residence. Green space exposure was quantified using the derived normalized difference vegetation index (NDVI) in a 250 m buffer around the mother’s residential addresses. NDVI, with values ranging from − 1 to 1 is an indicator of overall greenness based on land surface reflectance of visible and near infrared parts of spectrum. Time weighted averages across the 12 months postpartum were assigned and categorized into grey (− 1.0 to 0.2), moderate green (0.2 to 0.3), and green (0.3 to 1.0). Home environment characteristics including dust, moisture, and visible mould levels were determined based on questionnaires completed by mothers and a walk-through home assessment by trained research staff from the CHILD Study ; each exposure was categorized into high or low (above or below median). Bacterial V4 16S rRNA gene sequencing was previously performed .
Data analysis was conducted in R v. 3.5.2 . Presence/absence was tested by κ2 and logistic regression for categorical variables and analysis of variance (ANOVA) for continuous variables. Bacterial enrichment based on presence of fungi was assessed by linear discriminant analysis effect size (LEfSe) with default parameters and logarithmic LDA score threshold of 4 . Association of fungal structure with covariates was assessed using ANOVA after centre log-ratio transformation [61, 62]. P values were corrected with Benjamini-Hochberg’s false discovery rate (FDR) method . Milk microbiota outliers were defined as those contributing greater than the median plus twice the interquartile range of the sample variance to the total . LASSO was performed using glmnet package with default parameters to define the optimum lambda with 10-fold cross-validation. Using the optimum lambda, coefficients were estimated for model parameters . While LASSO does not necessarily provide evidence for biological associations, it is a good technique for feature selection as it shrinks the beta coefficients of variables with minimal contribution to zero. This approach minimally increases the bias while enhancing the interpretability of results. LASSO is most suitable when a small to moderate number of moderate-sized effects are expected  (as is the case in our study). Association of LASSO-selected variables with fungal presence was assessed using multivariable logistic regression and the performance was assessed by area under the receiver operation curve (AUC). AUC values correspond to the accuracy of a binary classification with higher numbers indicative of higher accuracy . Association of fungal alpha diversity with host and environmental factors were assessed by Wilcoxon rank sum or Kruskal-Wallis tests. Correlation of bacterial and fungal alpha diversity measures was assessed using Spearman rank correlation following log10 transformation. Fungi β diversity was assessed on Bray-Curtis dissimilarity using permutation ANOVA (PERMANOVA). Data analyses codes will be available upon request.
Inter-kingdom network analysis
Co-occurrence of the most abundant bacterial and fungal genera (> 1% mean relative abundance) were assessed by Spearman rank correlation using CoNet . One hundred bootstrap samples were used to infer pseudo p values. Only edges with correlation scores of > 0.5 and p < 0.05 (Bonferroni-corrected) were retained. The network was visualised in Cytoscape .
Availability of data and materials
The datasets generated and analysed during the current study are available in the Sequence Read Archive of NCBI repository (BioProject accession number PRJNA536254, https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA536254).
Amplicon sequencing variant
Area under the curve
Body mass index
Double stranded DNA
False discovery rate
Fucosyl transferase 2
Human milk oligosaccharide
Internal transcribed spacer
Least absolute shrinkage and selection operator
Linear discriminant analysis effect size
Normalized difference vegetation index
Polymerase chain reaction
Huseyin CE, O’Toole PW, Cotter PD, Scanlan PD. Forgotten fungi-the gut mycobiome in human health and disease. FEMS Microbiol Rev. 2017;41(4):479–511.
Laforest-Lapointe I, Arrieta MC. Microbial Eukaryotes: a missing link in gut microbiome studies. mSystems. 2018;3(2):e00201–17.
Ege MJ, Mayer M, Normand AC, Genuneit J, Cookson WO, Braun-Fahrländer C, Heederik D, Piarroux R, von Mutius E, Group GTS. Exposure to environmental microorganisms and childhood asthma. N Engl J Med. 2011;364(8):701–9.
Arrieta MC, Arevalo A, Stiemsma L, Dimitriu P, Chico ME, Loor S, Vaca M, Boutin RCT, Morien E, Jin M, et al. Associations between infant fungal and bacterial dysbiosis and childhood atopic wheeze in a nonindustrialized setting. J Allergy Clin Immunol. 2018;142(2):424–34.
Ward TL, Dominguez-Bello MG, Heisel T, Al-Ghalith G, Knights D, Gale CA. Development of the human mycobiome over the first month of life and across body sites. mSystems. 2018;3(3):e00140–17.
McGuire MK, McGuire MA. Got bacteria? The astounding, yet not-so-surprising, microbiome of human milk. Curr Opin Biotechnol. 2017;44:63–8.
Boix-Amoros A, Martinez-Costa C, Querol A, Collado MC, Mira A. Multiple approaches detect the presence of Fungi in human Breastmilk samples from healthy mothers. Sci Rep. 2017;7(1):13016.
Jimenez E, de Andres J, Manrique M, Pareja-Tobes P, Tobes R, Martinez-Blanch JF, Codoner FM, Ramon D, Fernandez L, Rodriguez JM. Metagenomic analysis of Milk of healthy and mastitis-suffering women. J Hum Lact. 2015;31(3):406–15.
Boix-Amoros A, Puente-Sanchez F, du Toit E, Linderborg KM, Zhang Y, Yang B, Salminen S, Isolauri E, Tamames J, Mira A, et al. Mycobiome profiles in breast milk from healthy women depend on mode of delivery, geographic location and interaction with bacteria. Appl Environ Microbiol. 2019;85(9):e02994–18.
Amir LH, Donath SM, Garland SM, Tabrizi SN, Bennett CM, Cullinane M, Payne MS. Does Candida and/or Staphylococcus play a role in nipple and breast pain in lactation? A cohort study in Melbourne, Australia. BMJ Open. 2013;3(3):e002351.
Heisel T, Nyaribo L, Sadowsky MJ, Gale CA. Breastmilk and NICU surfaces are potential sources of fungi for infant mycobiomes. Fungal Genet Biol. 2019;128:29–35.
Morrill JF, Pappagianis D, Heinig MJ, Lonnerdal B, Dewey KG. Detecting Candida albicans in human milk. J Clin Microbiol. 2003;41(1):475–8.
Derakhshani H, Plaizier JC, De Buck J, Barkema HW, Khafipour E. Association of bovine major histocompatibility complex (BoLA) gene polymorphism with colostrum and milk microbiota of dairy cows during the first week of lactation. Microbiome. 2018;6(1):203.
Bliss JM, Basavegowda KP, Watson WJ, Sheikh AU, Ryan RM. Vertical and horizontal transmission of Candida albicans in very low birth weight infants using DNA fingerprinting techniques. Pediatr Infect Dis J. 2008;27(3):231–5.
Darwazeh AM, Al-Bashir A. Oral candidal flora in healthy infants. J Oral Pathol Med. 1995;24(8):361–4.
Jimenez E, Arroyo R, Cardenas N, Marin M, Serrano P, Fernandez L, Rodriguez JM. Mammary candidiasis: a medical condition without scientific evidence? PLoS One. 2017;12(7):e0181071.
Sautour M, Sixt N, Dalle F, L'Ollivier C, Fourquenet V, Calinon C, Paul K, Valvin S, Maurel A, Aho S, et al. Profiles and seasonal distribution of airborne fungi in indoor and outdoor environments at a French hospital. Sci Total Environ. 2009;407(12):3766–71.
Ren P, Jankun TM, Belanger K, Bracken MB, Leaderer BP. The relation between fungal propagules in indoor air and home characteristics. Allergy. 2001;56(5):419–24.
Frey-Klett P, Burlinson P, Deveau A, Barret M, Tarkka M, Sarniguet A. Bacterial-fungal interactions: hyphens between agricultural, clinical, environmental, and food microbiologists. Microbiol Mol Biol Rev. 2011;75(4):583–609.
Hoarau G, Mukherjee PK, Gower-Rousseau C, Hager C, Chandra J, Retuerto MA, Neut C, Vermeire S, Clemente J, Colombel JF, et al. Bacteriome and Mycobiome interactions underscore microbial Dysbiosis in familial Crohn’s disease. MBio. 2016;7(5):e01250–16.
Mukherjee PK, Wang H, Retuerto M, Zhang H, Burkey B, Ghannoum MA, Eng C. Bacteriome and mycobiome associations in oral tongue cancer. Oncotarget. 2017;8(57):97273–89.
Subbarao P, Anand SS, Becker AB, Befus AD, Brauer M, Brook JR, Denburg JA, HayGlass KT, Kobor MS, Kollmann TR, et al. The Canadian healthy infant longitudinal development (CHILD) study: examining developmental origins of allergy and asthma. Thorax. 2015;70(10):998–1000.
Amend AS, Seifert KA, Samson R, Bruns TD. Indoor fungal composition is geographically patterned and more diverse in temperate zones than in the tropics. Proc Natl Acad Sci U S A. 2010;107(31):13748–53.
Moossavi S, Sepehri S, Robertson B, Bode L, Goruk S, Field CJ, Lix LM, de Souza RJ, Becker AB, Mandhane PJ, et al. Composition and variation of the human milk microbiome is influenced by maternal and early life factors. Cell Host Microbe. 2019;25:324–35.
Azad MB, Robertson B, Atakora F, Becker AB, Subbarao P, Moraes TJ, Mandhane PJ, Turvey SE, Lefebvre DL, Sears MR, et al. Human Milk oligosaccharide concentrations are associated with multiple fixed and modifiable maternal characteristics, environmental factors, and feeding practices. J Nutr. 2018;148(11):1733–42.
Kelly RJ, Rouquier S, Giorgi D, Lennon GG, Lowe JB. Sequence and expression of a candidate for the human secretor blood group alpha (1,2) fucosyltransferase gene (FUT2). Homozygosity for an enzyme-inactivating nonsense mutation commonly correlates with the non-secretor phenotype. J Biol Chem. 1995;270(9):4640–9.
Lewis ZT, Totten SM, Smilowitz JT, Popovic M, Parker E, Lemay DG, Van Tassell ML, Miller MJ, Jin YS, German JB, et al. Maternal fucosyltransferase 2 status affects the gut bifidobacterial communities of breastfed infants. Microbiome. 2015;3:13.
Verkhnyatskaya S, Ferrari M, de Vos P, Walvoort MTC. Shaping the infant microbiome with non-digestible carbohydrates. Front Microbiol. 2019;10:343.
Gonia S, Tuepker M, Heisel T, Autran C, Bode L, Gale CA. Human Milk oligosaccharides inhibit Candida albicans invasion of human premature intestinal epithelial cells. J Nutr. 2015;145(9):1992–8.
Moossavi S, Atakora F, Miliku K, Sepehri S, Robertson B, Duan QL, Becker AB, Mandhane PJ, Turvey SE, Moraes TJ, et al. Integrated analysis of human Milk microbiota with oligosaccharides and fatty acids in the CHILD cohort. Front Nutr. 2019;6:58.
Gloor GB, Wu JR, Pawlowsky-Glahn V, Egozcue JJ. It’s all relative: analyzing microbiome data as compositions. Ann Epidemiol. 2016;26(5):322–9.
Zhang E, Tanaka T, Tajima M, Tsuboi R, Nishikawa A, Sugita T. Characterization of the skin fungal microbiota in patients with atopic dermatitis and in healthy subjects. Microbiol Immunol. 2011;55(9):625–32.
Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, Huttenhower C. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12(6):R60.
Delavenne E, Mounier J, Asmani K, Jany JL, Barbier G, Le Blay G. Fungal diversity in cow, goat and ewe milk. Int J Food Microbiol. 2011;151(2):247–51.
Fernandez L, Langa S, Martin V, Maldonado A, Jimenez E, Martin R, Rodriguez JM. The human milk microbiota: origin and potential roles in health and disease. Pharmacol Res. 2013;69(1):1–10.
Moossavi S, Azad MB. Origins of human milk microbiota: new evidence and arising questions. Gut Microbes. 2019. https://doi.org/10.1080/19490976.19492019.11667722.
Bandara HMHN, Panduwawala CP, Samaranayake LP. Biodiversity of the human oral mycobiome in health and disease. Oral Dis. 2019;25(2):363–71.
Azevedo MM, Teixeira-Santos R, Silva AP, Cruz L, Ricardo E, Pina-Vaz C, Rodrigues AG. The effect of antibacterial and non-antibacterial compounds alone or associated with antifugals upon fungi. Front Microbiol. 2015;6:669.
Adair KL, Douglas AE. Making a microbiome: the many determinants of host-associated microbial community composition. Curr Opin Microbiol. 2016;35:23–9.
Callon C, Duthoit F, Delbes C, Ferrand M, Le Frileux Y, De Cremoux R, Montel MC. Stability of microbial communities in goat milk during a lactation year: molecular approaches. Syst Appl Microbiol. 2007;30(7):547–60.
Wojcik A, Blaszkowska J, Kurnatowski P, Goralska K. Sandpits as a reservoir of potentially pathogenic fungi for children. Ann Agric Environ Med. 2016;23(4):542–8.
Adams RI, Miletto M, Taylor JW, Bruns TD. Dispersal in microbes: fungi in indoor air are dominated by outdoor air and show dispersal limitation at short distances. ISME J. 2013;7(7):1262–73.
Delavenne E, Mounier J, Deniel F, Barbier G, Le Blay G. Biodiversity of antifungal lactic acid bacteria isolated from raw milk samples from cow, ewe and goat over one-year period. Int J Food Microbiol. 2012;155(3):185–90.
Allonsius CN, Vandenheuvel D, Oerlemans EFM, Petrova MI, Donders GGG, Cos P, Delputte P, Lebeer S. Inhibition of Candida albicans morphogenesis by chitinase from lactobacillus rhamnosus GG. Sci Rep. 2019;9(1):2900.
Bahram M, Hildebrand F, Forslund SK, Anderson JL, Soudzilovskaia NA, Bodegom PM, Bengtsson-Palme J, Anslan S, Coelho LP, Harend H, et al. Structure and function of the global topsoil microbiome. Nature. 2018;560(7717):233–7.
Rowan-Nash AD, Korry BJ, Mylonakis E, Belenky P. Cross-Domain and Viral Interactions in the Microbiome. Microbiol Mol Biol Rev. 2019;83(1):e00044–18.
Seo S, Choung JT, Chen BT, Lindsley WG, Kim KY. The level of submicron fungal fragments in homes with asthmatic children. Environ Res. 2014;131:71–6.
Sharpe RA, Bearman N, Thornton CR, Husk K, Osborne NJ. Indoor fungal diversity and asthma: a meta-analysis and systematic review of risk factors. J Allergy Clin Immunol. 2015;135(1):110–22.
Hytten FE. Clinical and chemical studies in human lactation. Br Med J. 1954;1(4855):175–82.
Hahn-Holbrook J, Saxbe D, Bixby C, Steele C, Glynn L. Human milk as “chrononutrition”: implications for child health and development. Pediatr Res. 2019;85(7):936–42.
Moraes TJ, Lefebvre DL, Chooniedass R, Becker AB, Brook JR, Denburg J, HayGlass KT, Hegele RG, Kollmann TR, Macri J, et al. The Canadian healthy infant longitudinal development birth cohort study: biological samples and biobanking. Paediatr Perinat Epidemiol. 2015;29(1):84–92.
Walker AW. A lot on your plate? Well-to-well contamination as an additional confounder in microbiome sequence analyses. mSystems. 2019;4:e00186–19.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13(7):581–3.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, Fierer N, Pena AG, Goodrich JK, Gordon JI, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7(5):335–6.
Nilsson RH, Tedersoo L, Ryberg M, Kristiansson E, Hartmann M, Unterseher M, Porter TM, Bengtsson-Palme J, Walker DM, De Sousa F. A comprehensive, automatically updated fungal ITS sequence dataset for reference-based chimera control in environmental sequencing efforts. Microbes Environ. 2015;30(2):145–50.
McMurdie PJ. Holmes S: phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8(4):e61217.
Davis NM, Proctor D, Holmes SP, Relman DA, Callahan BJ: Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 2018, 6(1):226. doi: 2https://doi.org/10.1186/s40168-40018-40605-40162.
Klopp A, Vehling L, Becker AB, Subbarao P, Mandhane PJ, Turvey SE, Lefebvre DL, Sears MR, Investigators CS, Azad MB. Modes of infant feeding and the risk of childhood asthma: a prospective birth cohort study. J Pediatr. 2017;190:192–9.
Takaro TK, Scott JA, Allen RW, Anand SS, Becker AB, Befus AD, Brauer M, Duncan J, Lefebvre DL, Lou W, et al. The Canadian healthy infant longitudinal development (CHILD) birth cohort study: assessment of environmental exposures. J Expo Sci Environ Epidemiol. 2015;25(6):580–92.
R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2017.
Palarea-Albaladejo J, Martin-Fernandez JA. zCompositions -- R package for multivariate imputation of left-censored data under a compositional approach. Chemom Intell Lab Syst. 2015;143:85–96.
Gloor GB, Reid G. Compositional analysis: a valid approach to analyze microbiome high-throughput sequencing data. Can J Microbiol. 2016;62(8):692–703.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Statist Soc B. 1995;57(1):289–300.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software. 2010;33(1):1–22.
Tibshirani R. Regression shrinkage and selection via the Lasso. J R Statist Soc B. 1996;58(1):267–88.
Hanley JA, McNeil BJ. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology. 1982;143(1):29–36.
Faust K, Raes J. CoNet app: inference of biological association networks using Cytoscape. F1000Res. 2016;5:1519.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
We are grateful to all the families who took part in this study, and the whole CHILD team, which includes interviewers, nurses, computer and laboratory technicians, clerical workers, research scientists, volunteers, managers, and receptionists. We acknowledge support from the Canadian Urban Environmental Health Research Consortium (CANUE, www.canue.ca) for providing NDVI metrics, indexed to Canadian postal codes. We thank Diana Lefebvre (McMaster University) for managing the CHILD database, Leah Stiemsma (University of British Columbia) for creating the antibiotic variables, and Lorena Vehling (University of Manitoba) for creating the breastfeeding variables. We also thank Anita Kozyrskyj (University of Alberta), James Scott (University of Toronto), and Mohammad Bahram (University of Tartu) for constructive discussions.
The Canadian Institutes of Health Research (CIHR) and the Allergy, Genes and Environment Network of Centres of Excellence (AllerGen NCE) provided core support for the CHILD Study. This research was specifically funded by Children’s’ Hospital Research Institute of Manitoba. Support was also provided by the Heart and Stroke Foundation and Canadian Lung Association, in partnership with the Canadian Respiratory Research Network and AllerGen NCE, the Natural Science and Engineering Research Council of Canada (NSERC) Discovery Program and internal grants from the University of Manitoba. Infrastructure at the Gut Microbiome Laboratory was supported by grants from the Canadian Foundation for Innovation. This research was supported, in part, by the Canada Research Chairs program. MBA is a Fellow of the Canadian Institutes for Health Research (CIFAR) Humans and the Microbiome Program. SM holds a Research Manitoba Doctoral Studentship. HS is supported by Michael Smith Foundation for Health Research Postdoctoral Trainee and Canadian Institute of Health Research Fellowship. These entities had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; and preparation, review, or approval of the manuscript.
Ethics approval and consent to participate
Written informed consent was obtained from the participants. The study was approved by the Human Research Ethics Boards at McMaster University, the Hospital for Sick Children, and the Universities of Manitoba, Alberta, and British Columbia.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Table S1. Characteristics of mother-infant dyads from the CHILD cohort included in this analysis (N = 271) in comparison with previous bacterial analysis (N = 393). Table S2. Milk characteristics of mother-infant dyads for fungal sequencing. Table S3. Comparison of characteristics of mother-infant dyads from the CHILD cohort included in this study across study cities. Table S4. Most abundant fungal genera (> 1% mean relative abundance) in the human milk microbiota among 58 mothers with detectable milk fungi in the CHILD cohort. Figure S1. Milk bacterial composition outlier. Figure S2. Association of milk bacterial clusters and bacterial composition outliers with fungal richness and diversity tested by ANOVA.
About this article
Cite this article
Moossavi, S., Fehr, K., Derakhshani, H. et al. Human milk fungi: environmental determinants and inter-kingdom associations with milk bacteria in the CHILD Cohort Study. BMC Microbiol 20, 146 (2020). https://doi.org/10.1186/s12866-020-01829-0
- Human milk
- CHILD cohort study