Bronchoalveolar lavage fluid was originally collected from PWCF and disease control (DC) subjects undergoing clinically indicated flexible bronchoscopy. The DC population consisted of children without a confirmed diagnosis of CF. BALF samples were collected under IRB-approved protocols at the University of Minnesota and Children’s Hospital of Colorado. Informed consent, adolescent assent, and parental permission along with Health Insurance Portability and Accountability Act (HIPPA) of 1996 authorization were obtained according to each institution’s IRB rules and regulations. Demographic and clinical data associated with the day of BALF sample collection, including past medical history, indication for bronchoscopy, lung function measurements, culture results, comorbidities, and cytology data were collected from the electronic medical record and the Cystic Fibrosis Foundation Patient Registry.

Flexible bronchoscopy was performed on children in accordance with each institution’s standard-of-care clinical guidelines. Excess BALF that was not needed for clinical testing was stored for future research. A portion of neat unprocessed sample was set aside. If there was sufficient sample volume remaining for processing, another portion was centrifuged at 250 × g for 10 min at 4°C. Following centrifugation, the pellet and supernatant were separated. Neat, pellet, and supernatant were then all stored at −80°C until analysis. Neat and pellet samples were transported overnight on dry ice to Children’s Hospital of Colorado for DNA extraction, quantitative PCR (qPCR) and SSU-rRNA gene sequencing. Both neat and pellet samples were included in the sequencing analysis as has been previously described (Harris et al., 2007; O’Connor et al., 2021).

Total fungal load (TFL) was measured using a published TaqMan qPCR assay (Liu et al., 2012), which targets approximately 350 nucleotides in the V6–V8 variable region of the SSU-rRNA. The nucleotide region corresponds to Saccharomyces cerevisiae positions 1,199–1,549 (Escherichia coli positions 974–1,295). The assay was based on parameters that were published using a cloned fungal fragment (S. cerevisiae) for standards (10³–10⁸ copies, 10-fold dilution series, as standards). TFL was quantified as the logarithmic transformation of the copies per reaction. For fungal profiling, we utilized TFL to adjust template amount (up to maximum volume) when attempting PCR. While all samples were attempted, any with no apparent band after PCR were deemed to have insuffient load for amplification. DNA sequencing was performed using primers based on the fungal load qPCR assay (PCR primer sequences provided in Supplementary Figure 1) and utilizing a modified PCR approach to enable amplification from a larger range of fungal load (Williamson et al., 2021). PCR was performed in 100 μl total reaction volume split across four wells of a 96-well plate. This provides three independent amplification reactions with the sample DNA and an index-specific negative control to monitor bench assembly of PCR. The independent reactions containing template are meant to average initial rounds of amplification to mitigate PCR drift (Wagner et al., 1994), and are reconstituted into a single well (75 μl) amplicon sample for confirmation of amplification and subsequent sequencing steps. One of the sequencing primers was shortened relative to the qPCR assay to generate a 3′ mismatch between the primer and the human SSU-rRNA gene to improve the ratio of fungal DNA amplification. Analyses were performed on negative control samples, which consisted of only assay reagents, to assess background as has been done previously (O’Connor et al., 2021; Williamson et al., 2021).

Illumina MiSeq paired-end reads were aligned to human reference genome hg19 with bowtie2 and matching sequences discarded (iGenomes¹; Langmead and Salzberg, 2012). As previously described, the remaining non-human paired-end sequences were sorted by sample via barcodes in the paired reads with a python script (Markle et al., 2013). The data presented in the study are deposited the NCBI Sequence Read Archive, accession number PRJNA922524 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA922524). The sorted paired reads were assembled using phrap (Ewing et al., 1998; Ewing and Green, 1998). Pairs that did not assemble were discarded. Assembled sequence ends were trimmed over a moving window of five nucleotides until average quality met or exceeded 20. Trimmed sequences with more than one ambiguity or shorter than 250 nt were discarded. Assembled sequences were aligned and classified with SINA (1.3.0-r23838) (Pruesse et al., 2012) using the 10,265 fungal sequences in Silva 138.1NR99 (Quast et al., 2012) as reference configured to yield the Silva taxonomy. Operational taxonomic units (OTUs) were produced by clustering sequences with identical taxonomic assignments. The FungiFQV2 primers generated 9,802,508 sequences for 168 samples (average sequence length: 315 nt; average sample size: 56,771 sequences/sample; minimum sample size: 5,282; and maximum sample size: 583,650). The median Goods coverage score was ≥ 99.65% at the rarefaction point of 5,282 sequences. The software package Explicet (v2.10.5, www.explicet.org; Robertson et al., 2013) was used for display, analysis (rarefied values for alpha diversity metrics).

Persons with cystic fibrosis (PWCF) and DC subjects were compared using Wilcoxon rank sum tests for continuous variables and chi-square or Fisher’s Exact test for categorical variables. To account for differences in sequencing depth, the relative abundance (RA) of each taxon was calculated (100*number of sequences for specific taxon/total number of sequences). The Shannon diversity index was calculated to measure the alpha diversity. Cubic B splines were used to model the relationships between microbial variables, including TFL and diversity, and age by disease group as has been previously described (O’Connor et al., 2021). Spearman correlations were used to assess the correlations between continuous variables. The Morisita-Horn beta diversity measure was used as a distance metric of samples sorted from the same 18S SSU-rRNA mycobiome data set for both Principal Coordinate Analysis (PCoA) and unsupervised hierarchical clustering analysis (HCA) with complete linkage.

One hundred ninety-one samples were collected from 128 (67%) DCs and 63 (33%) PWCF. Patient demographics are summarized in Table 1. CF subjects were significantly older at the time of BALF sample collection and had significantly increased cellular markers of inflammation. Indications for flexible bronchoscopy and primary medical diagnoses for the DC group were previously reported (O’Connor et al., 2021).

Total fungal load as measured by qPCR was higher in PWCF when compared to DC subjects (Wilcoxon, p < 0.01; Figure 1A). This was true across the age spectrum (Figure 2A), with increased fungal burden observed with increased age at collection in both CF (Spearman Coefficient = 0.25, p = 0.05) and DC (Spearman Coefficient = 0.32, p < 0.01) populations. Median TFL was slightly higher in samples with positive fungal cultures (Figure 1B), though it did not show statistical significance [median (IQR): 5.14 (4.55–5.41) in culture positive compared to 4.48 (3.64–5.38) in culture negative, Wilcoxon, p = 0.24]. Culture-negative samples exhibited a wider interquartile range (IQR) of TFL values compared to culture-positive samples. Additionally, a wider and lower IQR was observed for samples with no culture data available (Figure 1B). Background TFL was detected in very low amounts by qPCR, with an average of only 0.9 log10 copies per reaction in negative controls.

Of the 191 samples, 161 (84%) had sufficient load for sequencing, meaning they had an apparent band after PCR. TFL was higher in samples that were successfully sequenced [median (IQR): 4.70 (4.05–5.23) in sequenced samples compared to 2.57 (2.06–3.10) in non-sequenced samples, Wilcoxon, p < 0.01]. Patient demographics for samples that were successfully sequenced are summarized in Supplementary Table 1. CF samples exhibited increased frequency of sufficient amplification (Chi-squared, p < 0.01) with 60 (95%) being successfully sequenced.

Of the successfully sequenced samples, the median (IQR) Shannon alpha diversity measures were 2.64 (2.15–2.08) for CF samples and 2.68 (2.23–3.10) for DC Samples indicating no difference in mycobiome diversity (Wilcoxon, p = 0.29). However, PWCF exhibited slightly decreased diversity at older ages, though these findings were driven by a subset of individuals (Figure 2C). Shannon diversity was not associated with neutrophilic inflammation (Figure 2D). The relationship between Shannon diversity and TFL is shown in Supplementary Figure 2.

The BALF mycobiome compositions as measured by SSU-rRNA are shown in Supplementary Figure 3. Including only taxa with over 1% relative abundance, the most prevalent taxa included Malassezia, which was found in 133 samples (82%), undifferentiated fungi not categorized into specific taxa, which was found in 119 samples (74%), Cladosporium, which was found in 114 samples (71%), and Pleosporales, which was found in 82 samples (51%). Relative abundances of taxa separated by CF status and organized by age are shown in Figure 3. PWCF had higher relative abundance of Candida-Lodderomyces-clade [median (IQR): 1.18 (0.17–5.08)% in PWCF compared to 0.34 (0.12–2.22)% in disease controls, Wilcoxon, p = 0.02] and Aspergillus [median (IQR): 1.26 (0.27–6.91)% in PWCF compared to 0.62 (0.13–2.44)% in disease controls, Wilcoxon, p = 0.07; Supplementary Figure 4].

The mycobiome compositions for the low-diversity samples, defined as those with a Shannon Diversity index under one, are shown in Supplementary Figure 5. The 13 low-diversity CF samples consisted of samples dominated by Clavispora-Candida-clade (n = 1, 8%), Candida-Lodderomyces-clade (n = 4, 31%), Aspergillus (n = 4, 31%), Trichosporon (n = 1, 8%), Microascaceae (n = 1, 8%), Herpotrichiellaceae (n = 1, 8%), and Trichomonascaceae (n = 1, 8%), while the three low-diversity DC samples consisted of samples dominated by unclassified fungi (n = 1, 33%), Clavispora-Candida-clade (n = 1, 33%), and Taphrinomycotina (n = 1, 33%).

Samples with sufficient load for sequencing were clustered into hierarchical groups to evaluate potential associations between samples with similar mycobiomes and clinical variables, the annotated heatmap is shown in Supplementary Figure 6. Samples showed no clustering by CF status or other primary diagnoses. However, notable clusters included those containing high abundance of Aspergillus and Candida-Lodderomyces-clade, which were both associated with lower diversity (Spearman, p < 0.01 and p = 0.03) and higher TFL (Spearman, p = 0.11 and p < 0.01). Additionally, PCoA including both patient samples and background control samples showed no notable differences in the mycobiome composition between disease groups or between samples and negative controls (Supplementary Figure 7).

In the entire set of 191 samples, when compared to DCs, PWCF had higher white blood cells (WBC) in BALF (Wilcoxon, p < 0.01) and increased neutrophilic inflammation, as measured by percent neutrophils (Wilcoxon, p < 0.01). This was true across the range of TFL measured in samples by qPCR (Figure 2B). Percent neutrophils exhibited a slight positive correlation with TFL in DC samples (Spearman coefficient = 0.17, p = 0.07), but showed no correlation in PWCF (Spearman coefficient = 0.04, p = 0.76) which had heightened neutrophilic inflammation even at lower fungal burdens.

Fungal culture data were only available for 28 total BALF samples (15%) and was more frequently available in PWCF (Fisher’s exact, p < 0.01). When analyzed by study site, fungal culture results were more frequently available in samples from the University of Minnesota (Chi square, p < 0.01). When compared to bacterial culture results, fungal culture results were often accompanied by bacterial culture results (Fisher’s exact, p = 0.02) most of which were positive (Fisher’s exact, p < 0.01). Positive fungal cultures were also more frequently associated with positive bacterial cultures (Fisher’s exact, p = 0.046).

The 10 positive fungal cultures included Aspergillus fumigatis (n = 4, 40%), Candida albicans (n = 1, 10%), Scedosporium apiospermum (n = 1, 10%), an Acremonium species (n = 1, 10%), both Candida guilliermondi and Scopulariopsis brumptii (n = 1, 10%), Candida parapsilosis (n = 1, 10%), and an unspecified yeast fungus (n = 1, 10%). Each of the samples that were culture positive for A. fumigatis had DNA classified into the Aspergillus taxon detected in sequencing. All samples that were culture positive for Candida species had DNA classified into the Candida-Lodderomyces-clade taxon detected in sequencing. However, the sample that was culture positive for Acremonium species had no detectable amounts of DNA classified into the Acremonium taxa. Additionally, Scedosporium and Scopularopsis were not detected at all in the sequencing.