Bronchial epithelial cells (BECs) from children with asthma were obtained under an IRB approved study and were differentiated to an organotypic epithelium and primed via co-culture with healthy donor Th17 lymphocytes for 4 days prior to apical infection with human rhinovirus-16 (RV-16) and RNA harvest 48 hours later. RNA sequencing with WGCNA analysis was performed to identify modules of gene expression altered by Th17 priming or RV-16 infection in BECs or Th17 cells. Gene expression was correlated with viral copy number and with secreted protein levels.

BEC cultures and differentiation: BECs were obtained from children with asthma undergoing elective surgery under studies #12490 and #1596 approved by the Seattle Children’s Hospital Institutional Review Board as we have previously described.¹ Parents of subjects provided written consent and children over 7 years of age provided assent. Basal cells were selected and proliferated under submerged culture conditions using PneumaCult™EX-Plus medium (Stemcell™). Passage 2-3 BECs were differentiated for ≥ 21 days at an air-liquid interface (ALI) to generate organotypic, pseudostratified epithelial cultures.¹

BEC-Immune Cell Co-culture: Primary human Th17 cells (STEMCELL™) from a healthy donor were suspended (25,000 cells/mL) in a 50:50 mixture of PneumaCult™ ALI medium and RPMI 1640 with 10% FBS and added to the basolateral chamber of transwell cultures below differentiated BECs with media changed every 48 hours.

RV-16 infection and sample collection: BECs were infected with human rhinovirus A16 (RV-A16, Source: ATCC®)) at a multiplicity of infection (MOI) of 0.5 on the apical surface for 2 hours at 34°C as we have previously described.² For BEC-Immune cell co-culture, infection on the apical surface occurred on day 4 after the start of Th17-BEC co-culture. (Figure 1) At 24 hours following RV-16 infection (or 5 days after co-culture initiation for uninfected cultures), Th17 cells were removed for RNA isolation and human primary neutrophils from a healthy donor (iQ Biosciences®; 25,000 cells/mL) in a 50:50 mixture of PneumaCult™-ALI medium and RPMI 1640 with 10% FBS were added to the basolateral chamber below differentiated BECs. 24 hours later neutrophils and BECs were harvested independently for RNA isolation, and media was collected.

Flow Cytometry: Primary human Th17 cells (STEMCELL^TM) from a healthy donor were stained with fixable viability dye (Zombie Violet, Biolegend cat. 423114) before staining with antibodies against CD4 (APC-Cy7, clone OKT4, Biolegend cat 317417), CCR6 (BV711, clone G034E3, Biolegend cat 353435), and CXCR3 (PE/Dazzle^TM 594, clone G025H7, Biolegend cat 353735). For intracellular staining, cells were stimulated with Cell Stimulation Cocktail (ThermoFisher, cat 00-4970-93) and Protein Transport Inhibitor Cocktail (ThermoFisher, cat 00-4780-93) for 4 hours and viability stained before fixation and permeabilization with Foxp3/Transcription Factor Staining Buffer set according to the manufacturer’s instructions (ThermoFisher, cat 00-5523-00). Permeabilized cells were stained with antibodies against RORgt (PE, clone AFKJS-9, ThermoFisher cat 12-6988-82) and IL-17A (APC, clone BL168, Biolegend cat 512333). Data was collected using a BD FACSymphony^TM A5 Cell Analyzer. Data was analyzed using FlowJo  v10 (BD). Characterization confirmed that 86.2% of cells were CD4+/CCR6+/CXCR3- and 99.3% positive for RORgT and 12.9% positive for IL-17A after stimulation (FACS_Figure.png)

Protein measurement: In collected media, protein concentrations of IL-1β and IFN-α were measured via Human Luminex® Assay (R&D®).

RNA sequencing: RNA was sequenced on a NextSeq 2000 sequencer (Illumina) with paired-end 53-base reads at a target depth of 5 million reads per sample, as previously described.¹ Gene counts were generated as previously described using voomWithQualityWeights from the limma R package after alignment and quality-control.¹ Samples that had human aligned counts greater than 1 million mapped reads and a median coefficient of variation coverage less than 0.7 were kept for downstream analyses. Data from each cell type were processed separately. All raw bulk RNA-sequencing data available at NCBI GEO GSE314289

Rhinovirus viral load quantification: To quantify rhinovirus viral load, the raw sequencing reads were aligned to the human reference genome (GRCh38) using STAR (v2.7.0a).³ Reads that were unmapped to the human reference genome (GRCh38) were used for viral quantification using the Salmon (v1.10.2) quasi-mapping algorithm against a custom viral reference database of 42 complete viral genomes (Supplemental Table 1).⁴ Transcript-level abundance estimates of number of reads were used to quantify rhinovirus abundance in each sample.

Statistics: Differentially expressed genes were identified using mixed effect linear models comparing expression differences by cell culture condition and infection status using the R package “kimma”.⁵ The model syntax was: gene_expression ~ condition*infection and including a random effect for epithelial cell donorID. This identified 4,030 genes that reached a false discovery rate adjusted p-value (FDR) < 0.05 by the Benjamini-Hochberg correction procedure for cell culture condition. These genes were then utilized for supervised weighted gene co-expression network analysis (WGCNA) to identify genes with similar expression patterns and group them into modules. Module values were summarized by taking the mean of all genes (log2 transformed values) in the respective module. These modules were then modeled using the same linear model syntax to identify those showing differential expression patterns by cell culture condition and infection status at a stringent FDR < 0.05 for each pairwise comparison. Multiple hypothesis testing correction was performed using the Benjamini-Hochberg procedure. Pathway enrichment was performed using the clusterprofiler R package, which calculates a hypergeometric FDR corrected p-value for enrichment of public genesets.⁶ We used the gene ontology biological processes (GO_BP), KEGG, Reactome, Biocarta, and MSigDB Hallmark genesets for enrichment.^7–12 The modules presented in the main text were annotated based on manual inspection of the enrichment terms and module genes. A similar model as used for gene expression was used for protein levels with multiple hypothesis testing correction using the Benjamini-Hochberg correction procedure. Pearson correlation of module expression in BECs to module expression in Th17 cells and for protein level to module expression was performed.

1.         Doni Jayavelu N, Altman MC, Benson B, Dufort MJ, Vanderwall ER, Rich LM, White MP, Becker PM, Togias A, Jackson DJ, Debley JS. Type 2 inflammation reduces SARS-CoV-2 replication in the airway epithelium in allergic asthma through functional alteration of ciliated epithelial cells. J Allergy Clin Immunol. 2023 Jul;152(1):56–67. PMCID: PMC10052850

2.         Vanderwall ER, Barrow KA, Rich LM, Read DF, Trapnell C, Okoloko O, Ziegler SF, Hallstrand TS, White MP, Debley JS. Airway epithelial interferon response to SARS-CoV-2 is inferior to rhinovirus and heterologous rhinovirus infection suppresses SARS-CoV-2 replication. Sci Rep. 2022 Apr 28;12(1):6972. PMCID: PMC9048621

3.         Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013 Jan 1;29(1):15–21. PMCID: PMC3530905

4.         Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017 Apr;14(4):417–419. PMCID: PMC5600148

5.         Dill-McFarland KA, Mitchell K, Batchu S, Segnitz RM, Benson B, Janczyk T, Cox MS, Mayanja-Kizza H, Boom WH, Benchek P, Stein CM, Hawn TR, Altman MC. Kimma: flexible linear mixed effects modeling with kinship covariance for RNA-seq data. Bioinformatics. 2023 May 4;39(5):btad279. PMCID: PMC10182851

6.         Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, Fu X, Liu S, Bo X, Yu G. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). 2021 Aug 28;2(3):100141. PMCID: PMC8454663

7.         Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015 Dec 23;1(6):417–425. PMCID: PMC4707969

8.         Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011 Jun 15;27(12):1739–1740. PMCID: PMC3106198

9.         Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000 Jan 1;28(1):27–30. PMCID: PMC102409

10.       Gene Ontology Consortium, Aleksander SA, Balhoff J, Carbon S, Cherry JM, Drabkin HJ, Ebert D, Feuermann M, Gaudet P, Harris NL, Hill DP, Lee R, Mi H, Moxon S, Mungall CJ, Muruganugan A, Mushayahama T, Sternberg PW, Thomas PD, Van Auken K, Ramsey J, Siegele DA, Chisholm RL, Fey P, Aspromonte MC, Nugnes MV, Quaglia F, Tosatto S, Giglio M, Nadendla S, Antonazzo G, Attrill H, Dos Santos G, Marygold S, Strelets V, Tabone CJ, Thurmond J, Zhou P, Ahmed SH, Asanitthong P, Luna Buitrago D, Erdol MN, Gage MC, Ali Kadhum M, Li KYC, Long M, Michalak A, Pesala A, Pritazahra A, Saverimuttu SCC, Su R, Thurlow KE, Lovering RC, Logie C, Oliferenko S, Blake J, Christie K, Corbani L, Dolan ME, Drabkin HJ, Hill DP, Ni L, Sitnikov D, Smith C, Cuzick A, Seager J, Cooper L, Elser J, Jaiswal P, Gupta P, Jaiswal P, Naithani S, Lera-Ramirez M, Rutherford K, Wood V, De Pons JL, Dwinell MR, Hayman GT, Kaldunski ML, Kwitek AE, Laulederkind SJF, Tutaj MA, Vedi M, Wang SJ, D’Eustachio P, Aimo L, Axelsen K, Bridge A, Hyka-Nouspikel N, Morgat A, Aleksander SA, Cherry JM, Engel SR, Karra K, Miyasato SR, Nash RS, Skrzypek MS, Weng S, Wong ED, Bakker E, Berardini TZ, Reiser L, Auchincloss A, Axelsen K, Argoud-Puy G, Blatter MC, Boutet E, Breuza L, Bridge A, Casals-Casas C, Coudert E, Estreicher A, Livia Famiglietti M, Feuermann M, Gos A, Gruaz-Gumowski N, Hulo C, Hyka-Nouspikel N, Jungo F, Le Mercier P, Lieberherr D, Masson P, Morgat A, Pedruzzi I, Pourcel L, Poux S, Rivoire C, Sundaram S, Bateman A, Bowler-Barnett E, Bye-A-Jee H, Denny P, Ignatchenko A, Ishtiaq R, Lock A, Lussi Y, Magrane M, Martin MJ, Orchard S, Raposo P, Speretta E, Tyagi N, Warner K, Zaru R, Diehl AD, Lee R, Chan J, Diamantakis S, Raciti D, Zarowiecki M, Fisher M, James-Zorn C, Ponferrada V, Zorn A, Ramachandran S, Ruzicka L, Westerfield M. The Gene Ontology knowledgebase in 2023. Genetics. 2023 May 4;224(1):iyad031. PMCID: PMC10158837

11.       Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000 May;25(1):25–29. PMCID: PMC3037419

12.       Milacic M, Beavers D, Conley P, Gong C, Gillespie M, Griss J, Haw R, Jassal B, Matthews L, May B, Petryszak R, Ragueneau E, Rothfels K, Sevilla C, Shamovsky V, Stephan R, Tiwari K, Varusai T, Weiser J, Wright A, Wu G, Stein L, Hermjakob H, D’Eustachio P. The Reactome Pathway Knowledgebase 2024. Nucleic Acids Res. 2024 Jan 5;52(D1):D672–D678. PMCID: PMC10767911