Advertisement

Birch pollen allergen immunotherapy reprograms nasal epithelial transcriptome and recovers microbial diversity

Open AccessPublished:February 12, 2019DOI:https://doi.org/10.1016/j.jaci.2019.02.002
      To the Editor:
      Airway epithelial cells are known to have an important role in allergic rhinitis (AR).
      • Poole A.
      • Urbanek C.
      • Eng C.
      • Schageman J.
      • Jacobson S.
      • O'Connor B.P.
      • et al.
      Dissecting childhood asthma with nasal transcriptomics distinguishes subphenotypes of disease.
      • Joenvaara S.
      • Mattila P.
      • Renkonen J.
      • Makitie A.
      • Toppila-Salmi S.
      • Lehtonen M.
      • et al.
      Caveolar transport through nasal epithelium of birch pollen allergen Bet v 1 in allergic patients.
      • Toppila-Salmi S.
      • van Drunen C.M.
      • Fokkens W.J.
      • Golebski K.
      • Mattila P.
      • Joenvaara S.
      • et al.
      Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
      They constitute the first line of defense against inhaled aeroallergens and are active mediators of innate and adaptive immune responses.
      • Toppila-Salmi S.
      • van Drunen C.M.
      • Fokkens W.J.
      • Golebski K.
      • Mattila P.
      • Joenvaara S.
      • et al.
      Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
      Their aberrant functioning is linked with an intake of allergens,
      • Joenvaara S.
      • Mattila P.
      • Renkonen J.
      • Makitie A.
      • Toppila-Salmi S.
      • Lehtonen M.
      • et al.
      Caveolar transport through nasal epithelium of birch pollen allergen Bet v 1 in allergic patients.
      and their transcriptome is reprogrammed under exposure to pollens
      • Joenvaara S.
      • Mattila P.
      • Renkonen J.
      • Makitie A.
      • Toppila-Salmi S.
      • Lehtonen M.
      • et al.
      Caveolar transport through nasal epithelium of birch pollen allergen Bet v 1 in allergic patients.
      • Toppila-Salmi S.
      • van Drunen C.M.
      • Fokkens W.J.
      • Golebski K.
      • Mattila P.
      • Joenvaara S.
      • et al.
      Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
      as well as in AR
      • Toppila-Salmi S.
      • van Drunen C.M.
      • Fokkens W.J.
      • Golebski K.
      • Mattila P.
      • Joenvaara S.
      • et al.
      Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
      and atopic asthma.
      • Poole A.
      • Urbanek C.
      • Eng C.
      • Schageman J.
      • Jacobson S.
      • O'Connor B.P.
      • et al.
      Dissecting childhood asthma with nasal transcriptomics distinguishes subphenotypes of disease.
      Furthermore, epithelial cells interact with and are involved in generating an environmental niche for the respiratory microbiota, whose imbalance has been associated with seasonal AR
      • Choi C.H.
      • Poroyko V.
      • Watanabe S.
      • Jiang D.
      • Lane J.
      • deTineo M.
      • et al.
      Seasonal allergic rhinitis affects sinonasal microbiota.
      and childhood rhinitis and asthma.
      • Chiu C.Y.
      • Chan Y.L.
      • Tsai Y.S.
      • Chen S.A.
      • Wang C.J.
      • Chen K.F.
      • et al.
      Airway microbial diversity is inversely associated with mite-sensitized rhinitis and asthma in early childhood.
      However, the precise functions of epithelial host cells and respiratory microbes in AR are still largely elusive, especially during pollen allergen immunotherapy (AIT) that is associated with symptom reduction,
      • Bousquet J.
      • Khaltaev N.
      • Cruz A.A.
      • Denburg J.
      • Fokkens W.J.
      • Togias A.
      • et al.
      Allergic Rhinitis and its Impact on Asthma (ARIA) 2008 update (in collaboration with the World Health Organization, GA(2)LEN and AllerGen).
      decrease in allergen-specific biomarkers, and altered T- and B-cell responses.
      • Akdis C.A.
      • Akdis M.
      Mechanisms of allergen-specific immunotherapy.
      We collected nasal brushings for RNA sequencing from 5 healthy subjects and 3 birch pollen AR patients with and without AIT at 2 springs and winters and studied seasonal, AR, and AIT-related alterations in the nasal epithelial and microbial transcriptomes (Fig 1, A; see Fig E1 and Table E1 in this article's Online Repository at www.jacionline.org). Pollen count and AR symptom information were also assessed, revealing the presence of high amounts of birch pollen at spring samplings (Fig 1, B) and a marked improvement in quality of life in subjects with AR with AIT compared with controls (P < .005) and subjects with AR without AIT (P < .03) but not between other groups (Fig 1, C).
      Figure thumbnail gr1
      Fig 1Study overview. A, Study flowchart showing number of subjects, sampling points, and the start of the AIT. Samples were collected at 4 consecutive sampling points from 5 healthy control subjects and 6 subjects with AR. Three patients with AR started AIT. All subjects were without medication for at least 4 weeks before sampling. In springs, patients with symptomatic AR were without antihistamines for at least 3 days before sampling. B, Counts of birch and total pollen during the course of the study in gray and black, respectively. Counts of other pollens than birch were under detection level during sampling during the spring samplings. There were no counts of pollen in the air during the winter samplings. C, Total visual analogue scale (VAS) symptom score at the day of sampling. Control and AR-AIT groups (P < .005) as well as AR-AIT and AR-noAIT groups (P < .023) differed in interaction by 2-way repeated-measures ANOVA. Statistically significant interaction was not observed between control and AR-noAIT groups at the alpha level of 0.05.
      RNA sequencing resulted in 90 million mappable reads per sample on average. Of all the annotated human protein-coding genes, 17,347 were deduced expressed and 360 differentially expressed between different time points within groups and between different groups within time points (Fig 2, G). Identified were also 166 (Fig 2, A and B) and 17 (Fig 2, D and E) protein-coding genes with an altered expression between the consecutive springs and winters, respectively. Notably, we identified the greatest transcriptional reprogramming between springs in the AR-AIT group, indicating that AIT alters epithelial expression in the presence of allergens. Analyses also revealed 3 allergy-related pathways that were affected between the spring samplings. An asthma pathway was found to be altered in AR-noAIT subjects, whereas Toll-like receptor (TLR) and chemokine signaling pathways were both affected in AR-noAIT and AR-AIT subjects (Fig 2, C; see Fig E2 in this article's Online Repository at www.jacionline.org). Pathway enrichment analysis of winter data revealed pathways with coordinated expression change only in healthy controls (Fig 2, F). Analysis of expressed variants pinpointed in turn 8 variants expressed in 2 or more subjects with AR at some time point but in none of the healthy controls (see Fig E3 in this article's Online Repository at www.jacionline.org).
      Figure thumbnail gr2
      Fig 2Overview of RNA sequencing data. A and D, Protein-coding genes statistically differentially expressed (Q value ≤ 0.1 and absolute log2 fold-change ≥ 1.5) between spring (Fig 2, A) and winter (Fig 2, D). The heatmap was drawn using log2 (+1 offset) counts per million expression values, mean centered and scaled by gene averages. Red indicates upregulation of the gene and blue downregulation of the gene relative to the average. B and E, Venn diagrams show the total number of differentially expressed genes between spring (Fig 2, B) and winter (Fig 2, E) sampling points at each group. C and F, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched among differentially expressed genes per group. Shades of blue and red indicate significance of the enrichment, and size of the dot represents gene count. Listed at the bottom in parentheses is the total number of differentially expressed genes in each group with an association to some KEGG pathway. G, Number of differentially expressed genes in selected-pairwise comparisons. Comparisons not shown are as follows: AR-noAIT Nov2 vs AR-noAIT May1 (52), AR-AIT Nov2 vs AR-AIT May1 (1), AR-AIT May1 vs AR-noAIT May1 (7), AR-AIT Nov1 vs AR-noAIT Nov1 (0), AR-AIT May2 vs AR-noAIT May2 (6), AR-AIT Nov2 vs AR-noAIT Nov2 (1), Control May2 vs Control May1 (27), Control Nov2 vs Control May1 (7), and Control Nov2 vs Control Nov1 (17). H, Relative abundances of microbial (archaeal, bacterial, and viral) genera and average microbial load per group. Only genera accounting for more than 5% of the total microbial load in any group are shown. The line denotes the average number of microbe-classified reads within the group.
      Further analysis of the gene expression profiles of the 3 allergy pathways between the spring samplings highlighted marked similarities in the AR-AIT and control groups that were not seen in the AR-noAIT group (Fig E2). These results imply that AIT may restore epithelial gene expression toward normal and indicate that effectivity of AIT could be screened from nasal epithelium in addition to leukocytes. Specifically, the MHCII components were upregulated at the second spring in AR-AIT and control groups but not in the AR-noAIT group (Fig E2, A), indicating that AIT restores the compromised antigen-presenting capacity of epithelial cells in AR. We also found that genes that are downstream effectors of the chemokine signaling or pattern recognition and provide proinflammatory, antiviral, chemotactic, and T-cell stimulatory effects behaved alike between the AR-AIT and control groups (Fig E2, B and C). These findings are in line with the findings that changes in expression of TLR genes are associated with AR and suggest a role for TLR agonists in treatment of AR.
      • Toppila-Salmi S.
      • van Drunen C.M.
      • Fokkens W.J.
      • Golebski K.
      • Mattila P.
      • Joenvaara S.
      • et al.
      Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
      • Akdis C.A.
      • Akdis M.
      Mechanisms of allergen-specific immunotherapy.
      Notably, the expression of several asthma-related genes was found to be in opposite between the AR-AIT and AR-noAIT subjects (Fig E2).
      Microbial classification of sequencing data was performed to explore whether AR alters nasal microbiota (archaeal, bacterial, and viral) and whether AIT could restore microbial imbalances toward normal. On average, approximately 500 counts per million (∼16,340 read-pairs) per sample were assigned to microbial taxa, 98.13% of which received a genus-level classification (see Fig E4, A, in this article's Online Repository at www.jacionline.org). The classification showed that bacteria, archaea, and viruses were part of the active nasal microbiota, the most common genera being Bacillus (average abundance, 42.23%), Methanocaldococcus (average abundance, 35.72%), and Alpharetrovirus (average abundance, 4.32%). Similar to previous studies,
      • Lal D.
      • Keim P.
      • Delisle J.
      • Barker B.
      • Rank M.A.
      • Chia N.
      • et al.
      Mapping and comparing bacterial microbiota in the sinonasal cavity of healthy, allergic rhinitis, and chronic rhinosinusitis subjects.
      a large sample-to-sample variation was observed (Fig E4, A). Particularly, 6 samples taken at the second spring varied greatly from the rest (Fig E4, A and B) and were, for instance, the drivers of the greater abundance of viruses at the second spring compared with the other time points (Fig 2, H). Interestingly, examination of changes in species abundancies (Fig E4, C) pinpointed Pseudomonas aeruginosa to be more abundant in the first spring in comparison to the second spring in the AR-AIT group.
      We also computed alpha diversities to evaluate the effect of AR and AIT on the microbial diversity of nasal epithelia (see Fig E5, Fig E6, Fig E7, A-N, in this article's Online Repository at www.jacionline.org). This analysis revealed that control subjects primarily had the highest alpha diversity, differing from that seen previously in a study on seasonal AR
      • Choi C.H.
      • Poroyko V.
      • Watanabe S.
      • Jiang D.
      • Lane J.
      • deTineo M.
      • et al.
      Seasonal allergic rhinitis affects sinonasal microbiota.
      but similar to that focusing on children with asthma and rhinitis.
      • Chiu C.Y.
      • Chan Y.L.
      • Tsai Y.S.
      • Chen S.A.
      • Wang C.J.
      • Chen K.F.
      • et al.
      Airway microbial diversity is inversely associated with mite-sensitized rhinitis and asthma in early childhood.
      Interestingly, most diversity indices suggested an increase in diversity between the first and second winters in all groups. Most prominent was the increase in the AR-AIT group, whereas some increase was also detectable in the control and AR-noAIT groups (Fig E6, A-N). The diversity at the second winter in the AR-AIT group also changed more toward that of the control group than what was the corresponding change in the AR-noAIT group (Fig E6, A-N). These findings are largely in line with the previous studies noting that the bacterial diversity varies during the allergy season
      • Choi C.H.
      • Poroyko V.
      • Watanabe S.
      • Jiang D.
      • Lane J.
      • deTineo M.
      • et al.
      Seasonal allergic rhinitis affects sinonasal microbiota.
      and suggest that AIT may increase microbial diversity and restore it toward normal.
      Limitations of this study include the small subject number, lack of placebo group, differences in baseline allergic symptoms between the groups, differences in pollen seasons, differences in air quality, and technical differences in sampling, which may in part have compromised results. Yet, the study provided interesting insights into the epithelial transcriptome during AIT and revealed that AIT causes subtle but significant alterations in asthma, TLR signaling, and chemokine signaling-related genes and may as well recover microbiological diversity toward normal. Seasonal heterogeneity represented the largest source of variation in transcriptomes, indicating a need for novel biomarkers in AIT treatment monitoring that accommodate inherent heterogeneity and seasonal variation.
      We acknowledge Anne-Maria Konkola, Docent Sanna Korkonen, Docent Pekka Malmberg, Leena Petman, Mirja-Liisa Sipola, BDent Emma Terna, Tanja Utriainen, and Docent Jan Weckström, personnel at the Sequencing and Bioinformatics Units at FIMM Technology Centre supported by the Helsinki Institute of Life Science and Biocenter Finland , and the volunteer subjects and their family members for making this study possible.

      Methods

       Subjects

      Study subjects were recruited from the Skin and Allergy Hospital of Helsinki University Hospital. The study plan was approved by the ethical committee of Hospital District of Helsinki and Uusimaa, Finland (permission no. 19/13/003/00/11). Written informed consent was received from all subjects and their parents if the age of the participant was less than 18 years. The study has been registered in ClinicalTrials.com (no. NCT01985542). Baseline data of the study subjects are presented in Table E1. The total number of participants entering the study was 23 (Fig E1). The AR-AIT group received subcutaneous immunotherapy in November 2011 after the second sampling visit (Fig 1 and Fig E1), meaning that 2 samplings of the AR-AIT group were performed before and 2 during AIT.

       Nasal brushings and RNA extraction

      Nasal epithelial brushing was performed to middle meatus of both sides of nasal cavity after slight blowing of nose without local anesthesia as described.
      • Mattila P.
      • Renkonen J.
      • Toppila-Salmi S.
      • Parviainen V.
      • Joenvaara S.
      • Alff-Tuomala S.
      • et al.
      Time-series nasal epithelial transcriptomics during natural pollen exposure in healthy subjects and allergic patients.
      Epithelial cells were collected at 4 time points, washed once with ice-cold nuclease-free PBS, and resuspended immediately into RNAlater RNA stabilization reagent (Qiagen, Hilden, Germany) to preserve RNA profiles. The epithelial RNA isolation was done next day using Qiagen RNeasy Mini Kit with the optional DNAse treatment included.

       Library preparation and RNA sequencing

      Agilent Bioanalyzer RNAnano chip (Agilent, Santa Clara, Calif) was used to evaluate the integrity of RNA and Qubit RNA kit (Life Technologies, Carlsbad, Calif) to quantitate RNA in epithelial cell samples. If acceptable in quality (RIN value >7), 1.0 μg of total RNA sample was ribodepleted and prepared to RNA sequencing library by using ScriptSeq v2 Complete kit (Illumina, Inc, San Diego, Calif). RNA sequencing libraries were purified with SPRI beads (Agencourt AMPure XP, Beckman Coulter, Brea, Calif). The library quality control was evaluated on high-sensitivity chips by Agilent Bioanalyzer (Agilent). Paired-end sequencing of sequencing libraries with 100-bp read length was performed using Illumina HiSeq technology (HiSeq 2000, Illumina, Inc, San Diego, Calif). Planned read amount was 40 million reads per sample.

       RNA sequencing data processing

      RNA sequencing data were preprocessed as described previously.
      • Kumar A.
      • Kankainen M.
      • Parsons A.
      • Kallioniemi O.
      • Mattila P.
      • Heckman C.A.
      The impact of RNA sequence library construction protocols on transcriptomic profiling of leukemia.
      Briefly, Trimmomatics
      • Bolger A.M.
      • Lohse M.
      • Usadel B.
      Trimmomatic: a flexible trimmer for Illumina sequence data.
      was used to correct read data for low quality, Illumina adapters, and short read length. Filtered paired-end reads were aligned to the human genome (GRCh38) using the STAR
      • Dobin A.
      • Davis C.A.
      • Schlesinger F.
      • Drenkow J.
      • Zaleski C.
      • Jha S.
      • et al.
      STAR: ultrafast universal RNA-seq aligner.
      with the guidance of EnsEMBL v82 gene models. Default 2-pass per-sample parameters were used, except that the overhang on each side of the splice junctions was set to 99. The alignments were then sorted and PCR duplicates were marked using Picard, feature counts were computed using SubRead,
      • Liao Y.
      • Smyth G.K.
      • Shi W.
      The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote.
      feature counts were converted to expression estimates using Trimmed Mean of M-values normalization,
      • Robinson M.D.
      • Oshlack A.
      A scaling normalization method for differential expression analysis of RNA-seq data.
      and lowly expressed genomic features with counts per million (CPM) value less than or equal to 1.00 in less than half the controls or birch pollen patients were removed. Default parameters were used, with the exception that reads were allowed to be assigned to overlapping genome features in the feature counting.

       Host gene expression analysis

      Differential expression testing was performed using the edgeR
      • Robinson M.D.
      • McCarthy D.J.
      • Smyth G.K.
      edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.
      software and included testing of differential expression between and within groups at different sampling points. In the statistical testing, comparisons between subject groups used a combined factor of subject group and sampling point, whereas comparisons within subject groups also used a factor for the subject. The resulting P values were adjusted Storey's Q-value approach, with significance defined as a Q value of less than or equal to 0.10. A cutoff value of absolute log2 fold-change of greater than or equal to 1.5 and EnsEMBL v82 biotype annotations were used as additional filters to select differentially expressed genes (DEGs) with protein-coding annotation for the downstream analysis. Heatmaps of differentially expressed protein-coding genes were produced with pheatmap R package.
      • Kolde R.
      pheatmap: Pretty Heatmaps. R package version 1.0.10.
      Hierarchical clusters were generated using the spearman correlation and ward.D2 as the linkage method, with the exception of using ward.D2 and Euclidean distance for genes that were differentially expressed between different sampling years at springs and using complete linkage and spearman correlation for genes that were differentially expressed between different sampling years at winters. CPM data were used to generate heatmaps. Venn diagrams were generated using the VennDiagram R package.
      • Chen H.
      VennDiagram: generate high-resolution Venn and Euler plots. R package version 1.6.20.
      Functional profiles of DEGs were investigated with clusterProfiler
      • Yu G.
      • Wang L.G.
      • Han Y.
      • He Q.Y.
      clusterProfiler: an R package for comparing biological themes among gene clusters.
      using functions enrichGO and enrichKEGG. Outputs of enrichment analyses were visualized using dotplot function in clusterProfiler. Biologically relevant pathways found by clusterProfiler were visualized using pathview R package.
      • Luo W.
      • Brouwer C.
      Pathview: an R/Bioconductor package for pathway-based data integration and visualization.
      In the process, Kyoto Encyclopedia of Genes and Genomes (KEGG) gene IDs of the selected pathways were fed along with log2 fold-change values from relevant comparisons. Color codes on the pathway map were used to illustrate genes that were differentially expressed and the direction of their expression changes. Fold-change values beyond that range were truncated to the closest extreme; that is, values greater than 2 were truncated to 2, and values less than −2 truncated to −2. Downstream analyses were performed using R 3.3.1 with Bioconductor 3.0.

       Variant analysis

      Transcript variants were called from STAR alignments using the GATK best practices workflows for transcriptome data
      • McKenna A.
      • Hanna M.
      • Banks E.
      • Sivachenko A.
      • Cibulskis K.
      • Kernytsky A.
      • et al.
      The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.
      and then annotated using Annovar
      • Wang K.
      • Li M.
      • Hakonarson H.
      ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.
      as defined previously.
      • Kumar A.
      • Kankainen M.
      • Parsons A.
      • Kallioniemi O.
      • Mattila P.
      • Heckman C.A.
      The impact of RNA sequence library construction protocols on transcriptomic profiling of leukemia.
      Quality control analyses were performed as defined previously.
      • Kumar A.
      • Kankainen M.
      • Parsons A.
      • Kallioniemi O.
      • Mattila P.
      • Heckman C.A.
      The impact of RNA sequence library construction protocols on transcriptomic profiling of leukemia.
      Variant calls were further filtered by accepting only those that were present in 2 or more AR cases, not present in any control case, and predicted to be pathogenic by various pathogen prediction methods part of the Annovar
      • Wang K.
      • Li M.
      • Hakonarson H.
      ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.
      annotation tool. Heatmap was plotted using pheatmap. The functional effects of variants were taken from Annovar
      • Wang K.
      • Li M.
      • Hakonarson H.
      ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.
      outputs. In addition, we plotted barplot using CPM expression value of genes in healthy control and AR groups.

       Microbial community profiling

      Microbial community profiling was performed as previously described
      • Dufva O.
      • Kankainen M.
      • Kelkka T.
      • Sekiguchi N.
      • Awad S.A.
      • Eldfors S.
      • et al.
      Aggressive natural killer-cell leukemia mutational landscape and drug profiling highlight JAK-STAT signaling as therapeutic target.
      with some modifications. Specifically, RNA sequencing data were preprocessed for adapter trimming, low-quality bases filtering, and removal of reads less than 36 bp in length by using Trimmomatic.
      • Bolger A.M.
      • Lohse M.
      • Usadel B.
      Trimmomatic: a flexible trimmer for Illumina sequence data.
      Paired-end reads passing the preprocessing were mapped against rRNA sequences from RFAM
      • Nawrocki E.P.
      • Burge S.W.
      • Bateman A.
      • Daub J.
      • Eberhardt R.Y.
      • Eddy S.R.
      • et al.
      Rfam 12.0: updates to the RNA families database.
      v12.3 using the Burrows-Wheeler Aligner
      • Li H.
      • Durbin R.
      Fast and accurate short read alignment with Burrows-Wheeler transform.
      with default settings, and reads-matching rRNAs were filtered by using samtools.
      • Li H.
      • Handsaker B.
      • Wysoker A.
      • Fennell T.
      • Ruan J.
      • Homer N.
      • et al.
      The Sequence Alignment/Map format and SAMtools.
      Centrifuge
      • Kim D.
      • Song L.
      • Breitwieser F.P.
      • Salzberg S.L.
      Centrifuge: rapid and sensitive classification of metagenomic sequences.
      was then used to classify paired-end reads to microbial taxa. Alignment data were converted to kraken-style output. In the classification, reads were aligned against 27,127 known complete bacterial, archaeal, and viral genome assemblies, the human genome, and 10,615 technical artifact sequences that were available in the RefSeq
      • O'Leary N.A.
      • Wright M.W.
      • Brister J.R.
      • Ciufo S.
      • Haddad D.
      • McVeigh R.
      • et al.
      Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation.
      database at February 2018. Default parameters were used, with the exception that only 1 (ie, the lowest common ancestor) assignment was reported for read-pairs with multiple primary assignments. Taxa having less than 100 read-pairs assigned to them in any sample were removed. Pairwise comparisons between and within groups at different sampling points were performed by applying DeSeq2
      • Love M.I.
      • Huber W.
      • Anders S.
      Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.
      on the number of reads covered by the clade rooted at the given taxon level. In the analyses, size factors were estimated by using the poscounts method, comparisons between subject groups were done using a combined factor of subject group and sampling point, and comparisons within subject groups with a model where individuals were nested within subject groups. Each taxonomic level was analyzed separately and variance stabilizing transformation was used to generate expression estimates for heatmap visualizations. The Storey's Q-value adjustment
      • Storey J.D.
      The positive false discovery rate: a Bayesian interpretation and the q-value.
      was used to correct data for multiple hypothesis testing, with significance defined as a Q value of less than or equal to 0.05. Finally, alpha diversity (Observed, Chao1, ACE, Shannon, Simpson, InvSimpson, and Fisher), beta diversity (Bray-Curtis dissimilarity), and rarefaction analyses were done using the Phyloseq software
      • McMurdie P.J.
      • Holmes S.
      phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data.
      applied on the number of reads assigned directly to the given taxonomic level.

      Results

       General

      Subjects with AR, and especially AR-AIT cases, had higher median of total serum IgE, median of birch-specific serum IgE and skin prick test wheel diameter to birch and symptom scores during samplings (Table E1 and Fig 1). The AR-AIT group reported benefit (Fig 1) and reported no severe side effects at the end of subcutaneous immunotherapy 3 years after start of subcutaneous immunotherapy (data not shown).

       Transcriptome of nasal epithelium

      We generated in total 5164 million raw paired-end transcript reads. Manual inspection of the quality plots generated using FASTQ indicated that sequencing data were of excellent quality. On average, 90 million reads were mapped to the reference per subject. Mapped reads were then used to generate CPM expression estimates, revealing expression of 13,873 protein-coding genes among healthy controls and 17,347 across all the 44 samples. Altogether, expression of 34,896 genes was found. To gain insight into cellular processes dysregulated in AR and AIT, DEGs between sample groups were identified using edgeR.
      • Akdis C.A.
      • Akdis M.
      Mechanisms of allergen-specific immunotherapy.
      In this analysis, we identified altogether 360 genes to be differentially expressed with the Q value less than or equal to 0.1 and absolute log2 fold-change greater than or equal to 1.5 between different time points within group and between different groups within time points.

       Alterations in gene expression profiles in the 2 following springs and winters

      Comparison of the transcriptome profiles between the springs revealed 119 DEGs between the samplings in the AR-AIT group, 49 between the samplings in the AR-noAIT group, and 27 between the samplings in healthy controls (Fig 2, B), which suggests that the greatest transcriptional reprogramming took place in the AR-AIT group followed by the AR-noAIT group and healthy subjects. Comparison of the 2 consecutive winters revealed only 17 DEGs among healthy controls and none among patients with AR, suggesting that AIT alters epithelial expression only in the presence of allergens (Fig 2, E).

       Differential expression of immune response and signaling pathways

      We performed KEGG pathway enrichment analysis to discover functional themes shared by DEGs. This analysis revealed altogether 21 KEGG pathways with coordinated expression change between the spring samplings (Fig 2, C). Of these, 4 were associated with genes differentially expressed in the AR-AIT group, including the chemokine signaling and TLR signaling pathway (Fig 2, C). Genes differentially expressed in the AR-noAIT group were in turn associated with 8 pathways, including IL-17 signaling and asthma pathway that were discovered only in this comparison (Fig 2, C). The healthy group genes were enriched in 11 KEGG pathways (Fig 2, C). Altogether we detected 3 allergy-related pathways, of which asthma was discovered only in the AR-noAIT comparison and TLR signaling and chemokine signaling pathways were discovered in AR-noAIT and AR-AIT comparisons (Fig 2, C). Pathway enrichment analysis of winter comparison data revealed pathways with coordinated expression change only in healthy controls (Fig 2, F).

       In-depth analysis of allergy-related pathways

      Allergy-related pathways found in the pathway analysis (Fig E2) were analyzed more in-depth to study these mechanisms. First, the asthma pathway consisted in total 3 DEGs. The AR-AIT group and healthy controls displayed upregulation of MHCII and downregulation of FcεRI at the second spring (Fig E2, A). The expression of various other members of the pathway was altered, although unsignificantly, between spring samplings (Fig E2, A). These more borderline findings included IL-13, which is a T-cell–specific transcription factor, and interleukin
      • Pawankar R.
      • Mori S.
      • Ozu C.
      • Kimura S.
      Overview on the pathomechanisms of allergic rhinitis.
      IL-4, which is IgE synthesis switch factor,
      • Pawankar R.
      • Mori S.
      • Ozu C.
      • Kimura S.
      Overview on the pathomechanisms of allergic rhinitis.
      and IL-5, which is an eosinophil growth factor.
      • Pawankar R.
      • Mori S.
      • Ozu C.
      • Kimura S.
      Overview on the pathomechanisms of allergic rhinitis.
      Expression of genes of the asthma pathway between time points occurred to opposite directions in AR-AIT and AR-noAIT groups (Fig E2, A). The second pathway of interest was chemokine signaling pathway. This pathway consisted of 49 dysregulated genes (Fig E2, B), including various chemokines and chemokine receptors. Expression of genes between time points occurred to opposite directions in AR-AIT and AR-noAIT groups (Fig E2, B). The third pathway of relevance to AR was the the TLR signaling pathway (Fig E2, C) consisting altogether of 67 dysregulated genes, such as p38, TNF-α, IL-12, and INF-α genes. Expression of genes between springs happened to opposite directions in AR-AIT and AR-noAIT groups (Fig E2, C ).

       Transcript variants expressed in patients with AR

      The GATK best practice for RNA sequencing
      • McKenna A.
      • Hanna M.
      • Banks E.
      • Sivachenko A.
      • Cibulskis K.
      • Kernytsky A.
      • et al.
      The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.
      was used to identify expressed variants. This approach revealed altogether 3,268,177 (on average 74,277 per subject) variants passing GATK filters. Removal of intronic and silent variants and polymorphisms resulted in 8174 (on average 186 per subject) variants that were further narrowed down to 8 potential candidates expressed in at least 2 subjects with AR at any time point but in none of the healthy control samples (Fig E3).

       Functional characterization of microbiome showed no AR-related changes

      The effect of AR and AIT on the active nasal microbiota was studied by identifying microbial reads from the RNA sequencing data and performing microbial classification. On average approximately 16,340 read-pairs (500 CPMs) per sample were assigned to bacterial, archaeal, or viral taxa with a high sample-to-sample variation (minimum of 791 reads [24 CPMs] and maximum of 69,428 reads [1,791 CPMs]). Of these microbial reads, on average approximately 98.13% were classified at the genus level and 67.08% at the species level. The genus-level classifications (Fig 2, H, and Fig E4) showed that overall the most abundant genera were Bacillus, Methanocaldococcus, and Alpharetrovirus, with average relative proportions of approximately 42.23%, approximately 35.72%, and approximately 4.32%, respectively, and other genera having average relative proportions of approximately 1.57% (Acinetobacter) or less. The relative proportions of the detected genera varied greatly between the samples, with Bacillus demonstrating the greatest variation (relative proportions ranging from 0% to ∼63.76%). The relative proportion of viruses, particularly Alpharetrovirus, was greater in the second spring sampling point than in the other sampling points for all the groups, whereas the relative proportion of Bacillus was reduced (Fig 2, H, and Fig E4, A). This variation, however, was mainly driven by the second spring samples of 6 cases, who were distributed among all groups (Fig E4, A). Examination of the community compositional differences between the samples revealed these 6 samples to be the most dissimilar from the rest, which in turn had rather similar community compositions (Fig E4, B). To further investigate microbial richness, we estimated the alpha diversity measures for all the samples (Fig E5, Fig E6, Fig E7, A-N). Most indices demonstrated an increase in the diversity for the AR-AIT group when the second winter was compared with the first (Fig E6, A-N). Some change in diversity between second and first winters was also seen for the control group, whereas some indices also indicated an increase for the AR-noAIT group. The increase seen in the AR-noAIT group was, however, less than that observed for the AR-AIT group. The second winter diversity indices of the AR-AIT group were also mainly closer to those of the control group than were the indices of the AR-noAIT group (Fig E6, A-N). Finally, examination of significant (Q value ≤ 0.05) changes in species abundancies (Fig E4, C) revealed more changes in all groups between the 2 spring sampling points than between the winter sampling points. The winter comparisons revealed only increased species for the control and AR-AIT groups, whereas for the AR-noAIT group only decreased species were reported. The AR-AIT group had significantly less Pseudomonas aeruginosa in the second spring sampling point than the first. For the other groups, the comparisons revealed no significant changes in the abundance of this bacterium. >
      Figure thumbnail fx1
      Fig E1Flowchart of the study. The total number of participants entering the study was 23. Initially, 7 patients with AR were assigned to pollen AIT (AR-AIT group) and 9 others to conventional therapy (AR-noAIT group). The study also included 7 control subjects (control group) without allergy. One participant in the AR-AIT group discontinued the study before starting AIT because of a diagnosis of a cardiac disease requiring surgery. The other participants in the AR-AIT group fulfilled the whole AIT scheme of 3 years, with the total dose being more than 2000.000 standardized quality (SQ). AIT was performed according to standard 3-year protocol or normal 3-year scheme.
      • Arvidsson M.B.
      • Lowhagen O.
      • Rak S.
      Effect of 2-year placebo-controlled immunotherapy on airway symptoms and medication in patients with birch pollen allergy.
      Subcutaneous injections of birch pollen extract (Betula verrucosa, ALK-Abelló) were administered in the clinic and included an induction phase with increasing dosing starting with a dose of 20 SQ. The maintenance-phase dose was 100,000 SQ. Six subjects discontinued the study before the last follow-up. Moreover, 5 of the 16 subjects who completed the study had poor RNA quality in at least 1 nasal epithelial sample, leading to their exclusion. Thus, a total of 11 subjects (5 healthy controls, 3 AR-AIT, and 3 AR-noAIT) completed the study. These participants were adults and of white European ancestry except 1 Chinese male in the AR-noAIT group. NGS, Next Generation Sequencing; SPT, skin prick test.
      Figure thumbnail fx2
      Fig E2Allergy-related pathways identified in the pathway analysis. Gene expression profiles of (A) asthma, (B) chemokine signaling, and (C) TLR signaling pathways among patients with AR with pollen AIT, patients with AR without AIT, and control subjects. Boxes in the figure represent genes. The gradient colors indicate the log2 fold-change of the gene between the spring samplings in AR-AIT (I/left), AR-noAIT (II/middle), and control (III/right) groups. Fold-change values beyond the range (from −2 to 2) were truncated to the closest extreme; that is, values greater than 2 were truncated to 2, and values less than −2 truncated to −2. Asterisk indicates statistically significant difference with a Q value of less than or equal to 0.10. Green gradient colors indicate upregulation of the gene at the second spring, and red colors indicate upregulation at the first spring. AC, Adenyl cyclase; AKT, AKT Serine/Threonine Kinase; AP-1, Activator protein 1; APC, Antigen presenting cell; BAD, Bcl-2-antagonist of cell death; BCR, BCR, RhoGEF And GTPase Activating Protein; Ca2+, Calcium ion; CASP8, Caspase 8; CD14, Monocyte differentiation antigen CD14; CD40, Tumor necrosis factor receptor superfamily member 5; CD40L, Tumor necrosis factor ligand superfamily member 5; CD80, CD80 antigen; CD86, CD86 antigen; CDC42, Cell division control protein 42; Crk, CRK Proto-Oncogene, Adaptor Protein; DAG, Diacylglycerol; DC, Dendridic cell; DOCK2, Dedicator of cytokinesis protein 2; ECP, Eosinophil cationic protein; ELMO1, Engulfment and cell motility protein 1; EPO, Eosinophil peroxidase; ERK, Mitogen-activated protein kinase; FADD, FAS-associated death domain protein; FAK, Focal adhesion kinase; FCεR1, High affinity immunoglobulin epsilon receptor; FOXO3, Forkhead box protein O3; Gαi, G i alpha subunit; Gβγ, G beta-gamma complex; GRB2, Growth factor receptor-binding protein 2; GRK, Rhodopsin kinase; GSK3, Glycogen synthase kinase 3; IFNα, Interferon alpha; IFNαβR, Interferon receptor; IFNβ, Interferon beta; IgE, Immunoglobulin E; IkappaBα, NFKB Inhibitor Alpha; IKβ, Inhibitor of nuclear factor kappa-B; IKK, IκB kinase; IKKα, Inhibitor of nuclear factor kappa-B kinase subunit alpha; IKKβ, Inhibitor of nuclear factor kappa-B kinase subunit beta; IKKε, Inhibitor of nuclear factor kappa-B kinase subunit epsilon; IKKγ, Inhibitor of nuclear factor kappa-B kinase subunit gamma; IP10, Interferon Gamma-Induced Protein 10; IP3, Inositol trisphosphate; IRAK1, Interleukin-1 receptor-associated kinase 1; IRF, Interferon regulatory factor; I-TAC, Interferon-inducible T-cell alpha chemoattractant; ITK, IL2-inducible T-cell kinase; Jak, Janus kinase; JNK, Mitogen-activated protein kinase 8/9/10 (c-Jun N-terminal kinase); LBP, Lipopolysaccharide-binding protein; LTC4, Leukotriene C4; LTD4, Leukotriene D4; MBP, Bone marrow proteoglycan; MD2, Lymphocyte antigen 96; MEK1, Mitogen-activated protein kinase kinase 1; MEK1/2, Mitogen-activated protein kinase; MHCII, Major histocompatibility complex, class II; MIG, Monokine induced by gamma interferon; MIP1, MEKK2-interacting protein; MKK, Mitogen-activated protein kinase; MYD88, Myeloid differentiation primary response protein MyD88; NADPH, 47-kilodalton cytosolic subunit of the multi-protein complex known as NADPH oxidase; NFKB, Nuclear factor NF-kappa-B; NFKβ, Nuclear factor NF-kappa-B; NO, Nitric oxide; OPN, Secreted phosphoprotein 1; p105, NF-kappa-B p105 subunit; p130Cas, Cas Family Scaffold Protein; P13K, Phosphoinositide-3-kinase; P38, p38 MAP kinase; p47phox, NADPH oxidase; PAF, Platelet activating factor; PAK1, p21-activated kinase 1; PAR3, Partitioning defective protein 3; PGD2, Prostaglandin D2; PI3K, Phosphoinositide 3-kinase; PIP3, Phosphatidylinositol 3,4,5 trisphosphate; PKA, Protein kinase A; PKC, Protein kinase C; PLCβ, Phosphatidylinositol phospholipase C, beta; PREX1, Phosphatidylinositol-3,4,5-trisphosphate-dependent Rac exchanger 1 protein; PYK2, Protein Tyrosine Kinase 2; RAC, Ras-related C3 botulinum toxin substrate; RAC1, Ras-related C3 botulinum toxin substrate 1; RAF, RAF proto-oncogene serine/threonine-protein kinase; RANTES, C-C Motif Chemokine Ligand 5; Rap1, Ras-related protein Rap-1; Ras, Ras protein family member; RasGrP2, RAS Guanyl Releasing Protein 2; RHOA, Ras homolog gene family, member A; RIP1, Receptor-interacting serine/threonine-protein kinase 1; ROCK, Rho-associated protein kinase; ROS, Reactive oxygen species; Shc, SHC-transforming protein; SOS, Son of sevenless; SRC, Tyrosine-protein kinase; STAT, Signal transducer and activator of transcription protein; TAB, TAK1-binding protein; TAK1, Mitogen-activated protein kinase kinase kinase 7; TBK1, TANK-binding kinase 1; TH0, Naive T cell; TH2, T helper 2 cell, TIAM1, T-lymphoma invasion and metastasis-inducing protein 1; TIRAP, Toll-interleukin 1 receptor (TIR) domain-containing adaptor protein; TLR, Toll-like receptor; TNFα, Tumor necrosis factor alpha; TOLLIP, Toll-interacting protein; Tpl2, Tumor progression locus 2; TRAF3, TNF receptor-associated factor 3; TRAF6, TNF receptor-associated factor 6; TRAM, Translocating chain-associated membrane protein; TRIF, Toll-like receptor adapter molecule 1; WASP, Wiskott-Aldrich syndrome protein; VAV, Guanine nucleotide exchange factor VAV.
      Figure thumbnail fx3
      Fig E3Landscape of expressed variants in AR. Short nonsynonymous transcript variants identified in at least 2 subjects with AR but in none of the control cases. Variants are grouped by gene. Column annotations from top to bottom: subject group, sampling season, and sampling year. Barplot at the right indicate the average expression level of the gene in subjects with AR and control subjects expressed as CPM.
      Figure thumbnail fx4
      Fig E4Microbial variation. A, Relative abundances of microbial (archaeal, bacterial, and viral) genera and total microbial load per sample. Only genera accounting for more than 5% of the total microbial load are shown. The line denotes the total number of microbial reads per sample expressed as CPM. B, Principal coordinates analysis plot of microbial community structure based on Bray-Curtis distance. C, A heatmap of differentially abundant microbial species (Q value ≤ 0.05) between any season- and group-matched year comparison. The different shades of blue illustrate the variance stabilizing transformation values. Darker shades of blue indicate higher values. The species were hierarchically clustered using average linkage and distance metric Pearson correlation. Acronym S1 in the sample name stands for first spring (May), S2 for the second spring, W1 for the first winter (November), and W2 for the second winter. Comparisons that reached the statistical significance are indicated by black boxes next to the heatmap.
      Figure thumbnail fx5
      Fig E5Alpha diversity of epithelial microbiotas across study groups during winter and spring. Alpha diversity indices of different sample groups using various alpha diversity metrics. Shown in the figure are alpha diversities computed using Shannon (A and B), ACE (C and D), Chao1 (E and F), observed (G and H), Fisher (I and J), InvSimpson (K and L), and Simpson (M and N) metrics with (Fig E5, B, D, F, H, J, L, and N) and without (Fig E5, A, C, E, G, I, K, and M) rarefaction of samples to the minimum sampling depth. In the figure, Nov1 stands for the first winter, Nov2 for the second winter, May1 for the first spring, and May2 for the second spring sampling point. ACE, Abundance-based coverage estimators.
      Figure thumbnail fx6
      Fig E6Alpha diversity of epithelial microbiotas between groups during winter. Alpha diversity indices of different groups during 2 consecutive winters using various alpha diversity metrics. Shown in the figure are alpha diversities computed using Shannon (A and B), ACE (C and D), Chao1 (E and F), observed (G and H), Fisher (I and J), InvSimpson (K and L), and Simpson (M and N) metrics with (Fig E6, B, D, F, H, J, L, and N) and without (Fig E6, A, C, E, G, I, K, and M) rarefaction of samples to the minimum sampling depth. In the figure, Nov1 stands for the first winter and Nov2 for the second winter sampling point. ACE, Abundance-based coverage estimators.
      Figure thumbnail fx7
      Fig E7Alpha diversity of epithelial microbiotas between groups during spring. Alpha diversity indices of different groups during 2 consecutive springs sampling using various alpha diversity metrics. Shown in the figure are alpha diversities computed using Shannon (A and B), ACE (C and D), Chao1 (E and F), observed (G and H), Fisher (I and J), InvSimpson (K and L), and Simpson (M and N) metrics with (Fig E7, B, D, F, H, J, L, and N) and without (Fig E7, A, C, E, G, I, K, and M) rarefaction of samples to the minimum sampling depth. In the figure, May1 stands for the first spring and May2 for the second spring sampling point. ACE, Abundance-based coverage estimators.
      Table E1Baseline characteristics of the subjects
      CharacteristicAll subjectsPatients with AR
      ControlsPatients with ARP valueAR-noAITAR-AITP value
      Baseline characteristics
       No. of subjects5633
       Age (y), median (IQR)44 (39-48)39 (34-46).4334 (32-42)41 (39-43).70
       Men/women (n)1/43/3.542/11/21.00
      Spirometry values (% predicted), median (IQR)
       FEV1 baseline107 (105-109)95 (94-96).2396 (95-103)97 (94-95).30
       FEV1 bronchodilation110 (105-110)97 (95-99).3997 (96-104)97 (95-99)1.00
       FEV1/FVC ratio baseline104 (101-104)98 (98-103).33103 (101-107)95 (92-98).30
       FEV1/FVC bronchodilation105 (101-106)102 (98-104).55104 (101-108)99.5 (97-102).40
       PEF baseline109 (102-125)109 (101-110).57110 (106-116)98 (87-109).20
       PEF bronchodilation104 (92-125)112 (92-114).74112 (102-113)100.5 (86-115)1.00
      Symptoms during sampling in spring1, median (IQR)
       Total RQLQ score0 (0-6)48.5 (24-74).01940 (23-49)74 (49-87).40
       Total VAS score15 (10-36)744 (358-1352).009217 (358-471)1352 (1128-1512).10
      Serum values in spring1, median (IQR)
       Total IgE (kU/L)9 (7-14)61.5 (39-200).01539 (26-275)72 (62-136).70
       Birch specific IgE (kU/L)0 (0-0)16.4 (3.3-24.2).00412 (6-16)24 (14-28).40
       Timothy specific IgE (kU/L)0 (0-0)0 (0-0).850 (0-45)0 (0-0)1.00
      SPT wheal diameter (mm), median (IQR)
       Negative control0 (0-0)0 (0-0)1.000 (0-0)0 (0-0)1.00
       Histamine (positive control)5 (5-5)5 (5-5)1.005 (5-5)5 (5-5)1.00
       Birch pollen0 (0-0)5 (4-6).0044 (3.5-6)5 (5-5.5).60
       Timothy grass pollen0 (0-0)0 (0-0)1.000 (0-2.5)0 (0-0)1.00
      Festuca pratensis pollen0 (0-0)0 (0-0)1.000 (0-2.5)0 (0-0)1.00
       Mugwort pollen0 (0-0)0 (0-0)1.000 (0-2.5)0 (0-0)1.00
      Cladosporium herbarum0 (0-0)0 (0-0)1.000 (0-0)0 (0-0)1.00
       Cat dander0 (0-0)0 (0-5).460 (0-2.5)0 (0-3)1.00
       Dog dander0 (0-0)0 (0-4).460 (0-2)0 (0-2.5)1.00
       Horse dander0 (0-0)0 (0-0)1.000 (0-0)0 (0-0)1.00
      Dermatophagoides pteronyssinus0 (0-0)0 (0-0)1.000 (0-0)0 (0-0)1.00
      FVC, Forced vital capacity; IQR, interquartile range; PEF, peak expiratory force; RQLQ, Rhinoconjunctivitis Quality of Life Questionnaire; SPT, skin prick test; VAS, visual analogue scale.
      Diagnosis of AR was based on a typical history, SPT (ALK-Abelló, Hørsholm, Denmark), total serum IgE, and serum birch and timothy allergen specific IgE antibodies. Healthy volunteers did not have symptoms and were negative for SPT of common aeroallergens and serum birch and timothy allergen specific IgE antibodies. Exclusion criteria were age less than 12 y, use of tobacco products, non-AR, AR symptoms caused by other than seasonal allergens, asthma, and general disease requiring regular medication. Asthma was excluded by absence of typical symptoms and by normal values in spirometry with bronchodilation test.
      • Mattila P.
      • Renkonen J.
      • Toppila-Salmi S.
      • Parviainen V.
      • Joenvaara S.
      • Alff-Tuomala S.
      • et al.
      Time-series nasal epithelial transcriptomics during natural pollen exposure in healthy subjects and allergic patients.
      One subject starting pollen AIT was tested for bronchial hyperresponsiveness by histamine challenge test. He had normal 15th percentile density of FEV1 result instead of bronchodilation test. Questionnaire regarding baseline characteristics and symptoms was collected at each sampling visit. We used the 28-item RQLQ.
      • Juniper E.F.
      • Thompson A.K.
      • Ferrie P.J.
      • Roberts J.N.
      Development and validation of the mini Rhinoconjunctivitis Quality of Life Questionnaire.
      Subjects filled the VAS so that they were without medication. The questionnaire included 18 questions concerning airway symptoms and 23 questions concerning general health. Value 0 (mm) indicated no symptoms, and value 100 (mm) indicated the worst case. The total maximum score of the 41 questions was 4100. In the analysis, VAS scores less than or equal to 3 mm were regarded as 0. P values were computed by Kruskal-Wallis and Mann Whitney U tests (continuous variables) or by Fisher exact test (dichotomous variables).

      References

        • Poole A.
        • Urbanek C.
        • Eng C.
        • Schageman J.
        • Jacobson S.
        • O'Connor B.P.
        • et al.
        Dissecting childhood asthma with nasal transcriptomics distinguishes subphenotypes of disease.
        J Allergy Clin Immunol. 2014; 133: 670-678.e12
        • Joenvaara S.
        • Mattila P.
        • Renkonen J.
        • Makitie A.
        • Toppila-Salmi S.
        • Lehtonen M.
        • et al.
        Caveolar transport through nasal epithelium of birch pollen allergen Bet v 1 in allergic patients.
        J Allergy Clin Immunol. 2009; 124 (e1-21): 135-142
        • Toppila-Salmi S.
        • van Drunen C.M.
        • Fokkens W.J.
        • Golebski K.
        • Mattila P.
        • Joenvaara S.
        • et al.
        Molecular mechanisms of nasal epithelium in rhinitis and rhinosinusitis.
        Curr Allergy Asthma Rep. 2015; 15: 495
        • Choi C.H.
        • Poroyko V.
        • Watanabe S.
        • Jiang D.
        • Lane J.
        • deTineo M.
        • et al.
        Seasonal allergic rhinitis affects sinonasal microbiota.
        Am J Rhinol Allergy. 2014; 28: 281-286
        • Chiu C.Y.
        • Chan Y.L.
        • Tsai Y.S.
        • Chen S.A.
        • Wang C.J.
        • Chen K.F.
        • et al.
        Airway microbial diversity is inversely associated with mite-sensitized rhinitis and asthma in early childhood.
        Sci Rep. 2017; 7: 1820
        • Bousquet J.
        • Khaltaev N.
        • Cruz A.A.
        • Denburg J.
        • Fokkens W.J.
        • Togias A.
        • et al.
        Allergic Rhinitis and its Impact on Asthma (ARIA) 2008 update (in collaboration with the World Health Organization, GA(2)LEN and AllerGen).
        Allergy. 2008; 63: 8-160
        • Akdis C.A.
        • Akdis M.
        Mechanisms of allergen-specific immunotherapy.
        J Allergy Clin Immunol. 2011; 127 (quiz 28-9): 18-27
        • Lal D.
        • Keim P.
        • Delisle J.
        • Barker B.
        • Rank M.A.
        • Chia N.
        • et al.
        Mapping and comparing bacterial microbiota in the sinonasal cavity of healthy, allergic rhinitis, and chronic rhinosinusitis subjects.
        Int Forum Allergy Rhinol. 2017; 7: 561-569