Interaction between the nasal microbiota and S. pneumoniae in the context of live-attenuated influenza vaccine

Streptococcus pneumoniae is the main bacterial pathogen involved in pneumonia. Pneumococcal acquisition and colonization density is probably affected by viral co-infections, the local microbiome composition and mucosal immunity. Here, we report the interactions between live-attenuated influenza vaccine (LAIV), successive pneumococcal challenge, and the healthy adult nasal microbiota and mucosal immunity using an experimental human challenge model. Nasal microbiota profiles at baseline are associated with consecutive pneumococcal carriage outcome (non-carrier, low-dense and high-dense pneumococcal carriage), independent of LAIV co-administration. Corynebacterium/Dolosigranulum-dominated profiles are associated with low-density colonization. Lowest rates of natural viral co-infection at baseline and post-LAIV influenza replication are detected in the low-density carriers. Also, we detected the fewest microbiota perturbations and mucosal cytokine responses in the low-density carriers compared to non-carriers or high-density carriers. These results indicate that the complete respiratory ecosystem affects pneumococcal behaviour following challenge, with low-density carriage representing the most stable ecological state.


Taqman Fast virus 1-step master mix (influenza A and B and RSV) TaqMan 2X Universal PCR
Master Mix (all other viruses; Applied Biosystems), 1.25 µL primer/probe mix. Amplification and detection was performed on a StepOnePlus real-time PCR unit (Life Technologies) under the thermal cycling conditions described by the manufacturer. Swabs collected at baseline were subjected to the complete viral panel. We additionally tested samples collected at day 2, 7 and 29 for influenzavirus in participants who received LAIV.

Nasal lining fluid: Luminex analysis and viral detection at d0/d2
Nasosorption samples were collected at each time point. Each sample was centrifuged for 10 min at 3,220xg to separate the nasal lining fluid from the filter. Filter and fluid were then stored separately.
Then, cytokines were eluted from stored Nasosorption filters (Hunt developments) using 100 µL of assay buffer (ThermoFisher) by centrifugation for 10 min at 3,220xg. Samples were centrifuged for 10 min at 16,000xg to clear them prior to acquisition. Samples (~40 µL) were acquired using a 30-plex magnetic human Luminex cytokine kit (ThermoFisher) and analysed on a LX200 (Biorad) with xPonent3.1 software (Luminex Corp) following manufacturer's instructions. A representative subset of 12 cytokines were selected for further analyses. Redundant/co-clustering cytokines were excluded,

5.
whilst considering previous findings from our group 21 . Samples were analyzed in duplicates and Nasosorption samples with a CV > 25% were excluded. Next, eluted sample (~60 µL) and raw nasal lining fluid were pooled (resulting in 80-120µL of sample) and subjected to a viral qPCR panel as described above.

16S-rRNA sequencing
We selected baseline (day -4), day 2, 7 and 29 nasal wash samples (4 time points; see Figure 1) for microbiota analyses. After bacterial DNA isolation, amplicon libraries of the 16S-rRNA gene (V4 region) were generated. Sequencing was executed as previously described 6,22 . PCR amplicon libraries were generated by amplification of the 16S ribosomal RNA gene using barcoded primers directed at the V4 hypervariable region, as previously described 22 . Primer pair 533F/806R was used for amplification. Amplicon pools from samples and controls were sequenced in seven runs using an Illuminia MiSeq instrument (Illumina Inc., San Diego, CA, USA), resulting in paired-end 250 nucleotide reads. We applied an adaptive, window-based trimming algorithm (Sickle, version 1.33) 23 using a quality threshold of Q20 (as opposed to Q30 22 ) and a length threshold of 150 nucleotides to filter out low quality reads/nucleotides. We the number of sequence errors was further reduced by applying an error correction algorithm (BayesHammer, SPAdes genome assembler toolkit, version 3.5.0) 24 . Next, reads were assembled into contigs (PANDAseq, version 2.9) 25,26 and demultiplexed (Qiime version 1.9.1; split_libraries.py) 27 . We removed singleton sequences and identified chimeras using both de novo and reference chimera identification. After removal of chimeric sequences, VSEARCH abundance-based greedy clustering was used to pick OTUs at a 97% identity threshold 28 .
OTUs were then annotated by the Naïve Bayesian RDP classifier (version 2.2) 29 with a classification confidence of 50% (default) 30 and annotations were based on the 97% identity SILVA 119 release reference database 31 .

6.
All analyses were performed in the R version 3.3.0 within R studio version 1.0.136. All figures were created using the ggplot2 R-package and edited using Illustrator CC.
We also provided a detailed schematic on the research questions/associations explored and a data analysis flow chart depicting an overview of the methods used (Supplementary Figure 10).

Variable definitions
In the manuscript describing the initial results of the LAIV-EHPC project, focussing on the effect of LAIV on pneumococcal carriage, results based on both pneumococcal detection methods (i.e. conventional culture and molecular) were presented, underscoring the importance of the increased sensitivity of molecular techniques 1,2 . For this manuscript we therefore decided to test two carriage outcome variables on the basis of nasal washes from day 2, 7 and 9: 1) carriage2 outcome (based on pneumococcal detection using conventional culture only), 'carriers', with a culture positive sample at any point and 'non-carriers', who were culture-negative at all times; and 2) carriage3 outcome (combination of pneumococcal detection using both conventional culture and molecular techniques), coded as 'high-dense carriers' (culture-positive at any point), 'low-dense carriers' (qPCR-positive and culture-negative) and 'non-carriers' (qPCR-and culture-negative at every point). Initial explorative analyses demonstrated higher explanatory power of carriage3 outcome, i.e. the variable incorporating qPCR results. We therefore decided to use this outcome variable throughout the rest of the manuscript instead of carriage2 outcome.

Nasal wash sample selection
In conformity with the original study 1,2 , we excluded samples from volunteers who were randomized to receive LAIV, but presumably did not receive the required dose due a systematic medication dispensing error. In addition, we excluded two samples in which we detected >3,500 reads. A total of 451 samples from 117 volunteers were included. In 101/117 (86.3%) volunteers, all 4 samples were available; in 15/117 (12.8%) volunteers, samples from 3 time points were available and in 1/117 (0.9%) volunteers, samples from 2 time points were available.

Data normalization and filtering
We generated an abundance-filtered dataset by including only those OTUs that were present at or above a confident level of detection (0.1% relative abundance) in at least two samples, retaining 485 OTUs (0.3% of reads excluded) 32 . We generated a rarefied OTU-table at a sequence depth of 3,500   reads, calculated the relative abundance of OTUs and used this table as input for downstream analyses. a-diversity measures were calculated for 100 rarefactions at a sequencing depth of 3,500 reads and averaged. Analysis of composition of microbiomes (ANCOM) was performed on raw, non-rarefied data 33 . We selected OTUs that were present in >25% of the samples (at at least one time point) with a mean relative abundance if present of 0.1%. This way, OTUs were selected for the comparison of baseline microbiota profiles between carriage3 outcomes. For metagenomeSeq-analyses we used the same selection of OTUs, though normalisation was performed using cumulative sum scaling (CSS) 34 .
β-diversity was assed using the Bray-Curtis dissimilarity metric (dissimilarities based on relative abundance of species) and the Jaccard index ('distances' based on presence/absence of species).
Quality control of 16S-rRNA gene amplicon sequencing OTU-tables were corrected for environmental and procedural contaminants on the basis of environmental (n=5) and procedural/laboratory control samples (n=66) that were taken at the moment of nasal wash collection and bacterial DNA isolation, respectively. Within the procedural control samples, we identified and removed those OTUs with a relative abundance of >1% in >5% (n=4) samples 35 , removing 2.5% of all reads. These 28 OTUs included, 25 OTUs with at least a genus level annotation, of which 72.0% were previously reported as contaminating genera by Salter et al 36 .
Based on our environmental control samples, we identified a subset of potential contaminating environmental OTUs that were also known bacterial community members in the upper respiratory tract, notably several streptococcal species, precluding us from simply removing these OTUs from the dataset. Alternatively, we hypothesized and subsequently verified that these OTUs would demonstrate a strong negative association between (log10+1-transformed) raw read counts and (log10-transformed) bacterial density (linear model, beta-coefficient <-0.1 and p<0.05), indicating a larger impact of contaminating species on low density samples. We next screened the OTU-table using these criteria (considering only OTUs present in >30% of samples with a mean read count of 10) in an unsupervised manner, which resulted in a total of 24 unique, potentially contaminating OTUs. Plotting read counts of these OTUs in time suggested the existence of a batch effect, displaying several time intervals associated with varying degrees of contamination. We therefore ran a non-parametric change point analysis (changepoint.np package) to identify the cut-off points of these time intervals (i.e. the moments in time in which a shift in raw read counts was observed, hereafter referred to as 'change points'). Change points that were identified across ≥5 OTUs and thus were likely related to a consistent contaminating batch-effect, were selected and the OTUs in which any of these change points were identified (n=11) were adjusted. We first removed the negative correlation between raw read counts and bacterial density (selecting the residuals from a linear model fitted for each interval and OTU). These residuals (with mean=0) were then shifted to a new mean (y-axis translation) which was based on the (log10+1-transformed) read count mean within high density samples (upper quadrant), in which the amount of contamination was expected to be low/non-existent. This new baseline was adjusted for the total number of reads identified in each sample.

PERMANOVA and non-metric multidimensional scaling (nMDS)
Global microbiota differences between carriage2/3 outcomes at baseline and at subsequent time points were visualised using non-metric multidimensional scaling (nMDS; metaMDS function in the vegan package; trymax=100) 37 based on the Bray-Curtis dissimilarity matrix. Ellipses representing the standard deviation of data points were calculated using the internal veganCovEllipse function. Stressvalues, which indicate how well the ordination captured the high-dimensional data (i.e. a measure of goodness-of-fit), were reported. We tested whether a nMDS-visualisation in a higher dimensional space would decrease the stress of the ordination using a scree plot (1-6 dimensions tested). However, since LAIV was allocated in a randomized manner and administered after baseline, this effect should not be tested for at this time point. At time points following pneumococcal inoculation we assessed the impact of carriage3 outcome (i.e. pneumococcal exposure/colonization) on nasal microbiota, adjusting for vaccination group, presence of any virus at baseline, the interaction between carriage2/3 outcome and both vaccination and the presence of any virus at baseline. We ran these tests both with and without the OTU representing S. pneumoniae (which was excluded before rarefactions) to quantify whether changes in overall microbiota composition were mainly driven by introduction of pneumococcus. P-values and R 2 -values accompanying nMDS-plots were based on the PERMANOVA-tests described above.

a -diversity
We visualised and analyses three measures of a-diversity: the number of observed species, the Shannon index and the Simpson index. For each a-diversity measure, differences between carriage3 outcomes at baseline were assessed using Wilcoxon rank sum tests and at each time point using mixed linear models. For the unstratified analyses, mixed linear models were fitted including the interaction between carriage3 outcome and time point as fixed effects and subject as random effect. Post-hoc tests on contrasts of interest (i.e. differences between carriage3 outcomes at each time point) were performed using the multcomp package 38 . We adjusted for multiple testing using the 'single-step' procedure (multcomp default), except when stated otherwise. To stratify our results for the effect of vaccine, we ran a similar mixed linear model, although this time including the interaction between carriage3 outcome, time point and vaccination group as fixed effects and subject as random effect.
Within vaccine groups, we again extracted the contrasts of interest (i.e. differences between carriage3 outcomes at each time point).

Clustering
Especially since we observed a bimodal relative abundance distribution in some OTUs, we complemented our supervised analyses using metagenomSeq (see below) with a clustering approach.
Individuals were clustered using unsupervised average linkage hierarchical clustering based on the Bray-Curtis dissimilarity matrix, as described before 22,39 . The number of clusters was determined based on the Silhouette and Calinski-Harabasz indices (fpc package) 40 . Clusters consisting ≥10 samples were considered for subsequent analyses. Subsequently, using a random forest algorithm, we identified the OTUs that discriminated most between clusters (Supplementary Figure

Detection of differentially abundant OTUs: metagenomeSeq and ANCOM
Both the metagenomSeq-package 34 and ANCOM-package 33 were used to determine differentially abundant OTUs between carriage3 outcomes at baseline. Analyses were performed without accounting for vaccine, only including carriage3 outcome as predictor. For metagenomeSeq, we additionally added a normalisation factor as predictor (default). For both packages, information on differentially abundant OTUs between carriage3 groups (if applicable within vaccination groups) was extracted using contrast matrices.

Cytokine
Cytokine data were log2-transformed; missing values were imputed with the mean for that cytokine.
Stratified analyses were performed using a linear model including carriage3 outcome, vaccination group and the interaction between carriage3 outcome and vaccination group as independent variables.
Relevant contrasts were extracted using contrast matrices and the multcomp-package. We tested 1) cytokine levels at day 0 and 2) area-under-the-curves (AUCs; day 0 tot day 9) of cytokine levels. To study the relationship between nasal microbiota and the mucosal host immune response we used canonical correspondence analysis (CCA, cca function of the vegan package, 999 permutations) and distance-based redundancy analysis (dbRDA; capscale function of the vegan package, 999 permutations). For both functions, the log10+1-transformed rarefied OTU-table was used as outcome variable. Significant terms (i.e. cytokines) were determined using anova.cca (vegan package). Results of both CCA and dbRDA were visualised in a two-dimensional space (based on the first two axes), including samples (data points), significant terms (i.e. cytokines, arrows), bacterial species most

12.
strongly associated with the first two axes (n=10 species with the highest absolute scores).
Furthermore, carriage3 outcomes were visualised using confidence ellipses (as previously described).
To compare how well the data were separated by carriage3 outcome for ordination methods incorporating both microbiota and cytokine data, versus microbiota data alone, we regressed carriage3 outcome against X-and Y-coordinates for each given sample, stratified by method used (i.e. nMDS [ Figure 2], dbRDA and CCA [ Figure 4]). Beta-coefficients for levels of carriage3 outcome (i.e. highor low-dense carriers and non-carriers) correspond with data separation in X-and Y-directions within each ordination space tested. Beta-coefficients were mean-centered and scaled for comparability across models. The average of the absolute values of these standardized beta-coefficients is a measure of data separation driven by pneumococcal carriage outcome. 42 13.

Supplementary Figure 5 -Random forest classifier taxa of the clusters and cluster distribution in time. (A)
Heatmap of mean relative abundance of the 20 OTUs that most strongly discriminated between microbiota clusters. Cluster membership was determined using average linkage hierarchical clustering based on the Bray-Curtis dissimilarity matrix. Colours correspond with row wise normalized relative abundances (i.e. yellow indicates the maximum relative abundance of that OTU across clusters, deep purple indicates the minimum relative abun-   (4) were detected in low-dense carriers (and to a lesser extent high-dense carriers) compared to non-carriers. * indicates a significant difference after correction for multiple testing. Red, EC+, high-dense carriers, n=49; blue, EC-, non-carriers, n=40 and orange, EC+-, low-dense carriers, n=27).

Carriage 3 outcome
Prevotella (25) Corynebacterium propinquum (3) Eikenella (104) Prevotella (25) Campylobacter rectus (43) 24.   controls (C) were shown. P-values were derived from mixed linear models with subject as random effect. Significant differences between high-dense and non-carriers, between low-dense and non-carriers and between low-dense and high-dense carriers were denoted with *, # and $, respectively. One symbol, p<0.05; two symbols, 0.005 ≤ p < 0.01 and three symbols p<0.005. Red, EC+, high-dense carriers; blue, EC-, non-carriers and orange, EC+-, low-dense carriers. See Supplementary Table 7 for sample size per study day/carriage 3 outcome/vaccine (sample size is the same for each -diversity measure tested). Means (points) and standard errors of the mean were (whiskers) of log 2 -transformed cytokine levels in (pg ml -1 )

26.
were plotted at each time point and stratified by carriage 3 outcome. Separate plots were generated for the LAIV and control groups. For statistical assessment, see Supplementary Table 8. Red, EC+, high-dense carriers; blue (n=18 for LAIV and n=18 for controls for each cytokine and at each time point), EC-, non-carriers (n=7 for LAIV and n=12 for controls) and orange, EC+-, low-dense carriers (n=9 for LAIV and n=8 for controls).

32.
Supplementary Percentages correspond with the proportion of samples available at each time point per vaccine group, out of the total number of participants in that group.

33.
Supplementary The effect of pneumococcal density/carriage 2 outcome was adjusted for month, presence of any virus at baseline (day -4), month of sampling (i.e. seasonal effects) and the interactions between carriage 2 outcome/pneumococcal density and vaccine/presence of any virus at baseline. These interactions were included to properly assess the associations between baseline microbiota and carriage 2 outcome/pneumococcal density, the latter of which could have been impacted by viral presence. Analyses were performed using PERMANOVA. For a detailed assessment on the association between baseline nasal microbiota and carriage 3 outcome, see Table 1. Only OTUs with an adjusted P-value of 0.05 and a -1.5 < log 2 FC < 1.5 were shown. Asterisks denote those OTUs that were selected more than once (i.e. these OTUs were significantly different in more than one contrast). The 'group_sig' column denotes the carriage status outcome (EC+, high-dense carriers; EC-, non-carriers and EC+-, low-dense carriers) a specific OTU was positively associated with. Tests were not stratified for vaccine group. logFC, log 2 FC; adj.P.val, Benjamini-Hochberg adjusted p-values. See also Supplementary Figure 8. EC+, high-dense carriers; EC-, non-carriers and EC+-, low-dense carriers.