Baseline characteristics
The recruitment of the study began on October 8, 2018 and lasted for 2 weeks. We enrolled 48 healthy volunteers with an average age of 26.6 years old (Table 1). The follow-up visits took place from October 2018 to December 2018 (Fig. 1). In total, 221 induced sputum samples were collected (Table 2). All the healthy volunteers resided in Beijing, China and never left Beijing during the study period. Thirty-eight subjects (79.2%) completed the five timepoints. One subject was diagnosed with acute laryngitis after timepoint D0. Six subjects had common cold at timepoint D30, and three at D60. All dropouts were lost to follow-up right after the last visit. The induced sputum samples were analyzed by sequencing of the 16S rRNA genes, yielding 24,027,321 reads after quality filtering with 108,720±107,741 (12,525~719,237) reads per sample on average. Thirty-six negative control samples were sequenced alongside the sputum samples obtaining 4,149±6,288 reads per sample.
Table 1
Demographics of healthy volunteers.
Parameter
|
Azithromycin Group
(n=24)
|
Placebo
Group
(n=24)
|
P Value
|
Male sex, n (%)
|
12(50)
|
12(50)
|
|
Age, yr, median±IQR
|
25±3.00
|
25±3.25
|
0.60d
|
Height, cm, mean±SD
|
168.75±7.57
|
167.21±8.52
|
0.51e
|
Weight, kg, mean±SD
|
62.71±12.27
|
61.27±9.55
|
0.65e
|
Race, n (%)
|
|
|
|
Han
|
23(95.83)
|
22(91.67)
|
1f
|
Others
|
1(4.17)
|
2(8.33)
|
|
Residencea, n (%)
|
|
|
|
Chaoyang District
|
20(83.33)
|
19(79.17)
|
1f
|
Others
|
4(16.67)
|
5(20.83)
|
|
Occupation, n (%)
|
|
|
|
Students
|
15(62.5)
|
19(79.17)
|
0.20h
|
Othersb
|
9(37.5)
|
5(20.83)
|
|
Microbial diversityc
|
|
|
|
a-Diversity, mean±SD
|
|
|
|
Richness
|
414.42±89.63
|
408.5±88.13
|
0.91e
|
Shannon
|
5.99±0.43
|
6.00±0.32
|
0.80e
|
β-Diversity
|
|
|
|
Unweighted Unifrac distance
|
NA
|
NA
|
0.75i
|
Weighted Unifrac distance
|
NA
|
NA
|
0.75i
|
a: All the participants were resided in Beijing, China, most of them lived in Chaoyang District, others in Changping District, Shijingshan District, Haidian District, Dongcheng District and Fengtai District.
b: Other occupations included physicians, company employees, and civil servants.
c: Diversity of the baseline sputum microbiota of the healthy volunteers using samples collected before drug administration.
d: Compared by two-sided Wilcoxon rank sum test.
e: Compared by two-sided Student's t-test.
f: Compared by two-sided Chi-squared test with continuity correction.
h: Compared by two-sided Pearson's Chi-squared test.
i: Compared by PERMANOVA.
IQR is short for interquartile range.
SD is short for standard deviation.
Table 2
Data of sputum samples collected at different timepointsa.
Timepoints
|
Azithromycin Group
(n=114)
|
Placebo Group
(n=107)
|
D0
|
24
|
24
|
D4
|
24
|
23
|
D14
|
24
|
23
|
D30
|
21
|
20
|
D60
|
21
|
17
|
a: The sputum samples were collected before the day of administration of the drugs (D0), the day after the treatment course was completed (D4), 14 days (D14), 30 days (D30) and 60 days (D60) post-dosing. D: day
b: Thirty-eight of the all subjects completed the five timepoints. Six subjects had a cold at timepoint D30. Three subjects had common cold at timepoint D60. One subject was diagnosed with acute laryngitis after timepoint D0. All dropouts were lost to follow-up right after the last visit.
Minimal influence of procedural contaminations on sputum samples sequencing
It was important to identify the impact of contamination on sequencing results, because the consumables and reagents used in the experiments, such as DTT, DNA isolation and library preparation might contain bacterial DNA. The microbial community composition between sputum samples and control samples was quite distinct based on both unweighted UniFrac and weighted UniFrac distance (PERMANOVA, R2=0.38, 0.28; P=0.001,0.001 Additional file 1: Figure S1a-b). The five Zero-radius Operational Taxonomic Units (ZOTUs, 100% sequence similarity) detected in sputum samples with the highest abundance comprised 23.9% of all sequences, while comprising 9.8% of all sequences detected in control samples (Additional file 1: Figure S1c). The most abundant ZOTU detected in control specimens (ZOTU36) accounted for only 0.029% of all sequences detected in sputum specimens (Additional file 1: Figure S1d). We performed microbial ecology analysis using the decontam package whose function was to identify contaminants in sequencing data based on statistical approach[19]. We used both the frequency and prevalence methods with threshold 0.1 and 0.5 respectively and confirmed only 3 contaminants ZOTUs (ZOTU1714, ZOTU2282, ZOTU2462) in our data (Additional file 1: Figure S1e-f, Additional file 2: Table S1). The relative abundance of the 3 ZOTUs were very low in both sputum and negative control samples (comprising 1.96 × 10-4%, 6.24 × 10-5%, 4.58 × 10-5% of all sputum samples sequences and 1.34 × 10-3%, 1.87 × 10-2%, 1.34 × 10-2% of all negative control samples sequences). ZOTU1714, ZOTU2282 and ZOTU2462 only appeared in 2.26%, 3.62% and 1.81% of sputum samples and 2.78%, 8.33% and 2.78% of control samples, respectively. The rare and low-prevalence sequences often had no effect on the microbial diversity[20]. In conclusion, we found few evidences of procedural contaminations influencing the species detected in sputum. After removal of the 3 ZOTUs, a total of 2827 ZOTUs were analyzed.
Stable bacterial DNA burden in healthy volunteers’ airways after azithromycin administration
We quantified the bacterial DNA of every sputum sample using droplet digital
PCR system. There was considerable variation across individual sputum samples (the first quartile and the third quartile, 1.8 x 107~9.7 x 107 copies/g). However, we found no significant difference in sputum bacterial DNA burden across different timepoints in either azithromycin or placebo group (Wilcoxon signed-rank test, Bonferroni adjusted P values (q)>0.05 for all Fig. 2).
The shifts in sputum phylogenetic diversity after antibiotic exposure
The baseline sputum microbiota of healthy volunteers (D0) showed no significant difference between the two groups (two-sided Student's t-test and PERMANOVA, P>0.05 for all Table 1). Immediately after the azithromycin course was completed (D4), both the richness and Shannon index were dramatically reduced compared to D0 (two-sided paired t-test, q=9.799x10-6 and q=3.92x10-6, respectively Fig. 3a-b). By D14, the species richness still remained decreased, but Shannon index significantly increased compared to D4 (two-sided paired t-test, q=1, q=1.062x10-3 Fig. 3a-b), suggesting that the relative abundance of surviving microorganisms affected by antibiotic started to recover. During the two-month follow-up, the species richness significantly increased at D30 compared to D14 (two-sided paired t-test, q=4.366x10-4) and completely recovered at D60 compared to D0 (two-sided paired t-test, q=0.387 Fig. 3a). Meanwhile, we found no significant differences in richness and Shannon index among the healthy volunteers in the placebo group across the five timepoints. (Additional file 3: Figure S2a-b).
Although pairwise PERMANOVA testing using unweighted UniFrac distance showed significant compositional differences between timepoint D0 and timepoints D4, D14, D30 and D60 (R2=0.083, 0.105, 0.081, 0.054; q=0.01, 0.01, 0.01, 0.05), principal coordinate analysis (PCoA) demonstrated that the sputum microbial compositions gradually returned towards their initial composition after profound differences at D4 (Fig. 3c). The microbial profiles recovery was also reflected by a highest unweighted UniFrac distance observed at D4, along with significantly descending magnitude over time compared to the placebo group (Fig. 3d). When we used weighted UniFrac distance accounting for relative abundance to perform pairwise PERMANOVA testing in azithromycin group, the significant change in microbial community composition compared to D0 was identified at D4, but not D14 (R2=0.104, 0.036; q=0.01, 0.92 Fig. 3e).
Microbial taxonomic variation during 60 days’ follow-up
We investigated the microbial taxonomic variation across the five timepoints by comparing the detection rate and relative abundance of the ZOTUs. We involved 294 ZOTUs whose relative abundance was greater than 0.01% and detection rate was greater than 50% at baseline in the analysis. After exposing to azithromycin, the detection rate of 112 ZOTUs decreased at D4. There was a depletion of some families, such as Veillonellaceae (17.86%), Leptotrichiaceae (9.82%), Fusobacteriaceae (8.93%), Neisseriaceae (7.14%), Actinomycetaceae (7.14%), Pasteurellaceae (6.25%), Prevotellaceae (6.25%) and Lachnospiraceae (5.36% Figure 4, Additional file 4: Table S2). By D14, 74 (66.07%) of the ZOTUs whose detection rate decreased significantly between D0 and D4 exhibited significant differences. By contrast, only 34 (30.36%) of the above mentioned ZOTUs showed differences in their relative abundance, suggesting that the loss of species usually recovered later than the abundance of them (Figure 4, Additional file 4: Table S2). Most of the ZOTUs returned to D0 level during the follow-up, but several ZOTUs, such as Fusobacteriaceae ZOTUs, Pasteurellaceae ZOTUs, Neisseria ZOTU and Leptotrichiaceae ZOTUs remained low in detection rate at D30, and had not been back to the initial state by the end of observational period (Fig. 4, Additional file 4: Table S2). We observed an enrichment of the family Streptococcaceae ZOTUs in the relative abundance at D4, immediately post-administration. Most of them rapidly decreased and returned to the baseline level at D14 (Fig. 4, Additional file 4: Table S2).
Among 294 ZOTUs, the detection rate of 204 ZOTUs returned to the D0 level at D14 and only 2.5% of them exhibited significant differences compared to D0 again at D30 (Additional file 5: Figure S3a, Additional file 6: Table S3). The relative abundance of 227 ZOTUs was of no significant difference between D0 and D14 while 9.7% of those ZOTUs had novel variation in the relative abundance at D30 (Additional file 5: Figure S3b-c, Additional file 6: Table S3). At family and genus levels, the relative abundance of 72.9% (n=35) of families and 79.8% (n=71) of genera had been not significantly different compared to baseline at D14 while 7 families and 10 genera had a shift at D30 again (Additional file 5: Figure S3d).
We found no significant microbial variation in the detection rate and relative abundance in the placebo group across the five timepoints (Additional file 7: Figure S4).
Attenuated resilience in response to environmental factors disturbance
Higher biodiversity richness was observed in samples collected at D30 but there was no significant difference in Shannon index between D14 and D30 (two-sided paired t-test, q=1 Fig. 3b). Besides, a significantly compositional shift was also observed (pairwise PERMANOVA, R2=0.081, q=0.01 Fig. 3e), in line with an increased weighted UniFrac distance (0.287±0.068 Fig. 3f) at this timepoint, suggesting that some variation occurred in the abundance of sputum taxa and prevented the recovery. We found that the daily mean concentrations of PM2.5, PM10, AQI values, daily mean humidity and temperature were significantly different (Kruskal-Wallis rank sum test, P<0.05 Additional file 8: Figure S5) among the five periods. The distance-based redundancy analysis (db-RDA) was performed to explain the environmental stresses on the sputum microbiota variation. By D4, not only the antibiotic, but also humidity had a significant influence on the structure of the microbiota (R2=0.027, q=0.004; R2=0.052, q=0.012 Fig. 5a). After the variance decomposition, we found that the antibiotic could individually explain 2.7% of microbiota variance (P=0.001) and humidity could individually explain 2.2% of microbiota variance (P=0.001). Hence, we concluded that administration of azithromycin was related to sputum microbiota, but the environmental factor, humidity, might also shape the microbial communities. By D30, db-RDA showed that the high concentration of PM2.5 significantly contributed to the variability in microbiota composition, but not antibiotics (R2=0.026, 0.013, q=0.004; q=0.088 Fig. 5b), suggesting environmental factors served as a novel shock to disturb the sputum microbiota which had recovered to the near-baseline microbial community composition.
In placebo group, the environmental factors weren’t associated with the sputum microbiota in healthy volunteers (Additional file 9: Figure S6).
The effects of azithromycin on airway microbial interactions within the microbial community network
To explore the changes of microbial interaction after exposure to azithromycin, a bacterial community network analysis was performed. The azithromycin group networks became less complicated and stabilized after the antibiotic exposure regarding decreases in the total number of vertices, edges, connectance, average degree, average clustering coefficient and centralization closeness (Additional file 10: Figure S7). The co-occurrence patterns of microbial communities before and after azithromycin administration were quite distinct. We identified the importance of ZOTUs using closeness centralization scores and degree centralization scores. ZOTUs belonging to families Fusobacteriaceae and Porphyromonadaceae set a central position at D0 and those assigned to families Prevotellaceae, Lachnospiraceae and Veillonellaceae played important roles as well. After exposure to antibiotics, families Veillonellaceae and Leptotrichiaceae came to the first while families Fusobacteriaceae and Lachnospiraceae’s superiority in the network declined. By D60, families Fusobacteriaceae and Lachnospiraceae returned to their original position again (Fig. 6a). Along with the decrease in detection rate or relative abundance of the major species after administration of antibiotics, the synergistic effect in the co-occurrence networks occupied the main status (Fig. 6a). Despite having a few overlapped edges between timepoint D0 and other timepoints networks, the closeness of shared nodes was quite different. Though, by D14, the relative abundance of most ZOTUs showed no significant difference compared to D0, the closeness of them still remained decreased (Fig. 6b-c). In control group, families Fusobacteriaceae, Lachnospiraceae, Prevotellaceae and Veillonellaceae kept occupying the important positions within the framework along with time (Additional file 11: Figure S8a).The co-occurrence patterns of the microbial communities were quite similar, showing that they contained more overlapped edges and the closeness of shared nodes exhibited very small differences among the five timepoints networks (Additional file 11: Figure S8b-c). It was interesting to see that both positive and negative correlation existed in the networks, suggesting that the interactions among the species made the bacterial community stable.
The effects of azithromycin on oral cavity microbiota
The bacterial DNA burden in oral wash samples (the first quartile and the third quartile, 1.1 x 106~4.5 x 106 copies/g) was significantly lower than sputum samples (Wilcoxon signed-rank test, q<0.05 for all Additional file 12: Figure S9). Oral cavity microbiota had lower species richness at baseline compare to sputum microbiota (484.33 ± 89.57 vs 414.42 ± 89.63 two-sided paired t-test, q=0.0033 Additional file 13: Figure S10, Additional file 14: Figure S11). PERMANOVA testing using both unweighted UniFrac distance and weighted UniFrac distance found that the microbial community composition at D0 between the two niches was quite distinct (R2=0.04, 0.05; P=0.029, 0.033 Additional file 14: Figure S11). Azithromycin altered the alpha-diversity and the community composition of the oral microbiota (Additional file 14: Figure S11). However, compared to the sputum microbiota, oral cavity microbiota had a different pattern of change over time (Additional file 14: Figure S11). Bacterial community network analysis showed that the interactions among the species in oral cavity microbiota community recovered earlier than sputum microbiota (Additional file 15: Figure S12, Additional file 16: Figure S13). Detailed explanation can be seen in the supplementary information (Additional file 17: Figure S14, Additional file 18: Figure S15, Additional file 19: Figure S16, Additional file 20: Figure S17, Additional file 21: Supplementary text).