Sinonasal Microbiome Sampling: A Comparison of Techniques

Background The role of the sino-nasal microbiome in CRS remains unclear. We hypothesized that the bacteria within mucosal-associated biofilms may be different from the more superficial-lying, free-floating bacteria in the sinuses and that this may impact on the microbiome results obtained. This study investigates whether there is a significant difference in the microbiota of a sinonasal mucosal tissue sample versus a swab sample. Methods Cross-sectional study with paired design. Mucosal biopsy and swab samples were obtained intra-operatively from the ethmoid sinuses of 6 patients with CRS. Extracted DNA was sequenced on a Roche-454 sequencer using 16S-rRNA gene targeted primers. Data were analyzed using QIIME 1.8 software package. Results At a maximum subsampling depth of 1,100 reads, the mean observed species richness was 33.3 species (30.6 for swab, versus 36 for mucosa; p > 0.05). There was no significant difference in phylogenetic and non-phylogenetic alpha diversity metrics (Faith’s PD_Whole_Tree and Shannon’s index) between the two sampling methods (p > 0.05). The type of sample also had no significant effect on phylogenetic and non-phylogenetic beta diversity metrics (Unifrac and Bray-Curtis; p > 0.05). Conclusion We observed no significant difference between the microbiota of mucosal tissue and swab samples. This suggests that less invasive swab samples are representative of the sinonasal mucosa microbiome and can be used for future sinonasal microbiome studies.


Sample collection
All samples were collected intra-operatively. Endoscopically-guided swabs were taken from the anterior ethmoid sinus. We used flocked swabs (Copan Italia S.p.A., Brescia, Italy) to maximize bacterial yield. After swabbing, mucosal tissue biopsies were taken from the corresponding region in the anterior ethmoid sinuses using a standard-sized Blakesley forceps (Size 3, cup size 5×15mm; Karl Storz, Tuttlingen, Germany). Samples were always taken from this region, irrespective of the presence or absence of pus. To avoid inadvertent contamination, any swabs that could have come into contact with the nasal vestibule during sampling were discarded. The swab head was then immediately separated into a sterile container, placed on ice, and then transported to the laboratory for storage at -80°c.

DNA extraction
Swab heads were thawed, cut into small pieces and then placed in 180 μl of enzymatic lysis buffer (Qiagen, CA, USA) overnight at room temperature. Bead homogenization was then performed (5mm steel beads agitated for 20 seconds at 15 Hz, followed by 0.1mm glass beads for 5 minutes at 30 Hz). The same extraction procedure was carried out for the tissue biopsy samples. The remainder of the extraction protocol was performed as per the Qiagen DNeasy Blood & Tissue Kit instructions (Qiagen, CA, USA). Extracted DNA was stored at -80°c until sequencing.

PCR amplification of 16S rRNA gene and pyrosequencing
Tag-encoded FLX-Titanium amplicon pyrosequencing for bacterial organisms was performed as previously described. [20] Briefly, a selective panbacterial 16S-rRNA gene-directed primer set ("27Fmod" AGRGTTTGATCMTGGCTCAG; and "519Rmodbio" GWATTACCGCGGCK GCTG) was applied against the 16S rRNA gene for PCR amplification. PCR and pyrosequencing was performed by MR DNA Lab (Shallowater, TX). A total of 28 cycles of PCR were performed. DNA was normalized at 7ng/μl average. Sequencing was performed on a Roche 454 sequencer (F. Hoffmann-La Roche Ltd, Basel, Switzerland). All samples included in this study were sequenced in one run.

Bioinformatics pipeline
Raw pyrosequencing data files were then processed using the open-source software pipeline "Quantitative Insights into Microbial Ecology" (QIIME) version 1.8, [21] utilizing virtual machines on the Australian National eResearch Collaboration Tools and Resources (NeCTAR) research cloud. Sequences were trimmed of barcodes and primers and sequence quality control was performed using QIIME's default script settings (sequence length 200-1000 nucleotides; minimum average qualilty score 25; maximum length of homopolymer runs 6; maximum number of ambiguous bases 6). Sequences were then denoised using QIIME's built-in denoiser, set on the Titanium profile. [22] Operational taxonomic units (OTUs) were clustered (openreference OTU "picking" against the August 2013 reference Greengenes 16S rRNA database [23]) at 97% similarity using uclust, [24] and then singleton sequences were removed. Taxonomic assignment of OTUs was then performed using BLAST [25] against the Greengenes 16S rRNA database. [23] After examining read counts, rarefactions of the OTU table were performed to a chosen maximum subsampling depth of 1100 sequences and rarefaction curves were plotted. Summary of taxonomic assignments were plotted as bar charts using QIIME. Observed species richness, as well as phylogenetic and non-phylogenetic alpha diversity metrics (Faith's Phylogenetic Diversity index "PD_Whole_Tree" and Shannon's index, respectively) were recorded and compared at the 1100 rarefaction depth. Phylogenetic and non-phylogenetic beta diversity matrices (Weighted/Unweighted UniFrac, and Bray-Curtis, respectively) were calculated. Three-dimensional Principal Coordinate Analysis (PCoA) plots were generated using EMPeror software [26] bundled with QIIME. Using the principal coordinates of the PCoA plots, a Procrustes transformation was performed (over the first three principal components) of the swab samples against those of the tissue samples. This was done using the QIIME script transform_coordinate_matrices.py. Using this script, an m 2 value is calculated, and Monte-Carlo simulations (1000 permutations) were done to calculate a p-value.

Statistical analysis
All statistics were performed using R statistical software [27] (R Foundation for Statistical Computing, Vienna, Austria) through the IPython notebook interface. [28] Statistical significance was considered at the 0.05 level. Alpha diversity metrics were compared between the two sample types using Wilcoxon signed rank test. Beta diversity distances within-group were compared to between-group non-parametrically using a 999-permutations t-test. Testing for the presence of a significant effect of sample type on beta diversity metrics was also done using permutational multivariate analysis of variance [29] (PERMANOVA), through the function "adonis" present in the vegan R package. [30] To accommodate the paired design, the adonis function was employed using the strata parameter; this allowed the permutations to be done only within the patient variable, not across. We then investigated statistically significant differential relative abundance (MRA) of any bacterial species (of more than 1%) between the two groups using Wilcoxon signed rank tests. Correlation between the taxonomic assignment summaries of the two groups (Pearson's Correlation coefficient) was calculated non-parametrically using a twosided 999 permutations test, using the QIIME script compare_taxa_summaries.py. This script (in paired mode) compares pairs of samples across two groups (within patient) and additionally calculates a summarized value for the whole comparison.

Demographic and clinical data
In total, six CRS patients undergoing endoscopic sinus surgery were included in this study. Two of the patients had CRS without nasal polyposis (CRSsNP) and four had nasal polyposis (CRSwNP). Five patients had concomitant asthma. Two patients had aspirin sensitivity.

Taxonomic summary
The mean number of sequences per sample was 2865.333 sequences (SD 2805.204). Our final OTU table contained 1169 unique OTUs (at 97% similarity) in 12 samples (6 tissues, 6 swabs). Upon taxonomy assignment, these OTUs represented 312 unique bacterial genera across 24 phyla. Fig 1 shows the distribution of bacterial phyla in the studied samples. Only 10 genera, out of the 312, had a mean relative abundance of more than 1%, and these are listed in Table 1.

Effect of sample type on observed species richness and alpha diversity
Rarefactions were performed to a depth of 1100 reads. At the maximum subsampling depth of 1100 reads, the mean observed species richness was 33.3 species (36 species for mucosal tissue, versus 30.6 species for swab). Rarefaction curves for richness in all 12 samples plateau at the maximum depth, (Fig 2A) indicating an adequate sampling procedure. There was no significant difference in the observed species richness between both groups (Wilcoxon Signed Rank test, p = 0.44). There was no significant difference in the phylogenetic and non-phylogenetic

Effect of sample type on beta diversity
Beta diversity distance matrices were generated, using Weighted/unweighted UniFrac and Bray-Curtis distances between all samples. The mean distances between samples within the same sample type group (i.e. within-swab, as well as within-tissue) did not differ significantly from mean distances between samples across sample type groups (i.e. between swab and mucosa).
(Weighted UniFrac, unweighted UniFrac, Bray-Curtis; p > 0.05; Fig 3) This indicates that samples within each group were as similar to each other as to samples across the two groups. Moreover, the mean Bray-Curtis and Weighted UniFrac within-patient distances were significantly lower than between-patient distances (p < 0.05; Fig 3), indicating lower diversity and closer similarity within each patients paired samples (versus between-patient). Similarly, the withinpatient unweighted UniFrac distances were also lower than the between-patient distances but this was not statistically significant. (Fig 3) We then examined the effect of sample type on the Weighted and unweighted UniFrac (phylogenetic) and Bray-Curtis (non-phylogenetic) distance matrices using PERMANOVA, while constraining permutations within the patient variable to account for the paired design. We found no evidence of significant impact of sample type on the Unweighted UniFrac distance matrix (pseudo-F-statistic = 0.93, p = 0.35), Weighted UniFrac

Principal Coordinate Analysis (PCoA) and Procrustes Analysis
We then examined beta diversity distances between samples using Principal Coordinate Analysis (PCoA). The PCoA plots show good clustering of the pair of samples collected from the same patient, but shows less defined clustering of swab versus tissue samples. (Figs 4 and 5).
Procrustes analysis was then done on these PCoA plots, such that to compare the principal coordinates of the swab samples to those of the tissue samples. The m 2 statistic produced by the Procrustes analysis is a value that can range from 0 (in this case the matrices are identical/ Boxplots showing distances from Bray-Curtis, unweighted UniFrac, and Weighted UniFrac distance matrices. Boxplots show no significant difference between within-sampletype and between-sampletype distances, and a statistically significant difference between within-patient and betweenpatient Bray-Curtis and Weighted UniFrac distances. Whiskers extend to cover the whole range of values. * = p < 0.05 (two-tailed Student's t-test, computed through QIIME's script make_distance_boxplots.py).

Investigating significant differential abundance of species between both sampling methods
We then used Wilcoxon ranked sum tests to compare mean relative abundances (MRA) of bacterial genera in the two groups. We only found one statistically significant result   (Table 1) No other genera showed a statistically significant differential abundance between the two sampling methods. (p > 0.05) The largest difference in MRA between swab and tissue was found in Corynebacterium (higher abundance in swabs, 25.6% difference in MRA; Table 1) Pearson's correlation coefficient between the two sampling methods Taxa summaries in both types of sample were compared in each patient in a paired method, using Pearson's correlation coefficient. This was calculated using the QIIME script compare_ taxa_summaries.py. The summarized overall coefficient calculated by the script was 0.77 [95% CI:0.75 to 0.78] (p = 0.001), indicating a strong correlation of taxonomic assignment between both sampling methods. The detailed within patient paired correlation results are found in Table 2, which shows a strong correlation between all pairs except for one pair of samples which showed no or weak correlation.

Discussion
In this study, we investigated whether there was a significant difference in microbiome results between swab and mucosal tissue samples obtained from the same site of the same patients. We investigated this through alpha diversity, beta diversity, PCoA, as well as the taxonomy assignments. Our results suggest a good correlation between swabs and mucosal tissue samples, with no major significant differences in bacterial composition. We found that mucosal tissue biopsy samples had higher observed species richness (30.6 species for swab, versus 36 species for mucosa) as well as alpha diversity indices (Faith's and Shannon's indices) (Fig 2B, 2C and 2D). However, this difference in richness and diversity did not reach statistical significance between the two groups. The plots in Fig 2A show that rarefaction curves have plateaued and reached asymptote. This indicates good sampling at the chosen cutoff depth of 1100 reads, and that a higher depth of sequencing is unlikely to discover more unique species (and increased diversity) for the sampled environment. We then explored whether the type of sample would have an effect on beta diversity metrics (Bray-Curtis and Weighted and unweighted UniFrac), and this showed no significant effect of sample type on beta diversity. (Fig 3) PCoA plots suggested a better clustering of samples by patient, rather than the type of sample. (Figs 4 and 5) Moreover, we calculated a high Pearson's correlation coefficient between the taxa summaries of both sample types (r = 0.77, p = 0.001).
Previous studies showed considerable variability in sinus microbiota between individuals. We believe comparing the microbiota of tissue biopsy and swab samples obtained from the same patient controls for this inter-patient variability. Both sample types were taken from the same sites. This paired comparison design is one of the strengths of this study, as it increases its statistical power and minimizes the effect of confounding variables.
The alpha and beta diversity metrics (Faith's PD_Whole_Tree, Shannon's, Bray-Curtis, Weighted UniFrac) measured in this study were particularly chosen for the following reasons. First, these measures include both non-phylogenetic metrics (Shannon's for alpha diversity and Bray-Curtis for beta diversity) as well as phylogenetic metrics (Faith's PD_Whole_Tree for alpha diversity and UniFrac for beta diversity). Phylogenetic metrics include information about the evolutionary distance between taxa, since they account for the structure of the phylogenetic tree. Non-phylogenetic metrics do not account for this information and only depend on taxa abundance (or prevalence) and evenness. Second, we have also included metrics that account for the relative abundance of taxa (Bray-Curtis, Weighted UniFrac), contrary to using only metrics such as Jaccard's or unweighted UniFrac that depend on absence/presence (i.e. prevalence) of taxa, without accounting for their relative abundance. In this way, we could demonstrate that both methods of sampling (tissue and swab) are not only similar in regards to absence/presence of taxa, but also in regards to their relative abundance within each sample. Despite mucosal tissue samples having higher richness and within-community (α) diversity, we could not demonstrate that this was significantly higher than swab samples. Some analyses in our study approached statistical significance (for example, Shannon's index comparison p value was 0.09 and Weighted UniFrac Permanova p = 0.06), which may imply the presence of a significant difference between swab and tissue samples, especially given the small number of samples in our cohort. Nevertheless, we believe positive results of all analysis types in the current study, which includes other metrics (alpha diversity metrics, beta diversity distances, PCoA plots, Pearson's correlation), supports a general conclusion of tissue and swab similarity. This raises the question of power in the realm of microbiota comparison studies, which is a subject of ongoing research. [31,32] Moreover, what is the clinically-significant "effect size" calculated out of the commonly used diversity measures, and that would correlate with sinus health and/or disease? These questions need to be better defined for future research. Future studies confirming the findings of our study may benefit from including a larger cohort. Another caveat with our current study design is the possibility that taking the swab samples (first) somehow perturbed the surface microbiome of the tissue, which was biopsied from the same site after swabbing (middle meatus/anterior ethmoid region). Although this limitation seems unavoidable with the current study design, we believe this issue would not cause a significant change to the tissue specimens, since the tissue biopsied is typically larger than the surface area covered by the swab head, and thus is not limited only to the thin strip touched by the swab.
Although our rarefaction curves indicate a saturation of diversity at our chosen cut-off rarefaction depth (Fig 2C and 2D), future studies should investigate using the Illumina sequencing platform, which allows sequencing at ultra-high depths, with millions of reads generated per run. [33] The Illumina platform is thus better poised to examine the population belonging to the "rare biosphere", and thus the only limitation to accurately characterizing this rare population would be sequencing errors or limitations of the bioinformatics analysis such as imperfect OTU clustering. [34,35] In our study, it appears that the UniFrac metric (which does not take abundance of the taxa into account, giving abundant and rare taxa equal weight) is less able to show the within-patient similarity of swab and tissue samples, in contrast to Bray-Curtis or Weighted UniFrac (Fig 3). This may suggest either a less than perfect ability of swab samples to accurately characterize the rare biosphere, or is an artefact of the OTU clustering. This is an additional reason to confirm the findings of this experiment on the Illumina platform in the future. Some sinonasal microbiome studies also used lavage sampling. [2,6] Unfortunately, this type of sampling was not investigated in our study.
Although the aim of this study is not to describe the CRS sinonasal microbiome, we also report on the taxa discovered (Table 1). Our findings show that Corynebacteria, Staphylococci and Propionibacteria are the most abundant micro-organisms. (Table 1) These genera were also reported in previous studies as abundant in the sinuses. [4,36]. Interestingly, out of the 312 bacterial genera assigned in our study, only 10 (about 3%) had a mean relative abundance of more than 1% (Table 1). This suggests the presence of a "rare biosphere" in the sinuses. This rare biosphere may constitute then about 97% of all species present in the sinuses. However the role played by these rare taxa, and whether they contribute significantly to sinus health or disease, is still unknown; although other studies suggested they have great functional importance at other sites such as the gut and oral microbiota. [37,38] With this said, the small sample size and lack of a healthy control group makes it difficult to draw further conclusions on the sinus microbiome in CRS patients.

Conclusion
In summary, our results suggest that there is no significant difference between mucosal tissue and swab samples and both methods showed strong correlation. We therefore propose that swab samples are sufficiently representative of the sinonasal mucosa microbiome and therefore can be used for future sinonasal microbiome studies. This obviates the need for invasive mucosal biopsies and also means that sinus microbiome swabs can be obtained from healthy and diseased patients intra-operatively, as well as post-operatively in the clinic. Future studies confirming these findings should explore: investigating sinus lavage samples, including a higher number of individuals, and using sequencing platforms that allow ultra-high depths of sequencing.
Supporting Information S1 File. QIIME taxa plots at the level of genus. (ZIP)