Airway Smooth Muscle-Specific Transcriptomic Signatures of Glucocorticoid Exposure.

Glucocorticoids, commonly used asthma controller medications, decrease symptoms in most patients, but some remain symptomatic despite high-dose treatment. The physiological basis underlying the glucocorticoid response, especially in asthma patients with severe, refractory disease, is not fully understood. We sought to identify differences between the transcriptomic response of airway smooth muscle (ASM) cells derived from donors with fatal asthma and donors without asthma to glucocorticoid exposure and to compare ASM-specific changes with those observed in other cell types. In cells derived from nine donors with fatal asthma and eight donors without asthma, RNA sequencing was used to measure ASM transcriptome changes after exposure to budesonide (100 nM 24 h) or control vehicle (DMSO). Differential expression results were obtained for this dataset, as well as 13 publicly available glucocorticoid-response transcriptomic datasets corresponding to seven cell types. Specific genes were differentially expressed in response to glucocorticoid exposure (7,835 and 6,957 in ASM cells derived from donors with fatal asthma and donors without asthma, respectively; adjusted P value < 0.05). Transcriptomic changes in response to glucocorticoid exposure were similar in ASM derived from donors with fatal asthma and donors without asthma, with enriched ontological pathways that included cytokine- and chemokine-related categories. A comparison of glucocorticoid-induced changes in the nonasthma ASM transcriptome with those observed in six other cell types showed that ASM has a distinct glucocorticoid-response signature that is also present in ASM cells from donors with fatal asthma.

Clinical Relevance

Glucocorticoid treatment results in strong gene expression changes in airway smooth muscle, some of which are cell-type specific. Transcriptome changes in response to budesonide treatment are similar in airway smooth muscle derived from donors with fatal asthma versus donors without asthma.

Asthma is an episodic, inflammatory respiratory disease that is characterized by increased airway responsiveness to specific environmental stimuli and affects more than 25 million Americans (1). Glucocorticoids are medications that are commonly used to treat asthma. Inhaled glucocorticoids that act directly in the lung are prescribed to most individuals with persistent asthma to decrease symptoms and exacerbations; however, some patients require long-term use of oral glucocorticoids to control their disease (2). Individuals with severe, refractory disease remain symptomatic despite high-dose glucocorticoid treatment (3), and although these patients represent a heterogeneous group, features that are shared by most include glucocorticoid insensitivity and, in part, irreversible airflow obstruction (4). A better understanding of the physiological basis underlying glucocorticoid responsiveness would increase our ability to identify and improve therapeutic options for individuals with severe asthma, as well as to obtain insights that apply to patients with milder disease.

The airway smooth muscle (ASM) is a target of glucocorticoid medications (5) that regulates airway narrowing (6) and contributes to the airway remodeling observed in severe asthma by having increased mass and cell size (7). In addition to directly reducing inflammation (8), glucocorticoid treatment affects other asthma-related phenotypes involving ASM, including bronchodilation (9), airway hyperresponsiveness (10), and contractility (11). Studies of primary human ASM cells in vitro have shown retention of smooth-muscle–specific protein expression and agonist-induced calcium mobilization, force generation, and relaxation responses that are important in asthma (12, 13). Furthermore, cells derived from individuals with fatal asthma retain a unique phenotype with differences in proliferative, contractile, and transcriptomic outcomes that are sustained over at least five passages (13–17). Thus, in vitro studies of the glucocorticoid response in ASM from donors with fatal asthma and donors without asthma offer a unique opportunity to understand asthma-specific processes.

Glucocorticoids act intracellularly by binding to glucocorticoid receptors that then translocate to cell nuclei and modulate the transcription of various genes in a tissue-dependent fashion (18–23). The antiinflammatory action of glucocorticoids occurs in part by 1) homodimer glucocorticoid–glucocorticoid receptor complexes stimulating the transcription of antiinflammatory genes by directly binding to DNA at glucocorticoid response elements, and 2) glucocorticoid–glucocorticoid receptor complexes inhibiting proinflammatory transcription factors such as NFκB (19–22). Other investigators and we have described global glucocorticoid-induced ASM gene expression changes in cells derived from donors without asthma or chronic disease (15, 24–27), but differences in the transcriptome response to glucocorticoids in ASM from individuals with asthma remain inadequately understood. RNA sequencing (RNA-Seq) is a technique that permits comprehensive and in-depth quantification of transcriptomes in a cell or tissue (28) and has been widely used in respiratory research (15, 29, 30). Here, we used RNA-Seq to identify differences in the transcriptome response to glucocorticoid exposure in ASM cells derived from donors with fatal asthma and donors without asthma, and we compared ASM-specific changes with those observed in six other cell types using 13 publicly available transcriptomic datasets.

Methods

Details regarding the methods used in this work are provided in the data supplement.

Ethics Statement

Lung tissue was obtained from the National Disease Resource Interchange and the International Institute for the Advancement of Medicine, and its use was approved by the University of Pennsylvania Institutional Review Board and the Rutgers Institutional Review Board. Use of the cells does not constitute human subjects research because all donor tissues were harvested anonymously and deidentified.

ASM Cell Culture and Treatment

Primary ASM cells were isolated from 17 white, nonsmoking donors (nine who died of fatal asthma and eight with no chronic illness or medication use). ASM cells were cultivated as described previously (13, 31). Passages 3 and 4 cells from each donor were exposed to 100 nM budesonide or control vehicle (DMSO) for 24 hours.

RNA-Seq Library Construction and Sequencing

Total RNA was extracted from cells using the miRNAeasy mini kit (Qiagen Sciences, Inc.). The Illumina TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, Inc.) was used to prepare stranded RNA-Seq libraries. Ambion External RNA Controls Consortium RNA Spike-In Control Mix 1 (Life Technologies Corp.) was added to the samples. Sequencing was performed on an Illumina HiSeq 2500 at the University of Pennsylvania Next-Generation Sequencing Core.

RNA-Seq Data Analysis

Taffeta scripts were used to analyze RNA-Seq data and generate HTML summary reports (https://github.com/blancahimes/taffeta) (15, 16). Gene-level differential expression analyses comparing 1) budesonide versus control in ASM derived from donors without asthma, 2) budesonide versus control in ASM derived from donors with fatal asthma, 3) ASM from donors with fatal asthma versus ASM from donors without asthma exposed to control, and 4) ASM from donors with fatal asthma versus ASM from donors without asthma exposed to budesonide were performed with DESeq2 (32). The RNA-Seq data are available from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE94335.

Gene Set Enrichment Analysis

Gene set enrichment analysis was performed with the fast gene set enrichment analysis algorithm (33). Gene sets of KEGG and Reactome pathway annotations were downloaded from MsigDB collections (http://software.broadinstitute.org/gsea/msigdb/collections.jsp). P values were assessed based on the distribution of enrichment scores after 10,000 permutations, and q values were obtained using the Benjamini-Hochberg approach.

Analysis of Publicly Available Glucocorticoid-Response Transcriptomic Datasets

We sought publicly available transcriptomic datasets that measured the effect of glucocorticoid exposure on various human cell types to determine ASM-specific transcriptomic changes induced by glucocorticoids. Gene Expression Omnibus searches of terms related to glucocorticoids yielded 13 glucocorticoid-response datasets. Differential expression results from a glucocorticoid versus control comparison for individual studies were obtained using RAVED (https://github.com/HimesGroup/raved). Full differential expression results for each individual study are available from our previously developed app, REALGAR (http://realgar.org/) (34).

Tissue-Specific Integration of Glucocorticoid-Response Transcriptomic Data

For cell types that had more than one transcriptomic dataset, we obtained integrated results using effect-size– and rank-based approaches. Gene significance per cell type was determined based on having q value < 0.05 in both integration procedures.

Hierarchical Clustering

Hierarchical clustering was performed to compare cell-type–specific gene expression patterns across glucocorticoid versus control comparisons, using log2 fold changes in expression of 56 genes that were significant in nonasthma ASM and also had results available in all datasets considered.

Results

Budesonide Exposure Modulates the Transcriptome of ASM Cells Derived from Donors with Fatal Asthma and Donors without Asthma

The characteristics of the 17 ASM donors are provided in Table 1. Cells from each donor were exposed to budesonide or control vehicle, yielding 34 RNA-Seq samples that were deemed to be of sufficiently high quality to include in differential expression analyses after quality control measures were obtained (Figure E1 and Table E1 in the data supplement). Analyses of transcriptomic changes in ASM that had been exposed to budesonide versus control revealed 7,835 differentially expressed genes in ASM derived from donors without asthma (Figure 1A) and 6,957 in ASM derived from donors with fatal asthma (Figure 1B), according to a 5% false discovery rate threshold (i.e., q < 0.05). Most of the genes that were differentially expressed in response to budesonide exposure (i.e., 5,664) were statistically significant in ASM derived from both donors with fatal asthma and donors without asthma (Figure 1C), with a consistent direction of effect (Figure 1D). The top-ranked significantly differentially expressed genes (i.e., having a q value < 10−50 and an absolute log2 fold change > 3) for ASM derived from donors with fatal asthma and donors without asthma are provided in Tables 2 and 3, respectively. Among these glucocorticoid-responsive genes were well-known transcripts that were previously identified in various tissues (e.g., FKBP5 [35], GLUL [36], PER1 [37, 38], and TSC22D3 [35, 39]), as well as more recently investigated ones (CRISPLD2 [15] and KLF15 [24, 27]) in studies of ASM. A summary of published evidence for glucocorticoid responsiveness–related roles of the top-ranked genes listed in Tables 2 and 3 is provided in Table E2. A comparison of ASM from donors with fatal asthma versus donors without asthma yielded no statistically significant differences in gene expression among cells with control exposure (Figure E2A), and two differentially expressed genes (i.e., CCK and PMEL) among cells that had been exposed to budesonide (Figure E2B and Table E3).

Table 1.

Characteristics of the ASM Donors

	Fatal Asthma (n = 9)	Control (n = 8)
Sex
Male	4	5
Female	5	3
Age
Mean ± SD	25.1 ± 14.8	28.6 ± 13.6
[Range]	[9–48]	[17–52]
Body mass index
Mean ± SD	23.0 ± 6.8	27.1 ± 8.8
[Range]	[11.3–34.5]	[21.5–46.9]
Cause of death
Anoxia/overdose	—	2
Anoxic brain injury	—	1
Cerebrovascular accident	—	1
Motor vehicle accident	—	4
Medication
Albuterol	7	0

Definition of abbreviation: ASM = airway smooth muscle.

All donors were white nonsmokers. There were no significant differences in sex, age, or body mass index between donors with fatal asthma and donors without asthma. A medical examiner ruled that the cause of death for donors with fatal asthma was “asthma attack/anoxia,” or a significant asthma event was listed as preceding death. Medications listed as taken by donors near the time of death were provided by family members to medical examiners. Besides albuterol as listed in the table, only two other asthma-related entries were provided: one “fluticasone/salmeterol” and one “asthma meds (unspecified).”

Figure 1.

Effects of budesonide (BUD) exposure on the transcriptome of airway smooth muscle (ASM) derived from donors with fatal asthma and donors without asthma. (A and B) Volcano plots of differential expression results for BUD versus control in ASM derived from donors without asthma (A) and donors with fatal asthma (B). The y-axis corresponds to the negative log (base 10) of Benjamini-Hochberg–corrected P values (i.e., q values) and the x-axis corresponds to the log (base 2) of the fold change for differences in expression. Differentially expressed transcripts with q value < 0.05 are colored in red. (C) Comparison of genes with statistically significant expression changes in BUD versus control in ASM derived from donors with fatal asthma and donors without asthma. The number of genes without statistically significant results is shown on the bottom right. (D) Scatter plot of negative log (base 10) of q values obtained for glucocorticoid-exposed cells from ASM derived from donors with fatal asthma versus donors without asthma. Genes with a consistent direction of effect sizes are in red and those with opposite directionality are in yellow.

Table 2.

Top Differentially Expressed Genes in ASM Derived from Donors without Asthma that Were Exposed to Budesonide vs. Control

Gene Symbol	Gene Name	Ensembl ID	Genomic Coordinates	Log2 Fold Change	q Value	Mean Normalized Counts
						Budesonide	Control
SAMHD1	SAM domain and HD domain 1	ENSG00000101347	chr20:36890229-36951843	3.84	2.36E-295	85,030.38	5,898.05
GLUL	Glutamate-ammonia ligase	ENSG00000135821	chr1:182381704-182392206	3.41	2.36E-295	141,585.58	13,262.35
NKD1	Naked cuticle homolog 1 (Drosophila)	ENSG00000140807	chr16:50548330-50640739	3.19	2.36E-295	2,308.45	258.82
TSC22D3	TSC22 domain family member 3	ENSG00000157514	chrX:107713221-107777342	4.69	2.36E-295	27,896.26	1,082.04
FKBP5	FK506 binding protein 5	ENSG00000096060	chr6:35573585-35728583	4.85	1.07E-266	16,504.59	569.54
FAM107A	Family with sequence similarity 107 member A	ENSG00000168309	chr3:58564117-58627610	5.18	2.11E-198	1,723.13	51.13
SLC6A6	Solute carrier family 6 (neurotransmitter transporter), member 6	ENSG00000131389	chr3:14402576-14489349	−3.78	5.53E-159	410.99	5,659.88
WNT2	Wnt member 2	ENSG00000105989	chr7:117276631-117323289	−3.18	2.69E-136	183.45	1,798.26
ITGA10	Integrin subunit α 10	ENSG00000143127	chr1:145891208-145910189	4.08	4.76E-123	1,642.39	98.60
ANGPTL7	Angiopoietin-like 7	ENSG00000171819	chr1:11189341-11195981	6.00	9.62E-119	2,873.14	63.73
FAM196A	Family with sequence similarity 196 member A	ENSG00000188916	chr10:127135426-127196158	5.42	7.71E-118	1,145.58	28.05
TNFRSF11B	TNF receptor superfamily member 11b	ENSG00000164761	chr8:118923557-118952200	−3.30	8.00E-115	4,133.93	40,170.43
ADRA2A	Adrenoceptor α 2A	ENSG00000150594	chr10:111077163-111080907	−4.56	1.55E-101	55.60	1,172.94
VCAM1	VCAM 1	ENSG00000162692	chr1:100719742-100739045	−3.97	3.54E-101	189.64	2,807.40
GPX3	Glutathione peroxidase 3	ENSG00000211445	chr5:151020438-151028993	4.33	8.91E-101	92,436.17	4,971.85
TSPAN8	Tetraspanin 8	ENSG00000127324	chr12:71125085-71441898	4.97	2.45E-99	1,003.49	36.88
FMO1	Flavin-containing monooxygenase 1	ENSG00000010932	chr1:171248471-171285978	−3.13	3.21E-95	34.54	290.13
CRISPLD2	Cysteine-rich secretory protein LCCL domain containing 2	ENSG00000103196	chr16:84819984-84920768	3.82	7.83E-95	12,878.40	965.89
RTKN2	Rhotekin 2	ENSG00000182010	chr10:62183035-62268707	−3.36	2.85E-93	60.17	641.42
MOB3B	MOB kinase activator 3B	ENSG00000120162	chr9:27325209-27529781	3.84	2.68E-89	805.24	60.61
PDK4	Pyruvate dehydrogenase kinase, isozyme 4	ENSG00000004799	chr7:95583499-95596491	3.68	5.29E-83	32,731.61	2,553.39
CORO6	Coronin 6	ENSG00000167549	chr17:29614756-29622907	3.01	1.59E-80	2,561.57	314.22
CPA4	Carboxypeptidase A4	ENSG00000128510	chr7:130293134-130324180	−3.19	9.44E-80	23.75	215.28
LRRTM2	Leucine-rich repeat transmembrane neuronal 2	ENSG00000146006	chr5:138868923-138875368	−4.94	1.05E-74	7.72	241.24
MAP1LC3C	Microtubule-associated protein 1 light chain 3 γ	ENSG00000197769	chr1:241995490-241999073	3.08	4.27E-72	1,617.69	202.76
PER1	Period circadian clock 1	ENSG00000179094	chr17:8140472-8156506	3.19	3.84E-67	2,033.90	214.94
SLC2A5	Solute carrier family 2 (facilitated glucose/fructose transporter), member 5	ENSG00000142583	chr1:9035107-9088478	−3.52	1.37E-62	26.09	295.36
KLF15	Kruppel-like factor 15	ENSG00000163884	chr3:126342635-126357442	4.00	1.17E-58	3,186.94	211.51
LINC01088	Long intergenic non-protein-coding RNA 1088	ENSG00000249307	chr4:78971748-79308798	4.24	2.57E-56	259.70	12.71
MARCH10	Membrane-associated ring-CH-type finger 10	ENSG00000173838	chr17:62701314-62808344	4.97	2.66E-55	365.12	12.42
SLC14A1	Solute carrier family 14 (urea transporter), member 1 (Kidd blood group)	ENSG00000141469	chr18:45724127-45752520	−5.30	1.95E-54	4.25	182.86
TIMP4	TIMP metallopeptidase inhibitor 4	ENSG00000157150	chr3:12153051-12159351	3.57	4.14E-51	597.41	53.43

A total of 32 differentially expressed genes were selected based on having a q value < 10−50 and an absolute log2 fold change > 3.

Table 3.

Top Differentially Expressed Genes in ASM Derived from Donors with Fatal Asthma that Were Exposed to Budesonide versus Control

Gene Symbol	Description Gene Name	Ensembl ID	Genomic Coordinates	Log2 Fold Change	q Value	Mean Normalized Counts
						Budesonide	Control
FKBP5	FK506-binding protein 5	ENSG00000096060	chr6:35573585-35728583	4.94	2.92E-212	13,215.86	448.56
GLUL	Glutamate-ammonia ligase	ENSG00000135821	chr1:182381704-182392206	3.24	2.92E-212	99,440.24	10,315.31
TSC22D3	TSC22 domain family member 3	ENSG00000157514	chrX:107713221-107777342	4.43	2.92E-212	23,256.87	1,097.98
CRISPLD2	Cysteine-rich secretory protein LCCL domain containing 2	ENSG00000103196	chr16:84819984-84920768	3.62	3.50E-197	10,835.92	944.47
FAM107A	Family with sequence similarity 107 member A	ENSG00000168309	chr3:58564117-58627610	5.23	8.79E-151	1,834.58	56.25
SPARCL1	SPARC-like 1	ENSG00000152583	chr4:87473335-87531061	4.02	2.33E-133	17,328.36	1,241.11
TNFRSF11B	TNF receptor superfamily member 11b	ENSG00000164761	chr8:118923557-118952200	−3.43	1.43E-131	3,481.06	36,916.42
SAMHD1	SAM domain and HD domain 1	ENSG00000101347	chr20:36890229-36951843	3.73	2.89E-131	73,971.45	5,077.30
NKD1	Naked cuticle homolog 1 (Drosophila)	ENSG00000140807	chr16:50548330-50640739	3.02	3.05E-128	1,883.08	225.67
ZBTB16	Zinc finger and BTB domain containing 16	ENSG00000109906	chr11:114059593-114250676	6.75	1.48E-107	1,352.16	15.41
VCAM1	VCAM 1	ENSG00000162692	chr1:100719742-100739045	−3.97	4.03E-102	302.66	4,219.32
CORO6	Coronin 6	ENSG00000167549	chr17:29614756-29622907	3.41	1.00E-99	2,057.03	201.94
WNT2	Wnt member 2	ENSG00000105989	chr7:117276631-117323289	−3.20	2.69E-98	150.32	1,470.78
ALOX15B	Arachidonate 15-lipoxygenase, type B	ENSG00000179593	chr17:8039017-8049134	7.95	2.87E-96	537.57	2.49
ITGA10	Integrin subunit α 10	ENSG00000143127	chr1:145891208-145910189	4.43	4.35E-95	1,400.29	65.00
FAM196A	Family with sequence similarity 196 member A	ENSG00000188916	chr10:127135426-127196158	5.59	3.15E-92	1,030.13	20.33
KLF15	Kruppel-like factor 15	ENSG00000163884	chr3:126342635-126357442	3.93	2.39E-84	2,492.04	187.84
ADRA2A	Adrenoceptor α 2A	ENSG00000150594	chr10:111077163-111080907	−4.60	1.21E-78	22.81	563.49
TSPAN8	Tetraspanin 8	ENSG00000127324	chr12:71125085-71441898	4.86	9.95E-71	888.70	35.63
RTKN2	Rhotekin 2	ENSG00000182010	chr10:62183035-62268707	−3.29	1.87E-66	81.20	805.66
SLC6A6	Solute carrier family 6 (neurotransmitter transporter), member 6	ENSG00000131389	chr3:14402576-14489349	−3.31	2.39E-64	670.35	5,708.26
PDK4	Pyruvate dehydrogenase kinase, isozyme 4	ENSG00000004799	chr7:95583499-95596491	3.68	1.22E-59	24,018.27	1,859.39
MOB3B	MOB kinase activator 3B	ENSG00000120162	chr9:27325209-27529781	3.87	1.95E-56	569.75	53.36
LINC01088	Long intergenic non-protein-coding RNA 1088	ENSG00000249307	chr4:78971748-79308798	4.05	1.95E-54	178.61	10.58
LRRTM2	Leucine-rich repeat transmembrane neuronal 2	ENSG00000146006	chr5:138868923-138875368	−4.77	5.92E-54	5.82	157.80

A total of 25 differentially expressed genes were selected based on having a q value < 10−50 and an absolute log2 fold change > 3.

To gain a sense of the biological pathways represented by the large number of genes that were differentially expressed in ASM exposed to budesonide, we performed a gene set enrichment analysis of results for 1) ASM from donors with fatal asthma and 2) ASM from donors without asthma. Seventeen ontological categories were overrepresented in fatal asthma, and 14 were overrepresented in nonasthma, based on having a q value < 0.05. These enriched categories were similar in both fatal asthma and nonasthma, according to the enrichment scores and direction of the genes’ differential expression (Figure 2). Twelve of 17 significant ontological categories in fatal asthma, and nine of 14 significant ontological categories in nonasthma were enriched by genes with decreased expression in response to glucocorticoid exposure (i.e., had a normalized enrichment score (NES) < 0), including categories related to cytokines, chemokines, and nitric oxide. Pathways enriched by genes with increased expression in response to glucocorticoid exposure (i.e., with NES > 0) included metabolic pathways of lipids and lipoproteins and amino acid metabolism. Complete lists of the genes in the leading-edge subsets of these categories are provided in Table E4.

Figure 2.

Gene set enrichment analysis results corresponding to differential expression results of a BUD versus control comparison of ASM derived from donors with fatal asthma and donors without asthma. Pathways are sorted by normalized enrichment scores from results for ASM from donors with fatal asthma, and those with a q value < 0.05 are marked with an asterisk (*). GPCR = G protein–coupled receptor; KEGG = Kyoto Encyclopedia of Genes and Genomes; MHC = major histocompatibility complex; TGF = transforming growth factor

Global Transcriptome Response to Glucocorticoid Exposure across Multiple Cell Types

We analyzed 13 publicly available transcriptomic studies of various cell types that compared glucocorticoid versus control exposure in vitro (Table E5). For cell types with more than one available study, we obtained integrated differential expression results. Specifically, integration results were obtained for four ASM studies (including the results for ASM derived from donors without asthma presented above) consisting of 18 pairs of donors, three childhood acute lymphoblastic leukemia (chALL) studies consisting of 12 pairs of donors, and three macrophage studies consisting of 13 pairs of donors. These integrated analyses identified 154, 8, and 78 significant glucocorticoid-responsive genes in ASM, chALL, and macrophages, respectively (Tables E6–E8).

A single study of MCF10A-Myc (immortalized human mammary epithelial cell line MCF10A overexpressing c-Myc) cells showed no statistically significant genes, and a bronchial epithelium cell study showed two significant genes (CAPN15 and CCL20). A comparison of significant genes across all other cell types (Figure E3) revealed that only two well-known glucocorticoid-induced genes, TSC22D3 and FKBP5, had statistically significant changes in ASM, chALL, macrophages, U2OS (human bone osteosarcoma epithelial cell line) cells, and lymphoblastoid cell lines (LCLs), whereas other previously reported glucocorticoid-responsive genes differed according to cell type (GLUL and PER1 had increased expression in ASM and macrophages; CRISPLD2 and KLF15 had increased expression in ASM and U2OS cells) (Figure 3).

Figure 3.

Cell-type–specific differential expression results for known glucocorticoid-responsive genes (TSC22D3, FKBP5, GLUL, PER1, CRISPLD2, and KLF15). Differentially expressed genes with q-value < 0.05 in corresponding tissues are marked with an asterisk (*). BE = bronchial epithelium; chALL = childhood acute lymphoblastic leukemia; CRISPLD2 = cysteine-rich secretory protein LCCL domain containing 2; FKBP5 = FK506 binding protein 51; GLUL = glutamate-ammonia ligase; KLF15 = Kruppel-like factor 15; LCL = lymphoblastoid cell line; MCF10A-Myc = immortalized human mammary epithelial cell line MCF10A overexpressing c-Myc; PER1 = period circadian clock 1; TSC22D3 = TSC22 domain family member 3; U2OS = human bone osteosarcoma epithelial cell line.

ASM-Specific Transcriptome Response to Glucocorticoid Exposure

We reasoned that the set of genes that were responsive to glucocorticoid exposure across four ASM datasets would be more generalizable than that observed in any individual ASM study, and therefore we focused on the set of 56 genes that were significantly differentially expressed in response to glucocorticoid exposure in nonasthma ASM and also had results available in the glucocorticoid-response datasets of other cell types to create an ASM-specific gene expression signature. Hierarchical clustering using the log2 fold change in expression for these 56 genes tended to be similar according to cell type (Figure 4). Of note, the five ASM comparisons were distinct from other cell types, including the cells from donors with fatal asthma that were not used to identify the glucocorticoid-response signature. The three chALL and three macrophage studies clustered according to cell type, with LCLs most closely resembling chALL cells, cell types with shared origin.

Figure 4.

Heatmap comparing ASM-specific differentially expressed genes across seven cell types. A total of 56 significantly differentially expressed genes were identified from the integration of four nonasthma ASM studies that also had results available in the 13 publicly available datasets considered. Dendrograms correspond to hierarchical clustering based on the Euclidean distance between the log2 fold change in expression.

Although overall glucocorticoid-induced gene expression changes were similar in ASM derived from donors with fatal asthma and donors without asthma, we searched for the greatest differences in effect sizes between cells from donors with fatal asthma and donors without asthma among the 154 genes obtained via an integrated analysis of all nonasthma ASM studies. The greatest difference was observed in CEBPD, which had a log2 fold change of 0.48 in ASM from donors with fatal asthma versus a log2 fold change of 0.93 in ASM derived from donors without asthma in the present study or 1.43 across the four integrated nonasthma ASM studies (Figure E4 and Table E9).

Discussion

Glucocorticoids are a mainstay of asthma therapy that modulate the transcription of genes in a tissue-dependent fashion (18–23). Previous studies have sought to elucidate their mechanisms of action in prominent asthma tissues, including ASM (15, 24, 25). Because primary ASM cells derived from individuals with fatal asthma retain a unique phenotype that exhibits differences in proliferative and contractile outcomes compared with ASM cells from individuals without asthma, we hypothesized that the transcriptomic response to glucocorticoids would differ between ASM cells from donors with fatal asthma and those obtained from donors without asthma. Although we found strong gene expression changes in ASM cells due to glucocorticoid exposure, the changes observed were similar in cells derived from donors with fatal asthma and donors without asthma in terms of both statistical significance and direction of effect. Consistent with this, a comparison of the transcriptomes of ASM cells from donors with fatal asthma versus donors without asthma under control exposure yielded no statistically significant gene expression changes, and under budesonide exposure, only two differentially expressed genes (CCK and PMEL) were observed. CCK and its receptor (CCKAR), which are known to influence gallbladder contraction (40), were recently implicated in ASM contraction and obesity-induced airway hyperresponsiveness in mice (41), lending further support for the involvement of CCK in asthma. PMEL is known to play a role in pigmentation and has been linked to amyloid formation, with less evidence for a direct role in asthma (42). Further experimental studies are needed to verify the potential role of these two genes in modulating glucocorticoid response in persons with asthma.

Several of the glucocorticoid-responsive genes identified were consistent with those observed in previous studies, including FKBP5 (35), GLUL (36), PER1 (37, 38), TSC22D3 (35, 39), CRISPLD2 (15), and KLF15 (24, 27). Due in part to the larger sample size in the current study compared with previous ones, we identified many more statistically significant genes than previous ASM glucocorticoid-response transcriptomic studies (15, 24, 25). The large number of genes identified underscores the difficulty of understanding glucocorticoid responsiveness via mechanistic studies of individual genes. To identify biological functions represented by the more than 6,000 genes that were differentially expressed in response to glucocorticoid exposure, an ontological enrichment analysis was performed with the gene set enrichment algorithm, a method that does not require filtering of genes to a small subset, but rather ranks all genes within a dataset. Consistent with the similarities at the level of differentially expressed genes, many ontological categories were similarly enriched in ASM cells derived from both donors with fatal asthma and donors without asthma that were exposed to glucocorticoids. For example, categories that were significant in fatal asthma and nonasthma included Reactome: chemokine receptors bind chemokines driven by decreased expression of proinflammatory genes (e.g., CCL2, CX3CL1, CCL11) and KEGG: cytokine–cytokine receptor interaction driven by decreased expression of many cytokine and TNF-related genes, as well as the Reactome: nitric oxide stimulates guanylate cyclase category, consistent with the known role of glucocorticoids in decreasing inflammation (8) and influencing the nitric oxide pathway (43, 44). Although the NES values were similar for categories enriched among the differentially expressed genes in both fatal asthma and nonasthma, there were categories that were only statistically significant in one group. For example, there was negative enrichment of the categories Reactome: signaling by GPCR, Reactome: G alpha (i) signaling events, and Reactome: voltage gated potassium channels among the fatal asthma, but not the nonasthma, results. Whether these shifts between statistically significant categories represent biologically meaningful changes in gene expression remains to be determined by additional functional studies of genes in these pathways.

Limitations of our study could have hindered our ability to detect statistically significant gene expression changes between transcriptomes of ASM derived from donors with fatal asthma and donors without asthma either when exposed to vehicle control or budesonide. First, although our sample size was large enough for us to detect the strong effects of glucocorticoid exposure, detecting interindividual variability among donors with fatal asthma may require a considerably larger sample size than nine patients and eight control subjects. As previous studies of whole-blood and bronchial epithelium expression signatures in people with asthma have shown, asthma endotypes can be detected among hundreds of subjects (45, 46). Second, because we lacked clinical data on donors with fatal asthma beyond what was provided in medical examiners’ reports, we were unable to establish some important characteristics for each subject. For example, the only medication information available was provided by families to medical examiners in response to a query regarding known medications the deceased individuals were taking near the time of death. Seven of the nine donors with fatal asthma were reported to have taken albuterol, one took fluticasone/salmeterol, and one took unspecified asthma medications. Although we assumed that these medication reports were accurate, we were unable to determine the completeness, dosages, or history of medication use. Thus, even though the medical examiners’ reports are likely to be accurate in stating that asthma was the cause of death, and previous studies have shown that ASM cells derived from donors with fatal asthma have different proliferative and contractility characteristics compared with those obtained from donors without asthma, there was a clinical heterogeneity among the subjects beyond the demographic characteristics used to match patients and control subjects that we were unable to account for. Studies involving a larger number of subjects, as well as subjects with confirmed severe, steroid-resistant asthma, would be helpful to verify that ASM gene transcription changes are similar in people with and without asthma.

Considerations regarding the experimental design that may have influenced our results include 1) the provenance of the ASM and 2) additional factors that influence the glucocorticoid response at a cellular level. The ASM cells used in the present study were obtained from tracheae and large bronchi, and thus our results may not reflect the transcriptomic glucocorticoid response in distal airways. Glucocorticoid receptor number and glucocorticoid receptor gene isoforms influence glucocorticoid sensitivity (23, 47, 48). Because we did not measure glucocorticoid receptor number or characterize the levels of α and β receptor isoforms in individual donors, our results reflect only the overall downstream transcriptomic changes. Future studies to address these issues may shed further light on the mechanisms underlying the differential glucocorticoid response between persons with and without asthma.

We compared glucocorticoid-response gene expression signatures obtained in four studies of nonasthma ASM with those obtained in six other cell types using 13 publicly available glucocorticoid-response studies. Among the four ASM studies consisting of a total of 18 pairs of cells derived from donors without asthma, the present study contributed the largest sample size (n = 8 donors), and thus the overall results for ASM are driven by the current results. Nonetheless, selection of genes based on all available studies was likely to select a generalizable expression signature. The hierarchical clustering results based on the 56-gene signature showed that there were cell-type–specific gene expression patterns across various cell types, with some expected similarities (e.g., all ASM studies clustered together, and chALL and LCL clustered together). Surprisingly, only two genes (TSC22D3 and FKBP5) had statistically significant results across all cell types. Several other well-known glucocorticoid-responsive genes were selectively expressed across cell types in a differential manner. Although the integrated results were biased according to the current availability of publicly available transcriptomic data, which do not evenly represent all cell and tissue types or involve treatment with the same type of glucocorticoid or time frame, integration of publicly available data offers a cost-effective avenue to identify the most robust gene signatures possible, as well as cell-type–specific ones. Future efforts to perform integrative analyses with an increased number of datasets may reveal more cell-type–specific signatures of glucocorticoid responsiveness. Additionally, such signatures may be linked to glucocorticoid responsiveness among patients with asthma. For example, we found that CEBPD was differentially expressed in response to glucocorticoid exposure in ASM derived from both donors with fatal asthma and donors without asthma in the current study (q = 2.71 × 10−5 and 0.032, respectively), but its overall ranking did not make it a high priority for further study. The evidence provided by integrated studies (Figure E4 and Table E9) and published functional evidence for CEBPD (Table E10), however, suggest that it should be prioritized for further study. Thus, integration results offer a complementary approach to prioritize individual genes for functional studies.

In summary, we have expanded studies of the glucocorticoid response in ASM to include donors with fatal asthma and a greater number of subjects. We identified more differentially expressed genes than previous studies and an ASM-specific expression signature that distinguished the ASM glucocorticoid response from that of six other cell types. The similarity of glucocorticoid-response expression between ASM cells derived from donors with fatal asthma and donors without asthma suggests that differences in glucocorticoid responsiveness among asthma patients are not mediated by striking changes in ASM transcript expression signatures. Further studies are necessary to confirm this observation, determine whether it applies to ASM obtained from distal airways, and identify other mechanisms that may confer differences in glucocorticoid response, such as poor delivery of glucocorticoids to target tissues.