Citizen-science based study of the oral microbiome in Cystic fibrosis and matched controls reveals major differences in diversity and abundance of bacterial and fungal species

ABSTRACT

Introduction: Cystic fibrosis (CF) is an autosomal genetic disease, associated with the production of excessively thick mucosa and with life-threatening chronic lung infections. The microbiota of the oral cavity can act as a reservoir or as a barrier for infectious microorganisms that can colonize the lungs. However, the specific composition of the oral microbiome in CF is poorly understood.Methods: In collaboration with CF associations in Spain, we collected oral rinse samples from 31 CF persons (age range 7-47) and matched controls, and then performed 16S rRNA metabarcoding and high-throughput sequencing, combined with culture and proteomics-based identification of fungi to survey the bacterial and fungal oral microbiome.Results: We found that CF is associated with less diverse oral microbiomes, which were characterized by higher prevalence of Candida albicans and differential abundances of a number of bacterial taxa that have implications in both the connection to lung infections in CF, as well as potential oral health concerns, particularly periodontitis and dental caries.Conclusion: Overall, our study provides a first global snapshot of the oral microbiome in CF. Future studies are required to establish the relationships between the composition of the oral and lung microbiomes in CF.

KEYWORDS:

• Oral microbiome
• Cystic fibrosis
• oral mycobiome
• Candida
• Pseudomonas

Introduction

Cystic fibrosis (CF) is a severe autosomal recessive genetic disease caused by mutations in the CF transmembrane conductance regulator (CFTR) gene [1]. CF is one of the most common rare genetic disorders, particularly in the Caucasian population, affecting one in 2000–3000 newborns in the European Union [2]. The CFTR protein acts as a chloride channel that transports ions across the apical membrane of epithelial cells throughout the body [3]. This channel is involved in the production of several secretions including sweat, digestive fluids, and mucus. If the channel is impaired, these secretions increase their thickness, mostly affecting the function of the lungs, but also other organs such as the pancreas, liver, kidneys, and intestine. Lung infections are common and often develop into chronic and severe, life-threatening forms due to a deficient mucociliary clearance of the thick mucus [4]. Several bacterial and fungal species are commonly associated with chronic respiratory infection of the lower airways in CF, including Aspergillus fumigatus, Candida albicans, Haemophilus influenzae, Staphylococcus aureus, Burkholderia cepacia complex, Stenotrophomonas maltophilia, Achromobacter xylosoxidans and Pseudomonas aeruginosa, with P. aeruginosa playing a major role in the morbidity and mortality of patients [5,6]. Hence, P. aeruginosa is a key pathogen in CF lung disease and has been found to be involved in the progressive obstructive pulmonary disease and bronchiectasis resulting from chronic endobronchial infection [7]. These infections start early in life and progressively increase with patient age, with P. aeruginosa present in the lungs of up to 80% of patients over the age of 18 years [8].

Despite the fact that P. aeruginosa is one of the most widespread and destructive opportunistic pathogens, it does not colonize the airways alone. Microbes commonly present in the oral cavity are also present in sputum from CF patients [9]. In fact, the oral cavity has been proposed as a reservoir of bacteria, both commensal and pathogenic, that can colonize the lower airways due to micro-aspirations [10,11]. A model has been proposed comparing the biodiversity of the respiratory tract to island biogeography [12], in which the mouth and throat, much like the mainland, are sources of relatively high species diversity, whereas the diversity of the airways decreases with the distance from the mouth, just as distant islands display more specific subsets of the mainland’s diversity. An investigation of this notion showed that samples taken from distal lung sites were more distinct from the upper respiratory tract than proximal sites [13]. Thus, there might be a strong connection between the oral microbiome and many of the pulmonary pathogens acting in CF.

Beyond the mouth’s microbiota, the physiology of the oral cavity is also a factor in CF. Both salivary pH and the pH of the airways in CF are typically lower than those of non-CF individuals [14,15], due in part to malfunction of the CFTR protein, leading to defective secretion of the buffer molecule bicarbonate [16]. The airway microbiota in CF have also shown a decrease in alpha diversity, a measure of the number of organisms present in a sample, and a value that decreases more with diminishing lung function [17–21]. Taken together, these factors have important implications on the microbiome and oral health in general, though there is some debate in the literature. The acidic environment may leave CF individuals more susceptible to both dental caries and periodontitis [22,23]. The combination of low pH and low alpha diversity may predispose CF individuals to dental caries in particular, as according to the ‘ecological plaque hypothesis,’ acidogenic bacteria foster an acidic environment and the diversity drops as many species intolerant to the change are unable to grow [24]. Periodontitis, on the other hand, has been shown to present with greater alpha diversity in the oral cavity [25–28]. While some studies suggest a greater risk for dental caries in CF [14,29], a systematic literature review found that CF generally had lower incidence of caries, or no difference [30]. The CFTR protein’s role of pH regulation is vital in odontoblasts, cells which secrete dentin during tooth development, and ameloblasts, cells which deposit enamel, and in fact, abnormal enamel mineralization has been seen in CF as a result of the defective CFTR protein [31,32]. Conversely, it has been suggested that the CFTR protein may actually promote periodontitis, and indeed a study of gingival biopsies showed greater and more widespread CFTR expression in patients with periodontitis compared to healthy controls [33], suggesting that a mutation in the CFTR gene could predispose CF individuals to have better periodontal health. The review of oral health in CF also suggested that treatments for CF, such as antibiotics and inhaled anti-inflammatory medications, may protect CF patients from the colonization of early caries and periodontitis pathogens [30]. One study found that airway pH increased in CF patients after antibiotic treatment [15], and there have been instances of higher pH in CF as compared to controls, which have been attributed to treatment with supplementary pancreatic enzymes [34]. The evidence suggests that the combination of factors typical to CF may inhibit cariogenic and periodontal pathogens in the oral cavity, but still produces an environment with a low pH and low alpha diversity that promotes damage to the outer layers of the teeth.

Here we present the first metabarcoding study, to our knowledge, of the oral microbiome in cystic fibrosis with a comparison to matched control samples. We compare the overall composition of oral bacteria between the groups and explore the differential abundances of bacterial taxa, as well as the presence or absence of yeast and mold species. We also calculate co-occurrence networks to study the underlying ecologies present in the CF and non-CF groups and further analyze the results alongside metadata collected through a citizen science approach.

Materials and methods

Sample collection

We contacted CF individuals and their relatives, through local associations of CF families integrated in the Spanish Federation for CF (https://fibrosisquistica.org/). Additionally, we analyzed samples from matched groups of individuals (see below) from a larger study targeting the general population (www.sacalalengua.org). All participants signed an informed consent form allowing the use of their saliva samples for microbiological research. For participants under the age of 18, the consent form was also signed by one of the parents or a legal guardian. This project was approved by the ethics committee of the Barcelona Biomedical Research Park (PRBB). Samples were collected from January to November 2017. Participants were asked not to ingest any food or beverage (except water) for 1 h before collecting the sample. All donors received clear indications about the sample collection procedure in person, and the collection of the samples was carried out with the assistance of a researcher involved in the project, following a demonstration. All participants responded to a uniform questionnaire (see below), which was adapted for CF in collaboration with CF partner associations. Before collection of the oral rinse, the pH of the saliva was measured using pH test strips (MColorpHast, Merck, range 5.0–10.0; 0.5 accuracy units), the accuracy of which have been previously validated [35]. Saliva samples were collected using a mouthwash as described earlier [35]. In brief, the protocol is as follows: participants rinsed their mouth with 15 mL of sterile phosphate-buffered saline (PBS) solution, for 1 min. Then, they returned the liquid into a 50 mL tube. The samples were then centrifuged at 4,500 g for 12 min at room temperature (r.t.) in an Eppendorf 5430 centrifuge equipped with an Eppendorf F-35-6-30 rotor. The supernatant was discarded and the pellets were resuspended with the remaining PBS, transferred to 1.5 ml tubes and centrifuged at 4,500 g for an additional 5 min at r.t. using an Eppendorf FA-45-24-11-HS rotor. Supernatants were discarded, and pellets were frozen and stored at − 80°C until further analysis.

A total of 31 oral rinse samples were collected from individuals with CF (ages 7–47) during the second edition of ‘Stick out your tongue’, a citizen science project in Spain (SLL2 – from ‘Saca La Lengua’ in Spanish, see http://www.sacalalengua.org [35]), and in collaboration with CF family associations in Spain. Citizen science aims to involve the public in large-scale investigations to obtain a wide range of data as well as to increase public understanding of relatively esoteric scientific undertakings [36,37]. Sample collection was coupled with science communication activities with CF individuals and their relatives, aiming to raise awareness about the microbiome, its role on health and disease, and its potential particularities in CF. The SLL2 project questionnaire about health and lifestyle was adapted to CF with the help of CF associations.

In order to determine the effects of familial relations in the context of this study, we also collected samples from some relatives. Some of the 31 CF individuals were related to others also with CF, including 8 siblings and 2 individuals that were partners. There were 36 other relatives without CF that also participated in the study, which primarily consisted of parents, but also included 2 other partners, 4 siblings, and 1 grandmother. The ages of the siblings with CF ranged from 7 to 41 years old, while the siblings without CF ranged from 22 to 32 years old. The parents ranged from 41 to 71 years old.

DNA extraction and sequencing

The DNA extraction and amplification and sequencing of the V3-V4 region of the 16S ribosomal RNA gene were performed as previously described [38]. Briefly, for sample DNA extraction we used the ZR-96 Fungal/Bacterial DNA kit (Zymo research Ref D6006), following manufacturer’s instructions. Two DNA samples derived from bacterial mock communities from the BEI Resources of the Human Microbiome Project were included as controls: ‘HM-782D’ and ‘HM-783D’. V3-V4 16S primers, PCR conditions, and library preparation for further Illumina MiSeq sequencing in multiplex with 2 × 300 bp reads using v3 chemistry were done following our previous protocols already described [38]. We also included negative controls for the DNA extraction and PCR amplification steps, which provided no visible band or quantifiable amount of DNA by gel visualization or Bioanalyzer, whereas all samples resulted in clearly visible bands after 20 cycles. Twelve such controls were subjected to library preparation and sequenced. Expectedly, these sequenced non-template controls systematically yielded very few reads (a range of 155–1005 reads per sample), in contrast to an average of ~64,000 reads/sample in sample-derived libraries.

Fungal composition analysis

To assess the fungal composition in our samples, we used traditional culture-based methods to enrich for possible fungal species instead of fungal metagenomics technologies (such as ITS amplification), which produced unsatisfactory results in this type of samples in previous studies [38]. In this way, we could overcome the limitations of the low presence of fungi in oral rinse samples and the difficulty to break the fungal cell wall to access the DNA in comparison to bacteria. We used previously optimized protocols in our group [38], which mainly consists of plating 6% of the original sample onto a YPD sterile plate with antibiotics (100 μg/ml of chloramphenicol and 100 μg/ml of ampicillin). After 7 days of incubation at 30°C, we assessed the number of colonies, their phenotypes, and the presence of bacteria. A maximum of 10 colonies were randomly selected per sample and re-grown onto a fresh plate for 24 h in the same conditions. We used MALDI-TOF (Matrix-Assisted Laser Desorption/Ionization-Time Of Flight) mass spectrometry analysis for fungal identification of each colony and, if the results were inconclusive after two attempts, we complemented this analysis with colony PCR to amplify the Internal Transcribed Spacer (ITS) hypervariable region of the 5.8S ribosomal gene (fungal marker) and further Sanger sequencing.

MALDI-TOF analysis was performed with a MALDI Biotyper (Bruker Daltonik MALDI Biotyper) in the Centre for Omics Sciences (COS) in EureCat (Centre Tecnològic de Catalunya, Reus, Spain), following previously described protocols [38]: total proteins from the fungal colonies were extracted in our lab following standard protocols and then samples were sent to COS for the analysis. Samples were deposited in duplicate on a Polished Steel Target Plate (Bruker), coated with the matrix α-Cyano-4-hydroxycinnamic (HCCA) and analyzed with MALDI-TOF/TOF (MALDI ultrafleXtreme, Bruker, Germany). Spectrum identification was performed using the real-time classification software MALDI Biotyper (Bruker, Germany). Thus, for each colony sample, we obtained two alternative scores, which were interpreted as follows: i) 2.3–3.0, highly probable species identification; ii) 2.0–2.299, secure genus identification, probable species identification; iii) 1.7–1.999, probable genus identification; iv) <1.7, not reliable information. Samples identified with scores lower than 2.0 were manually re-analyzed, and its spectrum was classified using the Biotyper Offline Classification software. We considered that the results were consistent with species identification when the best match had a score >2.0 and the second best match was at least >1.7 with at least the same genus as the first one. The remaining samples were re-grown, and the whole process was repeated. If after two more experiments with MALDI-TOF analysis the identification of the fungal species was still inconclusive, we performed ITS amplification and sequencing.

For ITS-Sequencing, we performed a colony PCR from fresh colonies, directly using biomass material from a recently replated colony (maximum of 12–24 hours old). To do this colony PCR we used DongSheng Biotech (DSBio) Taq mix (#2012) and 20 pmol of each primer (ITS1: 5ʹ-TCCGTAGGTGAACCTGCGG-3ʹ; NL4: 5ʹ-GGTCCGTGTTTCAAGACGG-3ʹ) in a total volume of 40 μl, and the following conditions: 5 min at 94°C, then 30 cycles of 30 s at 94°C, 30 s at 55°C, 1.5 min at 72°C, and a final extension of 5 min at 72°C. PCR products were purified with QIAquick PCR Purification Kit (Qiagen), following the manufacturer’s instructions. Samples were sent for Sanger Sequencing to Eurofins Genomics sequencing service, following SUPREMERUN recommendations. The resulting sequences were evaluated with the webtool Blast [39]. Only samples with a high score (>95%) of identity were considered as correctly identified.

In summary, yeast colonies analyzed from a total of 75 samples could be properly identified by MALDI-TOF. For the rest of the samples (17) further analysis with ITS amplification and Sanger sequencing was required: 11 of them could be identified with certainty (percentage of identity >95% in Blast). Additionally, for some samples, we were able to identify the predominant species responsible for bacterial growth on the YPD plates with antibiotics.

Pre-processing of 16S rRNA sequence reads and taxonomy assignment

Sequence reads from fastq files were filtered using the ‘dada2’ R package (version 1.10.1) [40] to produce counts of amplicon sequence variants (ASVs). Low quality reads were first removed by applying the filterAndTrim function with the following parameters: forward and reverse reads were trimmed to lengths of 275 and 230 nucleotides, respectively (truncLen = c(275,230)); the leading 10 nucleotides were trimmed in both reads (trimLeft = c(10,10)); reads with maximum expected errors greater than 5 in both reads were discarded (maxEE = c(5,5)); all other parameters used the default values. The remainder of the pipeline followed the suggestions in the tutorial from the authors of the tool ([41], https://benjjneb.github.io/dada2/tutorial.html). Taxonomy was assigned using the dada2-formatted database of SILVA version 132 [42]. A phylogenetic tree for use in UniFrac distance calculations was generated by following a protocol [43] that uses the ‘DECIPHER’ (version 2.10.2) [44] and ‘phangorn’ (version 2.5.5) [45] R packages. After processing reads with the dada2 pipeline, only those samples with at least 1,000 reads were retained. At the end of this process, there were a total of 1,648 samples, though only a portion of those were used for the analyses of this study, as described below in the section ‘Statistical analyses’.

For analyses regarding the abundances of taxa, a centered log ratio transformation was applied to the ASV counts. Zeros were first replaced with the ‘count zero multiplicative’ method in the cmultRepl function from the ‘zCompositions’ R package (version 1.3.4) [46]. Then centered log ratios were calculated using the codaSeq.clr function from the ‘CoDaSeq’ R package [47,48].

Diversity measures

Alpha diversity measures were calculated using the estimate_richness function from the ‘phyloseq’ R package (version 1.30.0) [49]. For beta diversity measures, both the weighted and unweighted UniFrac distances, which weight dissimilarity between samples by phylogenetic distances between taxa, were calculated using the UniFrac function from the ‘phyloseq’ package. The weighted UniFrac distances give additional weight based on taxa abundances. Bray-Curtis and Jaccard distances were calculated using the vegdist function from the ‘vegan’ R package (version 2.5–6) [50]. As the Jaccard distance is based on the presence or absence of taxa, the decostand function, also from the ‘vegan’ package, was applied to the ASV counts table, using the method ‘pa’ for presence/absence, before the vegdist function was applied. The Aitchison distance was calculated using the aDist function from the ‘robCompositions’ R package (version 2.2.1) [51].

Statistical analyses

When running statistical tests, we first randomly selected representative matched non-CF samples as controls 100 times to ensure consistency in the results. These same 100 sub-samples were used for each of the relevant tests and were matched for geographical location, age, and gender by the following process. Of the 1,617 non-CF samples in the SLL2 dataset, we removed those samples with any other chronic disorder, leaving 1,335 samples. The CF samples came from seven autonomous communities in Spain (Andalucia, Aragón, Cantabria, Catalunya, Madrid, Galicia, and the Basque Country), so from the healthy controls, we first randomly selected two times the proportion of CF samples from each of those locations (i.e. 2x 4/31, 1/31, 11/31, 6/31, 6/31, 2/31, and 1/31, respectively). To ensure a comparable age range, we determined rough age brackets of youth (under 20), adult (20–60), and senior (60 and over), and randomly selected from the geographically matched samples the same proportions of each age group from CF samples (i.e. 10/31 youths, 21/31 adults, 0/31 seniors). Among the CF samples, there were 14 females and 17 males. Thus, a given sub-sampling was finally rejected and reselected if the proportions of males and females were not similar to that of CF samples (i.e. (14±2)/31 females and (17±2)/31 males). From all of the 100 sub-samplings, a total of 352 samples were used as matched controls.

For each of these sub-samples, a number of statistical tests were run with the CF and matched controls together. First, we performed a permutational multivariate analysis of variance (permanova) based on each of the five distance metrics mentioned above using the adonis function from the ‘vegan’ package. The model included the following fixed effects: CF/non-CF, reported use of antibiotics, gender, age, and population of the city/town from which the sample came (as a generalized proxy of both location and lifestyle). There were 21 CF samples that reported antibiotic use and 10 that did not.

Then, to determine differential abundances of taxa and variation in other variables like alpha diversity and pH, we performed a linear model using the function lm from the base R package ‘stats’ (version 3.6.3) [52], again using the same fixed effects as for the permanova test. The abundance values used for these tests were the centered log ratios of the ASV counts, as described above. The Anova function from the ‘car’ R package (version 3.0–7) [53] was used to calculate type-II anova tables, from which p-values were taken for each fixed effect in the model. These p-values were corrected for multiple testing with the p.adjust function from the ‘stats’ package, using the ‘fdr’ method.

Inferred co-occurrence networks

To produce co-occurrence networks, we first filtered out very rare taxa to avoid spurious associations in taxa that do not appear regularly, by using the filterTaxonMatrix function from the ‘seqtime’ R package (version 0.1.1) [54]. We retained those taxa that had at least 15 counts in at least 20 of the 383 samples that included the CF and matched control samples. Then, we calculated the networks for the CF samples and each of the 100 matched control sub-samplings using the spice.easi function from the ‘SpiecEasi’ package (version 1.0.7) [55]. To produce the network figures, we used a tutorial from the authors of the ‘SpiecEasi’ package as a guide [56]. To calculate the Hamming distances between networks, we used the netdist function from the ‘nettools’ package (version 1.1.0) [57].

Results

Increased abundances of airway pathogens and decreased abundances of some periodontal pathogens in CF

Our analyses identified 26 bacterial genera that consistently differed significantly in abundance between CF and controls among the 100 sub-samples (Table 1, see Methods section for explanation of this process) and many of these organisms have implications in the pathogenesis of CF and a number of oral health conditions, which will be examined in the discussion section. The genera found at higher abundance in CF as compared to matched non-CF sets included Chryseobacterium, Microbacterium, Brevundimonas, Blvii28 wastewater-sludge group, Stenotrophomonas, Streptococcus, Rothia, Staphylococcus, Delftia, Comamonas, Scardovia, Desulfobulbus, genera of the family Clostridiales_vadinBB60_group, Mobiluncus, Sphingobacterium, and Mogibacterium (Table 1, Figure 1(b), Supplementary Figure 1). On the other hand, the genera Peptostreptococcus, genera of the Clostridiales family Family_XIII, Alloprevotella, Treponema, Aggregatibacter, Parvimonas, Bergeyella, genera of the order Saccharimonadales, and Fusobacterium were found at lower abundance in CF, as well as, to a lesser extent, Haemophilus (Table 1, Figure 1(b), Supplementary Figure 1). At the phylum level, Firmicutes and Actinobacteria, and unclassified sequences, were found at higher abundances in CF, while Spirochaetes, Fusobacteria and Patescibacteria presented lower abundances in CF (Supplementary Figure 2). There were no taxa from the genus to the phylum level that differed significantly between the non-CF relatives of CF individuals and the same 100 sub-samplings of matched controls. All these differences were not affected by the fixed effects controlled for in the calculations (antibiotic use, gender, age, and population).

Table 1. Significance of differentially abundant taxa and other variables between CF and matched controls. Columns indicate, in this order, the taxonomic level or the type of variable considered, the organism name or the variable name, the tendency of the difference in CF (↗: higher in CF, ↘: lower in CF, permanova results are not directional), the mean adjusted p-value of the statistical comparison between CF and matched controls, and the numbers of matched controls sub-samples for which the test is significant. Rows are ordered first by the tendency in CF samples, with organisms/variables that were greater in CF first, and then by mean adjusted p-value within each variable group

Taxonomic level/variable	Organism/variable	Tendency in CF	Mean adjusted p value	Number of significant sub-sample tests
Genus
	Chryseobacterium	↗	0.00199	100
	Microbacterium	↗	0.00231	100
	Brevundimonas	↗	0.00384	100
	Blvii28 wastewater-sludge group	↗	0.00847	98
	Stenotrophomonas	↗	0.00975	95
	Streptococcus	↗	0.0121	97
	Rothia	↗	0.0129	95
	Staphylococcus	↗	0.0130	92
	Delftia	↗	0.0138	93
	Comamonas	↗	0.0158	91
	Unclassified at phylum level	↗	0.0173	94
	Scardovia	↗	0.0207	90
	Desulfobulbus	↗	0.0257	87
	F.Clostridiales_vadinBB60_group.UCG	↗	0.0311	84
	Mobiluncus	↗	0.0318	83
	Sphingobacterium	↗	0.0325	82
	Mogibacterium	↗	0.0376	82
	Peptostreptococcus	↘	0.00853	97
	F.Family_XIII.UCG	↘	0.0106	97
	Alloprevotella	↘	0.0119	94
	Treponema	↘	0.0221	88
	Aggregatibacter	↘	0.0238	85
	Parvimonas	↘	0.0258	82
	Bergeyella	↘	0.0300	85
	O.Saccharimonadales.UCG	↘	0.0337	83
	Fusobacterium	↘	0.0436	69
	Haemophilus	↘	0.0611	51
Phylum
	Firmicutes	↗	0.00575	100
	Actinobacteria	↗	0.00998	97
	Spirochaetes	↘	0.0157	93
	Fusobacteria	↘	0.0259	89
	Patescibacteria	↘	0.0796	49
Fungi
	Candida albicans	↗	0.0310	82
	Yeasts detected (Yes/No)	↗	0.0765	61
Sialochemistry
	pH	↘	0.0249	86
Alpha diversity
	Shannon diversity index	↘	0.00199	100
	Faith’s phylogenetic diversity	↘	0.0000147	100
	Species richness	↘	0.0000145	100
Permanova (Beta diversity)
	Aitchison	-	0.001	100
	Jaccard	-	0.001	100
	Unweighted UniFrac	-	0.001	100
	Bray-Curtis	-	0.00118	100
	Weighted UniFrac	-	0.0258	90

Figure 1. Cystic fibrosis differs in factors affecting both oral and lung health. (a) Mean relative abundances of 15 most abundant genera in CF samples and matched controls. The remaining genera are grouped together and colored in white. (b) Two of the significantly differentially abundant genera are shown (centered log ratio values of Rothia and Treponema), as well as alpha diversity as calculated by the Shannon diversity index and the salivary pH. In the header of each boxplot, ‘Other effects’ refers to the significance of the fixed effects included in the calculations (antibiotic use, gender, age, and population). The label ‘none’ indicates that none of these had a significant effect on average for the given taxon or variable

Candida albicans more prevalent in CF

We used a culture-based approach paired with proteomics to determine the presence of yeast and mold species in our oral rinse samples (Table 2, see Methods section). Candida albicans was consistently significantly more prevalent CF samples than matched controls (Table 1, mean adjusted p among 100 sub-samples = 0.0310, significant in 82 of 100 sub-sample tests). C. albicans was detected in 17 of the 31 CF samples, while it was found in 50 of the 352 matched controls that were used in all of the sub-sampling tests (p = 0.0000875 for multinomial log-linear model including all CF and matched controls). We also detected yeast, regardless of species, more frequently in CF samples than matched controls, though less consistently than C. albicans in particular. Yeasts were found in 18 of 31 CF samples and in 74 of 352 matched controls (Table 1, mean adjusted p among 100 sub-samples = 0.0765, significant in 61 sub-sample tests, p = 0.000636 for multinomial log-linear model including all CF and matched control samples). These differences were not significantly affected by the fixed effects included in the calculations (antibiotic use, gender, age, and population).

Table 2. Analysis of colonies grown onto YPD + antibiotics plates. The table summarises the number (n) and frequency (%) of samples which formed colonies for a) yeasts, – with b) indicating the mean and the range of the number of colonies for the yeast positive samples – c) mold, d) bacteria, and e-p) identified fungal species as determined by MALDI-TOF per group: subjects with CF (Cystic Fibrosis) and matched control individuals

Colony analysis	CF (n = 31)	Controls (n = 352)
a) Yeasts	18 (58.06)	74 (21.02)
b) Number of yeast colonies, mean (min-max)	106.44 (1–519)	29.80 (1–578)
c) Mold	7 (22.58)	48 (13.64)
d) Bacteria	19 (61.29)	203 (57.67)
e) Candida albicans	17 (54.84)	50 (14.20)
f) Candida guilliermondii	2 (6.45)	4 (1.14)
g) Candida glabrata	1 (3.23)	2 (0.57)
h) Candida parapsilosis	1 (3.23)	5 (1.42)
i) Debaryomyces hansenii	1 (3.23)	3 (0.85)
j) Candida dubliniensis	0 (0)	4 (1.14)
k) Rhodotorula mucilaginosa	0 (0)	2 (0.57)
l) Candida lusitaniae	0 (0)	1 (0.28)
m) Candida tropicalis	0 (0)	1 (0.28)
n) Candida zeylanoides	0 (0)	1 (0.28)
o) Candida intermedia	0 (0)	1 (0.28)
p) Candida spp.,	0 (0)	2 (0.57)

Lower alpha diversity and salivary pH in CF

CF samples consistently had lower alpha diversity than the matched controls, as determined by the Shannon diversity index, Faith’s phylogenetic diversity, which incorporates weights based on phylogenetic distances between taxa, and species richness, which is a count of the unique taxa identified in each sample. The CF samples also consistently had lower pH than matched controls (Figure 1(b)). There was no difference in these alpha diversity values or in pH between the non-CF relatives of CF individuals and the same 100 sub-samplings of matched controls. The fixed effects included in the calculations (antibiotic use, gender, age, and population) did not have significant effects in any of these differences. There was also a significant difference in the overall composition of CF samples compared to matched controls based on a permanova test of the five distance metrics described in the methods, each of which measures distances in a different way (Table 1).

Networks of co-occurring taxa differ significantly between CF and controls

We inferred taxon co-occurrence networks among the 31 CF samples as well as each of the 100 sub-samples to explore underlying differences in the ecology of the oral microbiome in these conditions (Figure 2(a-b)). Using a Hamming distance calculation, which measures the degree to which the connections within two networks differ from each other, we found that the networks of the 100 matched control sub-samples were significantly more similar to each other than they were to the CF network (Kruskal–Wallis p-value = 2.2e-16, Figure 2(c)). The mean Hamming distance between matched control networks was 0.018 ± 0.005, while the mean distance between the CF network and those of the matched controls was 0.025 ± 0.004.

Figure 2. Networks of co-occurring taxa differ significantly between cystic fibrosis and matched control samples. Co-occurrence networks of taxa in (a) the 31 CF samples and (b) all 352 matched control samples. Only those vertices representing taxa mentioned in the text are labeled, as these were the relevant connections that differed between CF and controls. Vertices are colored by the phylum to which they belong. Vertex sizes are proportional to the abundances of those taxa. Edges are colored according to the trend of the association between indicated taxa, where blue edges are positive and red edges are negative. Edge widths are proportional to the strength of the associations. For the matched control network in (b), this comes from a network calculated for all 352 controls together, merely for the sake of visualization. All statistics included in the text are based on the networks from the 100 matched control sub-samples. (c) Distributions of Hamming distances between the CF network and the 100 matched control networks (CF vs mC – yellow) and between each of the matched control networks (mC vs mC – blue). The p-value indicates the significance of the difference between these distributions

There were a number of specific correlations between taxa that differed between the CF network and those of the matched controls, which have potential implications for the oral cavity as a reservoir of microorganisms for the lower airways, as well as for oral health conditions. Some connections that occurred exclusively in the CF network included negative associations between Alloprevotella and Brevundimonas, Lactobacillus and both Fusobacterium and Gemella, and Aggregatibacter and Mesorhizobium, as well as positive associations between Campylobacter and Leptotrichia, Megasphaera and Kingella, and Corynebacterium and Selenomonas. These connections are summarized in Table 3, in the ‘CF networks’ section. There were also some connections that did not occur in the CF network, but did occur consistently in the matched control networks. These are also highlighted in Table 3, in the ‘Matched control networks’ section, which also shows the number of times a given association was significant among the 100 control networks.

Table 3. Associations of co-occurring taxa that differ between CF and matched control networks. The column ‘Association’ lists the co-occurring taxa. The column ‘Direction of correlation’ indicates the direction of that association. The final column shows the number of matched control networks (of the 100 sub-samplings that were performed) in which the correlation was significant

	Association	Direction of correlation	Significant co-occurrences among 100 control sub-samples
CF networks
	Alloprevotella +Brevundimonas	Negative	0
	Lactobacillus + Fusobacterium	Negative	0
	Lactobacillus + Gemella	Negative	0
	Aggregatibacter + Mesorhizobium	Negative	0
	Campylobacter + Leptotrichia	Positive	0
	Megasphaera + Kingella	Positive	0
	Corynebacterium + Selenomonas	Positive	0
Matched control networks
	Ralstonia + Variovorax	Positive	96
	Streptococcus + Gemella	Positive	92
	Prevotella + Veillonella	Positive	92
	Prevotella + Solobacterium	Positive	61
	Prevotella + Atopobium	Positive	50
	Bradyrhizobium +Hyphomicrobium	Positive	86
	Gemella + Granulicatella	Positive	72
	Actinomyces + Stomatobaculum	Positive	62
	Actinomyces + Rothia	Positive	58
	Actinomyces + Atopobium	Positive	46
	Rothia + Granulicatella	Positive	55
	F.Family_XIII.UCG + Parvimonas	Positive	50
	Treponema + Filifactor	Positive	47
	Treponema + Dialister	Positive	33

Discussion

Our results provide further evidence for the relationship between the microbiota of the oral cavity and the lower respiratory tract and suggest a potential reservoir function of the oral microbiome [10–13]. We found a number of organisms that have been associated with lung infections in CF at higher abundance and prevalence in the oral rinse samples of CF individuals as compared to matched controls (Table 1), including the bacterial genera Chryseobacterium [58–60], Microbacterium [61,62], Brevundimonas [60,63,64], Stenotrophomonas [19,21,65–67], Streptococcus [21,68–71], Rothia [68,71,72], Staphylococcus [21,65,68,73], Delftia [60,64,74–76], Comamonas [19,59,64,76], Scardovia [64,76,77], Mobiluncus [76,78], Sphingobacterium [60,79], Mogibacterium [76,80], and the fungal species Candida albicans [68,81–84]. In addition, we found significant differences in the overall composition of the oral microbiome based on CF using a permanova test on five different distance metrics, each of which focuses on different aspects of the composition. Taken together, this highlights the strong differences across the entirety of the oral microbiomes of CF and non-CF individuals. Although we did not find a significant difference in the oral cavity between CF and non-CF in the abundance of Pseudomonas, the genus of the primary infective agent in CF, P. aeruginosa, the differences in these other taxa highlight the dramatic shift in the microbial equilibrium in CF. There is, of course, the possibility that the increased abundances of these organisms results from colonization of bacteria originating in the lower airways, but we feel that it is more likely that the mouth acts as a source of bacteria for the lungs, which become more susceptible to infection under the conditions of CF. For one, all of the organisms mentioned here are detected in both CF and non-CF oral cavities, which would preclude their appearance as a result of CF lung infections. Moreover, the absence of some of the primary infectors in CF in both sample groups, such as Burkholderia or Achromobacter, and the equivalent abundances of Pseudomonas would support the directionality we have suggested.

In fact, P. aeruginosa may be one of the few pathogens prominent in CF lung infections for which the abundance in the oral cavity is not an indicative biomarker, though it interacts with a number of oral species in the lung. Along with some of the other most significant CF pathogens like Staphylococcus and Stenotrophomonas, Pseudomonas is actually a late colonizer of the lungs. These taxa dominate the lung microbiome in CF after the age of 6, while the lungs of CF patients under the age of 2 are primarily composed of oral commensal species from the genera Streptococcus, Prevotella and Veillonella [21]. When Streptococcus colonizes before Pseudomonas, the oral species S. salivarius is able to inhibit the growth of gram-negative bacteria, such as P. aeruginosa, by the production of lactic acid which disrupts their cell membranes [85,86]. This effect may inhibit the initial spread of some CF pathogens to the lower airways, but it primarily acts within the mouth and throat [87], and so other mechanisms allow growth in the lungs as the patient ages. Even as some Streptococcal species inhibit Pseudomonas in the upper airways, others can promote its growth in the lower airways. One model that has been proposed for the physiology of CF lung infections suggests that 2,3-butanediol produced by oral and airway Streptococcus is abundant due to increased fermentation under more anaerobic conditions. The 2,3-butanediol can act as a buffer for the decreased pH, but then becomes a carbon source for both P. aeruginosa and Rothia mucilaginosa, and adds to the virulence of P. aeruginosa by allowing it to produce more reactive oxygen species and to provide additional electron acceptors to other anaerobic CF pathogens [71].

Rothia mucilaginosa is an oral commensal species and perhaps the most closely tied to P. aeruginosa in CF lung infections as it forms similar biofilms [88], may do so alongside P. aeruginosa [72], and can act as a source of metabolites necessary for the production of glutamate by P. aeruginosa, which is used as a component of its cell wall and may increase its virulence [89]. R. mucilaginosa has been associated with a decline in lung function in CF [68] and has been suggested to adapt very efficiently to the CF lung environment, to the point that individual CF patients may have unique strains of the species, as seen with P. aeruginosa and Staphylococcus aureus [72]. It is able to do so by the use of extracellular lactate, which is present at higher concentrations in CF individuals with poor lung function [90], in order to undergo fermentation under anaerobic conditions, as well as adapting to resist antibiotics, protect against foreign nucleic acids, and to take advantage of the Fe²⁺ produced by the reduction of Fe³⁺ by both P. aeruginosa and S. maltophilia in the low-pH CF lung [72].

Desulfobulbus is one of the only taxa found at higher abundance in the oral cavity of CF individuals in this study that, to our knowledge, has not been reported as a CF lung pathogen, though the evidence suggests that it may be well suited to a CF environment. Desulfobulbus is a sulfate-reducer [91], which is relevant because P. aeruginosa is able to use sulfate from mucins in the airways [92], and the resulting volatile sulfur compounds may be utilized by the fungus Aspergillus fumigatus, a common CF pathogen [93]. In the same vein, it has been shown that there is increased sulfation of mucus glycoproteins in CF [94,95].

The frequent treatment of CF with antibiotics may result in the development of antibiotic resistance in a number of pathogens, with evidence that rates of resistance are increasing in some species [66], including Candida albicans [96], also found at greater prevalence in CF samples in our study. C. albicans has been implicated in CF lung infections [81,82,84], and in particular, it has been associated with a substantial decline in lung function [68,83], as well as pancreatic insufficiency, exacerbation of lung infection, lower BMI and lower percentage of predicted forced expiratory volume in 1 second (FEV₁) [83].

Although in this study we did not have specific information regarding the oral health status of the donors, as they were taken from a larger exploratory study that was not focused explicitly on oral diseases (SLL2 – from ‘Saca La Lengua’ in Spanish, see http://www.sacalalengua.org [35]), the results presented here broadly support the notion seen throughout the literature that, in CF patients, there is a lower incidence of periodontitis, though they develop an acidic oral environment prone to harm the enamel and dentin of the teeth. Factors that may indicate lower incidence of periodontitis in CF include low alpha diversity [25–28], as seen in our data, malfunction of the CFTR protein [33], and generally higher use of anti-inflammatory and antibiotic medications [30]. Indeed, nearly all of the genera found at lower abundance in CF samples here have been associated with periodontal disease. Many of these are members of the various bacterial complexes which form biofilms at different stages of the disease [97–103]. Some of the organisms found at higher abundance in CF samples in this study have been implicated in dental caries and carious lesions on the teeth [104–116] and, conversely, others found at lower abundance in CF samples have been linked to the absence of caries and tooth damage [109,117–121].

Alloprevotella presents an interesting case that may merit further study on its own regarding CF. Species in this genus reduce nitrate in the saliva which contributes to the anti-inflammatory response to periodontitis [122–124], so its lower abundance in CF here could be linked in part to lower incidence of periodontitis. However, the nitrate reduction in CF patients seems to be a more complicated process. When it is metabolized to nitric oxide (NO), it has an anti-inflammatory effect in the airways of CF individuals. Some studies have shown that precursors to NO, like nitrate and nitrite, are higher in the saliva and exhaled breath condensate (EBC) of CF patients than in those of controls, but nonetheless, the amount of exhaled NO was lower in CF [125,126]. One proposed explanation is that NO may be produced normally in CF, but its diffusion is inhibited in the thick mucus produced in CF airways [127], though these first two studies suggest that there is an impairment in the formation of NO in CF patients. Another study found that increasing the intake of dietary nitrate led to an increase in exhaled NO as compared to placebo treatments [128]. From this information, we cannot extrapolate to determine the exact mechanism in the impairment of the NO cycle in CF, but the low abundance of the nitrate-reducing genus Alloprevotella forms a link to this process that may warrant a deeper investigation.

Significant differences in many co-occurrences of taxa among the CF samples and the 100 matched control groups suggest underlying ecological differences in these two conditions. The Hamming distance calculations showed that the 100 matched control networks were more similar to each other than to the CF network (Figure 2(c)), and we also highlighted a number of particular co-occurrences that were specific to either the CF or control networks, and these may have implications for the connections between the mouth and lung in CF, as well as for oral health conditions. In the CF network (Figure 2(a), Table 3), the negative association between the genera Alloprevotella and Brevundimonas follows with the results from our study and the literature, wherein Alloprevotella was lower and Brevundimonas was higher in CF [60,63,64]. It is unclear whether this connection is related to the NO cycle in CF, but at least some strains of Brevundimonas species are nitrate reducers [129], so it may be that Brevundimonas is able to outcompete Alloprevotella in the CF environment. There may be a more concrete explanation for the negative associations between Lactobacillus and both Fusobacterium and Gemella, which are periodontal pathogens [97,130,131]. Many studies explore the use of Lactobacillus species as probiotic treatments and have shown that they can inhibit periodontal pathogens [132] by co-aggregating with them, promoting immune responses and interrupting biofilm formations [133]. Campylobacter and Leptotrichia also have a meaningful connection in the CF oral microbiome, as they can be found in association in dental caries [134], where they can metabolize sugars to lactic acid, aiding in the process of cariogenesis [135,136], and we have already mentioned that extracellular lactate levels are higher in CF sputum [90]. Megasphaera is sometimes associated with CF [137,138] and has been considered a member of the core microbiome of the lower lung, regardless of CF status [139]. This genus can work to neutralize acidic conditions that damage teeth [140] and Kingella is acidogenic [141], so their co-occurrence in the CF network may indicate that they complement each other in the acidic CF oral cavity. The association between Corynebacterium and Selenomonas in CF is reasonable, as both have been implicated in CF [142–144].

Among the matched control samples, the co-occurrence networks primarily support our speculations on the underlying mechanisms that might lead to lower incidence of periodontitis and greater damage to enamel in CF patients as compared to non-CF controls. Prevotella and Veillonella are among the four most abundant genera in our dataset and were significantly associated in 92 of the matched control networks. Prevotella, Veillonella, and Solobacterium are all periodontal pathogens [98,145,146], and Prevotella melaninogenica and Solobacterium moorei, whose genera were significantly associated in 61 matched control networks, utilize cysteine to produce hydrogen sulfide (H₂S), resulting in halitosis [147,148]. Cysteine is a precursor to glutathione [149], a peptide that helps to protect the lung from oxidants, and is found at lower levels in the lungs of CF patients [150], and acetylcysteine, a prodrug to cysteine, has been used to improve lung function in CF patients [151]. So perhaps there is a connection between greater overall availability of cysteine in non-CF individuals and the production of H₂S by organisms like P. melaninogenica and S. moorei. Veillonella species also produce volatile sulfur compounds that cause halitosis [146]. Streptococcus and Gemella are among the most abundant taxa present in the oral microbiome (Figure 1(a)) so it is not strange that they would have a significant association in 92 of the 100 control networks, though it is curious that they did not in the CF network. As mentioned above, Streptococcus was significantly more abundant in our CF samples and has been associated with CF lung infections [21,68–71], and Gemella may play a role in exacerbations of CF infections [152]. Nonetheless, they both have some involvement in periodontitis as well [130,131,153,154], and species of Streptococcus in particular make up the ‘yellow-complex’ of periodontal pathogens, which are early colonizers in that disease [145,155]. Somewhat unexpected associations also occurred exclusively in the matched control networks between Gemella and Granulicatella (occurs in 72 of the 100 matched control networks), and between Rothia and Granulicatella (55 of the 100 control networks), as Granulicatella has also been associated with CF [73,156,157] and caries due to its acidogenic nature [158]. The significant associations seen exclusively in our matched control networks involving Ralstonia and Variovorax [159–162], Atopobium [163,164], Bradyrhizobium [165], Hyphomicrobium [166], Actinomyces [167], Stomatobaculum [168], Rothia [102,163,164], Clostridiales Family_XIII [101], Parvimonas [97,101,169], Treponema [97,101], Filifactor [101,170,171], and Dialister [172,173], all relate to organisms implicated in periodontitis.

Conclusion

In this study, we found significant differences in the oral microbiomes of CF and non-CF individuals from all around Spain, which have implications for the potential of the oral cavity to act as a reservoir of microorganisms for the lower airways, as well as for oral health in CF. Differential abundances of bacteria and fungi follow trends in the literature regarding CF lung infections, presenting similar microbial activity in the oral cavity. We highlighted underlying physiological differences that are apparent in the co-occurrences of taxa among CF samples and matched control samples, which add greater evidence to those trends discussed here. These results provide a snapshot of the unique composition of the CF oral microbiome and how it relates to oral health and the lower airways.

Author Contributions

Conceptualization, T.G. and J.R.W.; formal analysis, J.R.W.; validation, J.R.W.; visualization J.R.W., T.G.; methodology, E.S., S.I.G., E.C., E.K., L.C. L.A.B., N.A.S., M.A.T., A.P.S., A.B., C.C., J.P., T.G., M.A.; data curation, L.C., J.P.; writing original draft, J.R.W., S.I.G., E.S., E.C., T.G.; supervision, T.G.; project administration, T.G.; funding acquisition, T.G.; resources T.G., J.P., J.H. All authors have read and agreed to the published version of the manuscript.

Disclosure of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Acknowledgments

We are thankful to all citizens that participated in the second edition of the “Saca La Lengua’’ project by contributing samples and sharing ideas. In particular, for the work described here we are extremely thankful to the Spanish Federation of Cystic Fibrosis (www.fibrosisquistica.org) CF and their affiliated local associations around Spain that actively participated in the project. Only with their effort are studies like this possible. The authors acknowledge the CRG Genomics Core Facility, CRG Bioinformatics Core Facility, CRG biomolecular screening and protein technologies unit, CRG communication and public relationships department, and UCT ICTS High Performance Computing unit for providing access to the computing facilities. CRG authors acknowledge the Spanish Ministry for Economy, Industry and Competitiveness (MEIC) for the EMBL partnership, and Centro de Excelencia Severo Ochoa.

Data availability statement

The fastq files for the samples used for the analyses in this study were uploaded to the Sequence Read Archive (SRA) with the BioProject accession number PRJNA667,146 and can be found here: http://www.ncbi.nlm.nih.gov/bioproject/667146.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

The project was financed by the CRG through Genomics and Bioinformatics Core Facilities funds, and by the EduCaixa programme through funds from the Fundación Bancaria ‘La Caixa’, with the participation of the Center for Research into Environmental Epidemiology (CREAL), and the ‘Center d’Excellència Severo Ochoa 2013-2017’ programme (SEV-2012-02-08) of the Ministry of Economy and Competitiveness. Eppendorf, Illumina, and ThermoFisher sponsored the research by donating some materials and reagents. TG group acknowledges support from the Spanish Ministry of Science and Innovation for grant PGC2018-099921-B-I00, cofounded by European Regional Development Fund (ERDF); from the CERCA Programme/Generalitat de Catalunya; from the Catalan Research Agency (AGAUR) SGR423; from the European Union’s Horizon 2020 research and innovation programme under the grant agreement ERC-2016-724173; from Instituto de Salud Carlos III (INB Grant PT17/0009/0023 - ISCIII-SGEFI/ERDF).