Birch pollen allergen immunotherapy reprograms nasal epithelial transcriptome and recovers microbial diversity
Authors: Tanzeela Hanif, Kishor Dhaygude, Matti Kankainen, Jutta Renkonen, Pirkko Mattila, Teija Ojala, Sakari Joenväärä, Mika Mäkelä, Anna Pelkonen, Paula Kauppi, Tari Haahtela, Risto Renkonen, Sanna Toppila-Salmi
Journal: Journal of Allergy and Clinical Immunology (2019)

Airway epithelial cells are known to have an important role in allergic rhinitis (AR).1Poole 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.e12Abstract Full Text Full Text PDF PubMed Scopus (161) Google Scholar, 2Joenvaara 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-142Abstract Full Text Full Text PDF PubMed Scopus (38) Google Scholar, 3Toppila-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: 495Crossref PubMed Scopus (33) Google Scholar They constitute the first line of defense against inhaled aeroallergens and are active mediators of innate and adaptive immune responses.3Toppila-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: 495Crossref PubMed Scopus (33) Google Scholar Their aberrant functioning is linked with an intake of allergens,2Joenvaara 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-142Abstract Full Text Full Text PDF PubMed Scopus (38) Google Scholar and their transcriptome is reprogrammed under exposure to pollens2Joenvaara 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-142Abstract Full Text Full Text PDF PubMed Scopus (38) Google Scholar, 3Toppila-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: 495Crossref PubMed Scopus (33) Google Scholar as well as in AR3Toppila-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: 495Crossref PubMed Scopus (33) Google Scholar and atopic asthma.1Poole 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.e12Abstract Full Text Full Text PDF PubMed Scopus (161) Google Scholar Furthermore, epithelial cells interact with and are involved in generating an environmental niche for the respiratory microbiota, whose imbalance has been associated with seasonal AR4Choi 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-286Crossref PubMed Scopus (42) Google Scholar and childhood rhinitis and asthma.5Chiu 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: 1820Crossref PubMed Scopus (47) Google Scholar 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,6Bousquet 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-160Crossref PubMed Scopus (3693) Google Scholar decrease in allergen-specific biomarkers, and altered T- and B-cell responses.7Akdis C.A. Akdis M. Mechanisms of allergen-specific immunotherapy.J Allergy Clin Immunol. 2011; 127 (quiz 28-9): 18-27Abstract Full Text Full Text PDF PubMed Scopus (340) Google Scholar 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). 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). 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.3Toppila-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: 495Crossref PubMed Scopus (33) Google Scholar, 7Akdis C.A. Akdis M. Mechanisms of allergen-specific immunotherapy.J Allergy Clin Immunol. 2011; 127 (quiz 28-9): 18-27Abstract Full Text Full Text PDF PubMed Scopus (340) Google Scholar 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,8Lal 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-569Crossref PubMed Scopus (50) Google Scholar 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 AR4Choi 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-286Crossref PubMed Scopus (42) Google Scholar but similar to that focusing on children with asthma and rhinitis.5Chiu 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: 1820Crossref PubMed Scopus (47) Google Scholar 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 season4Choi 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-286Crossref PubMed Scopus (42) Google Scholar 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. 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 epithelial brushing was performed to middle meatus of both sides of nasal cavity after slight blowing of nose without local anesthesia as described.E1Mattila 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.Allergy. 2010; 65: 175-183Crossref PubMed Scopus (29) Google Scholar 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. 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 were preprocessed as described previously.E2Kumar 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.BMC Genomics. 2017; 18: 629Crossref PubMed Scopus (32) Google Scholar Briefly, TrimmomaticsE3Bolger A.M. Lohse M. Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data.Bioinformatics. 2014; 30: 2114-2120Crossref PubMed Scopus (28414) Google Scholar 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 STARE4Dobin A. Davis C.A. Schlesinger F. Drenkow J. Zaleski C. Jha S. et al.STAR: ultrafast universal RNA-seq aligner.Bioinformatics. 2013; 29: 15-21Crossref PubMed Scopus (19483) Google Scholar 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,E5Liao Y. Smyth G.K. Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote.Nucleic Acids Res. 2013; 41: e108Crossref PubMed Scopus (1507) Google Scholar feature counts were converted to expression estimates using Trimmed Mean of M-values normalization,E6Robinson M.D. Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data.Genome Biol. 2010; 11: R25Crossref PubMed Scopus (4122) Google Scholar 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. Differential expression testing was performed using the edgeRE7Robinson M.D. McCarthy D.J. Smyth G.K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics. 2010; 26: 139-140Crossref PubMed Scopus (21096) Google Scholar 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.E8Kolde R. pheatmap: Pretty Heatmaps. R package version 1.0.10.2018Google Scholar 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.E9Chen H. VennDiagram: generate high-resolution Venn and Euler plots. R package version 1.6.20.2018Google Scholar Functional profiles of DEGs were investigated with clusterProfilerE10Yu G. Wang L.G. Han Y. He Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS. 2012; 16: 284-287Crossref PubMed Scopus (12094) Google Scholar 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.E11Luo W. Brouwer C. Pathview: an R/Bioconductor package for pathway-based data integration and visualization.Bioinformatics. 2013; 29: 1830-1831Crossref PubMed Scopus (946) Google Scholar 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. Transcript variants were called from STAR alignments using the GATK best practices workflows for transcriptome dataE12McKenna 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.Genome Res. 2010; 20: 1297-1303Crossref PubMed Scopus (14912) Google Scholar and then annotated using AnnovarE13Wang K. Li M. Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res. 2010; 38: e164Crossref PubMed Scopus (7948) Google Scholar as defined previously.E2Kumar 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.BMC Genomics. 2017; 18: 629Crossref PubMed Scopus (32) Google Scholar Quality control analyses were performed as defined previously.E2Kumar 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.BMC Genomics. 2017; 18: 629Crossref PubMed Scopus (32) Google Scholar 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 AnnovarE13Wang K. Li M. Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res. 2010; 38: e164Crossref PubMed Scopus (7948) Google Scholar annotation tool. Heatmap was plotted using pheatmap. The functional effects of variants were taken from AnnovarE13Wang K. Li M. Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res. 2010; 38: e164Crossref PubMed Scopus (7948) Google Scholar outputs. In addition, we plotted barplot using CPM expression value of genes in healthy control and AR groups. Microbial community profiling was performed as previously describedE14Dufva 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.Nat Commun. 2018; 9: 1567Crossref PubMed Scopus (78) Google Scholar 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.E3Bolger A.M. Lohse M. Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data.Bioinformatics. 2014; 30: 2114-2120Crossref PubMed Scopus (28414) Google Scholar Paired-end reads passing the preprocessing were mapped against rRNA sequences from RFAME15Nawrocki 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.Nucleic Acids Res. 2015; 43: D130-D137Crossref PubMed Scopus (725) Google Scholar v12.3 using the Burrows-Wheeler AlignerE16Li H. Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform.Bioinformatics. 2009; 25: 1754-1760Crossref PubMed Scopus (26904) Google Scholar with default settings, and reads-matching rRNAs were filtered by using samtools.E17Li H. Handsaker B. Wysoker A. Fennell T. Ruan J. Homer N. et al.The Sequence Alignment/Map format and SAMtools.Bioinformatics. 2009; 25: 2078-2079Crossref PubMed Scopus (31884) Google Scholar CentrifugeE18Kim D. Song L. Breitwieser F.P. Salzberg S.L. Centrifuge: rapid and sensitive classification of metagenomic sequences.Genome Res. 2016; 26: 1721-1729Crossref PubMed Scopus (531) Google Scholar 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 RefSeqE19O'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.Nucleic Acids Res. 2016; 44: D733-D745Crossref PubMed Scopus (2648) Google Scholardatabase 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 DeSeq2E20Love M.I. Huber W. Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.Genome Biol. 2014; 15: 550Crossref PubMed Scopus (33027) Google Scholar 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 adjustmentE21Storey J.D. The positive false discovery rate: a Bayesian interpretation and the q-value.Ann Stat. 2003; 31: 2013-2035Crossref Scopus (1470) Google Scholar 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 softwareE22McMurdie P.J. Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data.PLoS One. 2013; 8: e61217Crossref PubMed Scopus (8002) Google Scholar applied on the number of reads assigned directly to the given taxonomic level. 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). 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.7Akdis C.A. Akdis M. Mechanisms of allergen-specific immunotherapy.J Allergy Clin Immunol. 2011; 127 (quiz 28-9): 18-27Abstract Full Text Full Text PDF PubMed Scopus (340) Google Scholar 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. 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). 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). 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 interleukinE23Pawankar R. Mori S. Ozu C. Kimura S. Overview on the pathomechanisms of allergic rhinitis.Asia Pac Allergy. 2011; 1: 157-167Crossref PubMed Google Scholar IL-4, which is IgE synthesis switch factor,E23Pawankar R. Mori S. Ozu C. Kimura S. Overview on the pathomechanisms of allergic rhinitis.Asia Pac Allergy. 2011; 1: 157-167Crossref PubMed Google Scholar and IL-5, which is an eosinophil growth factor.E23Pawankar R. Mori S. Ozu C. Kimura S. Overview on the pathomechanisms of allergic