Vitamin D Modulates Expression of the Airway Smooth Muscle Transcriptome in Fatal Asthma

Globally, asthma is a chronic inflammatory respiratory disease affecting over 300 million people. Some asthma patients remain poorly controlled by conventional therapies and experience more life-threatening exacerbations. Vitamin D, as an adjunct therapy, may improve disease control in severe asthma patients since vitamin D enhances glucocorticoid responsiveness and mitigates airway smooth muscle (ASM) hyperplasia. We sought to characterize differences in transcriptome responsiveness to vitamin D between fatal asthma- and non-asthma-derived ASM by using RNA-Seq to measure ASM transcript expression in five donors with fatal asthma and ten non-asthma-derived donors at baseline and with vitamin D treatment. Based on a Benjamini-Hochberg corrected p-value <0.05, 838 genes were differentially expressed in fatal asthma vs. non-asthma-derived ASM at baseline, and vitamin D treatment compared to baseline conditions induced differential expression of 711 and 867 genes in fatal asthma- and non-asthma-derived ASM, respectively. Functional gene categories that were highly represented in all groups included extracellular matrix, and responses to steroid hormone stimuli and wounding. Genes differentially expressed by vitamin D also included cytokine and chemokine activity categories. Follow-up qPCR and individual analyte ELISA experiments were conducted for four cytokines (i.e. CCL2, CCL13, CXCL12, IL8) to measure TNFα-induced changes by asthma status and vitamin D treatment. Vitamin D inhibited TNFα-induced IL8 protein secretion levels to a comparable degree in fatal asthma- and non-asthma-derived ASM even though IL8 had significantly higher baseline levels in fatal asthma-derived ASM. Our findings identify vitamin D-specific gene targets and provide transcriptomic data to explore differences in the ASM of fatal asthma- and non-asthma-derived donors.


Introduction
Asthma, a chronic inflammatory respiratory disease, is defined by increased airway responsiveness and inflammation to specific environmental stimuli [1]. Clinical therapy following established guidelines successfully controls asthma symptoms in most patients [2]. Such therapy includes use of glucocorticoids, anti-inflammatory medications [3], and β 2 -agonists that promote bronchodilation [4]. Up to 5% of asthma patients, however, are estimated to inadequately respond to conventional approaches, and consequently, are at greater risk of suffering lifethreatening asthma exacerbations [5]. While patients with severe, refractory disease represent a heterogeneous group, features shared by most are glucocorticoid insensitivity and, in part, irreversible airflow obstruction [6].
Low vitamin D levels are associated with decreased lung function and increased asthma exacerbations and medication use, although it is unclear whether there is a causal relationship among these observations or whether asthma patients can benefit from vitamin D supplementation [7]. Mechanisms by which vitamin D may ameliorate the asthma diathesis include modulating inflammation by increasing CD4+ T cell synthesis of IL-10 [8], suppressing Th17 cytokines in PBMCs [9], and increasing the number of regulatory T cells [10]. These antiinflammatory effects of vitamin D extend to structural cells, including the airway smooth muscle (ASM). For example, vitamin D reduced TNFα-induced CCL5 (RANTES) and CXCL10 (IP-10) secretion in ASM cells, while secretion of CX3CL1 (fractalkine), which is glucocorticoid-insensitive, was also decreased by vitamin D treatment [11]. Others found that vitamin D decreased expression of KPNA4 (importin α3), a signaling molecule for NFκB-induced inflammation in the ASM [12]. Vitamin D anti-inflammatory effects are, in part, mediated by enhancing the action of glucocorticoids, suggesting that vitamin D supplementation may benefit patients with glucocorticoid-insensitive asthma [7].
Airway remodeling, an important feature of severe persistent asthma that includes ASM hyperplasia, is modulated by vitamin D since growth factor-induced ASM proliferation was inhibited by vitamin D treatment [13,14]. Children with severe, therapy-resistant asthma also manifested lower levels of vitamin D than children with moderate asthma; endobronchial biopsies showed an inverse relationship of vitamin D levels and ASM mass, but not epithelial shedding, reticular basement membrane thickness, or airway inflammation measures [15]. Accordingly, the beneficial vitamin D effects on the ASM may extend beyond the enhancement of glucocorticoid actions since glucocorticoids have little effect on airway remodeling or ASM proliferation [11].
Vitamin D (1,25(OH) 2 D 3 ) plays a pivotal role in the regulation of gene transcription [16]. Consequently, the study of vitamin D effects on gene transcription in multiple cell types may reveal changes in signaling pathways that are relevant to specific diseases. For example, a microarray study that used vitamin D-treated ASM cells from a single male donor without asthma identified transcription changes in gene categories that might influence airway remodeling, including morphogenesis, cell growth and cell survival [17]. Recent advances in sequencing technologies provide comprehensive and in-depth characterization of transcriptomes via RNA-Seq [18]. Compared to the use of microarrays, RNA-Seq offers advantages by examining greater numbers of RNA species over a wider dynamic range of signal [19].
In this study, we used RNA-Seq to identify differences between the transcriptomes of ASM derived from fatal asthma and non-asthma donors to vitamin D treatment, and we also specifically examined four cytokines (i.e. CCL2, CCL13, CXCL12, IL8) whose expression differed between fatal asthma-and non-asthma-derived ASM at baseline and with vitamin D treatment. Table 1. Characteristics of ASM donors. All donors were white non-smokers. There were no significant differences between sex, age, or BMI of fatal asthma-vs. non-asthma-derived donors. Medical examiner ruled cause of death for fatal asthma donors was "Asthma Attack/Anoxia," or a significant asthma event is listed as preceding death.
Fatal Asthma (N = 5) Non-Asthma (N = 10)   A subset of 307 genes was vitamin D responsive in both fatal asthma-and non-asthma-derived ASM. As for the fatal asthma-vs. non-asthma-derived comparison, differentially expressed genes were prioritized by (1) selecting the genes with the lowest q-value of 1.92E-03 (n = 394 for non-asthma-derived ASM; n = 345 for fatal asthma-derived ASM), (2) selecting the subset of these genes with magnitude of Log 2 fold-change >1.3 (n = 105 for nonasthma-derived ASM; n = 214 for fatal asthma-derived ASM), and (3) dropping genes for which at least two individual samples did not have FPKM values >5 [Tables 3 and 4].
To gain a better sense of biological functions represented by the differentially expressed genes, we performed an ontological category enrichment analysis using the NIH DAVID tool [29]. For the non-asthma-derived ASM at baseline vs. vitamin D treated comparison, there with 24 annotation clusters with enrichment scores >2.50 [S3 Table]. For the fatal asthmaderived ASM at baseline vs. vitamin D treated comparison, there with 21 annotation clusters with enrichment scores >2.50 [S4 Table]. Some of the individual asthma-related enrichment categories within clusters with Benjamini-Hochberg adjusted p-values < 0.05 were also present in the fatal asthma-vs. non-asthma-derived comparison: extracellular matrix, immunoglobulin domain, response to steroid hormone stimulus, and response to wounding [S5 and S6 Tables], although some genes within categories differed. The oxidative stress and lung development categories were still significant for the non-asthma-derived ASM treated with vitamin D vs. baseline but not for the fatal asthma-derived ASM treated with vitamin D vs. baseline. Two asthma-related categories that were not present in the fatal asthma-vs. non-asthma-derived ASM comparison were significant in both the vitamin D treated vs. baseline comparisons: cytokine activity and chemokine activity [S5 and S6 Tables]. TNFα Responsiveness of Select Cytokines Measured by q-PCR Four genes (i.e. CCL2, CCL13, CXCL12, IL8) that (1) were members of the GO:0008009~chemokine activity ontological category, which was significantly over-represented among genes differentially expressed in response to vitamin D treatment in both fatal asthma-and nonasthma-derived ASM [S5 and S6 Tables], and (2) were differentially expressed in fatal asthma vs. non-asthma-derived ASM at baseline [S2 Table] were selected for detailed examination [Fig 1; Table 5]. Differential expression for these four genes at baseline and with vitamin D treatment was measured via qRT-PCR. No statistically significant differences in mRNA expression between fatal asthma-and non-asthma-derived ASM at baseline were found (data not shown). Fold changes for CCL2 and CXCL12 expression with vitamin D treatment vs. baseline were consistent with RNA-Seq results for most samples, while CCL13 and IL8 had consistent effects for only some samples [ Fig 2]. We found that TNFα treatment induced CCL2, CCL13, and IL8 mRNA expression in most fatal asthma-and non-asthma-derived ASM samples, but TNFα treatment did not induce CXCL12 [Fig 2]. Vitamin D decreased the mean TNFαinduced CCL2, CCL13, and IL8 mRNA expression in most ASM groups, but this change was not statistically significant or consistent across all samples [ Fig 2].

Characterization of Cytokine Secretion and Sensitivity to Vitamin D
Our RNA-Seq results suggested that mRNA for four cytokines were differentially expressed in ASM derived from fatal asthma vs. non-asthma donors at baseline and with vitamin D treatment. To address whether these differences in mRNA extended to changes in protein secretion, individual analyte ELISAs for the four selected cytokines (i.e. CCL2 (aka MCP-1), CCL13 (aka MCP-4), CXCL12 (aka SDF-1) and IL8) were performed. At baseline, CCL2 and CXCL12 were    secreted differently from fatal asthma and non-asthma derived ASM in directions consistent with RNA-Seq results, but only those of CXCL12 were statistically significant [Fig 3]. IL8 secretion was not detected at baseline, while that of CCL13 was low and not different between fatal asthma and non-asthma derived ASM [Fig 3]. Fatal asthma-derived ASM secreted greater TNFα-induced IL8 and lower TNFα-induced CXCL12 as compared to non-asthma-derived ASM, while TNFα-induced CCL2, and CCL13 levels between fatal asthma-and non-asthmaderived ASM were not significantly different [ Fig 3]. Overall, TNFα-induced secretion of CCL13 was low, while that of CCL2 and IL8 were high. Because IL8 was secreted at high levels that were significantly different between fatal asthma and non-asthma-derived ASM, we tested the dose-dependent effect of vitamin D on its TNFα-induced secretion. Additionally, CCL5 (aka RANTES), CXCL10 (aka IP-10) were also chosen for comparison experiments because these chemokines were known mediators in asthma that are differentially sensitive to vitamin D [11]. As shown in Fig 4, vitamin D inhibited three cytokine levels to a comparable degree in fatal asthma-vs. non-asthma-derived ASM, even for IL8, whose TNFα-induced baseline  Table 5. Differential expression results for all comparison groups corresponding to four cytokine genes that were selected for further study.
Genes were selected based on being (1) members of the GO:0008009~chemokine activity ontological category, which was significantly over-represented among genes differentially expressed in response to vitamin D treatment in both fatal asthma-and non-asthma-derived ASM, that were (2) differentially expressed in fatal asthma-derived vs. non-asthma-derived ASM at baseline. FPKM = fragments per kilobase of transcript per million mapped reads. secretion was significantly higher in fatal asthma-vs. non-asthma-derived ASM. Although vitamin D was an effective inhibitor of TNFα-induced cytokine secretion, none of the TNFαinduced cytokine levels were completely abrogated. Further, the amount of TNFα-induced cytokine inhibition by vitamin D differed among cytokines, with IL8 levels being inhibited the least in fatal asthma-and non-asthma-derived ASM cells.

Discussion
Previous studies have focused on understanding the potential therapeutic role of vitamin D in asthma because lower vitamin D levels are associated with increased asthma exacerbations and anti-inflammatory medication use, and decreased lung function [7,15]. A major mode of vitamin D (1,25(OH) 2 D 3 ) action is to directly alter the expression of genes [16]. Understanding the role of vitamin D in modifying the ASM transcriptome of fatal asthma patients may thus provide insights into the mechanisms by which vitamin D influences severe asthma. Our RNA-Seq results found that a large number of genes are differentially expressed in the ASM of fatal asthma vs. non-asthma donors, and that the responsiveness to vitamin D treatment differs between these two groups. Vitamin D treatment of both fatal asthma-and non-asthmaderived ASM resulted in the differential expression of many cytokines and chemokines, but corresponding ontological categories were not found to be over-represented in the baseline comparison between fatal asthma-and non-asthma-derived ASM. Vitamin D receptor (VDR) mRNA expression levels were not different between fatal asthma-and non-asthma derived ASM [ S3 Fig], suggesting that gene expression differences between these groups in response to vitamin D treatment did not simply reflect vitamin D sensitivity differences proportional to vitamin D receptor availability. We chose to focus our functional studies on four cytokines that were differentially expressed between fatal asthma-and non-asthma-derived ASM and were also modified by vitamin D treatment in fatal asthma-derived ASM (i.e., CCL2, CCL13, CXCL12, IL8). While each of these genes was differentially expressed in fatal asthma-vs. non-asthma-derived ASM at baseline and with vitamin D treatment in fatal asthma-derived ASM, expression variability across samples was evident [Fig 1]. Subsequent qPCR results did not replicate RNA-Seq differential expression between fatal asthma-and non-asthma-derived ASM at baseline, and were only partly consistent with vitamin D effects: consistent with RNA-Seq results, CCL2 had decreased and CXCL12 had increased mRNA with vitamin D treatment for most samples, while CCL13 and IL8 manifested opposite mean effects [Fig 2]. These disparate findings could be due to inter-individual variability, unaccounted for experimental error by either RNA-Seq or qPCR, limitations in the methods used to analyze RNA-Seq data, and/or the intrinsic differences in the experimental techniques. In the case of CCL13 and IL8, mRNA expression levels were low as evidenced by having most RNA-Seq FPKM values <10, which may have contributed toward inconsistency in measures between RNA-Seq and qPCR, as results from lowly expressed genes are known to be more error prone [30]. There is no consensus about the optimal way to analyze RNA-Seq data and considerable variability in results obtained based on platform and data analysis pipeline has been observed [31]. While RNA-Seq measures for each condition were normalized by taking some global measures into account (e.g., the total number of reads that mapped to the reference genome), qPCR measures were obtained in reference to a Although the amount of inhibition differed by cytokine, vitamin D inhibited three cytokine levels to a comparable degree in fatal asthma-derived vs. non-asthma-derived ASM, even for IL8, whose TNFα-induced baseline secretion was significantly higher in fatal asthma-vs. non-asthma-derived ASM. Data represent means ± standard error of the mean for ASM cells derived from 5 fatal asthma donors and 10 non-asthma donors. Each condition was measured in triplicate. Statistical significance was determined using Student's two-tailed t-test with a p<0.05 threshold. Bottom panels compare baseline TNFα-induced cytokine levels to those obtained with a maximal inhibitory vitamin D concentration (I max ). We selected the housekeeping gene with lowest amount of variation to minimize the impact that its variability of expression would have on our ability to measure changes in the expression of other genes of interest. However, because a housekeeping gene with even greater coefficient of variability yielded similar results, the housekeeping gene choice may not have been the main source of inconsistency between RNA-Seq and qPCR results. Further sequencing and qPCR experiments using a greater number of individuals may clarify the source of inconsistency between techniques.
Comparison of RNA and protein expression for the four cytokines in Figs 1 and 3 shows (1) decreased CXCL12 mRNA, protein secretion at baseline, and TNFα-induced protein secretion in fatal asthma-vs. non-asthma-derived ASM, (2) increased CCL2 mRNA, protein secretion at baseline, and TNFα-induced protein secretion in fatal asthma-vs. non-asthma-derived ASM, although the protein secretion results are not statistically significant, (3) increased IL8 mRNA and TNFα-induced protein secretion in fatal asthma-vs. non-asthma-derived ASM, but undetectable levels of protein secretion at baseline, and (4) increased CCL13 mRNA in fatal asthmavs. non-asthma-derived ASM, but unchanged protein secretion levels. Comparison of TNFαinduced mRNA and protein expression levels of for the four cytokines in Figs 2 and 3 shows that CCL2, CCL13, and IL8 had increased levels in both fatal asthma-and non-asthma-derived ASM, while CXCL12 had decreased levels in all but a few qPCR samples. Thus, consistency between mRNA and protein secretion levels is high for TNFα-induced changes, and present to varying degrees among the four cytokines for changes between fatal asthma-and non-asthmaderived ASM, with CXCL12 and IL8 showing the greatest consistency and CCL13 the weakest. While sometimes the levels of mRNA correlate with those of protein secretion, there are various stages between transcription and protein secretion that could bar a simple relationship, including protein translation, post-translational modifications that activate proteins and/or avoid degradation, and processes that successfully export protein from cells [32,33]. Based on our results, CCL2, CXCL12, and IL8 may have correlated mRNA and protein secretion levels that differ between fatal asthma-and non-asthma-derived ASM, but the difference in IL8 may only be noticeable with TNFα induction due to its low baseline levels. The mRNA differences of CCL13 between fatal asthma-and non-asthma-derived ASM did not extend to protein secretion levels. Further studies that address in more depth the mechanisms by which cytokine mRNA changes lead to protein level and secretion differences are needed to support our observed correlations for CCL2, CXCL12, and IL8.
The CCL2, CCL13, CXCL12, and IL8 cytokine genes selected for experimental follow-up are known to be involved in asthma related processes. A previous study found that CCL2 (aka MCP1) concentration was elevated in primary ASM supernatant and blood of asthma patients, its sputum levels were increased in asthma patients with bronchial wall thickening, and chemotaxis assays suggested that it mediates fibrocyte migration toward ASM [34]. CCL2 has also been found to promote mast cell recruitment [35] and Th17 cell migration [36] to the lung in mouse models of asthma. CCL13 (aka MCP4) is a monocyte chemoattractant that has been proposed as a biomarker of asthma because its levels were higher in plasma of (1) asthma patients vs. non-asthma controls and (2) asthma patients with acute exacerbations vs. stable asthma [37], and its receptor, CCR3, has been shown to induce migration of ASM cells in vitro [38]. CXCL12, a chemotactic factor with pleiotropic effects that are modulated by cofactors, recruits cell types including leukocytes [39] and mast cells [40]. A potential asthma therapeutic based on blocking of CXCL12 has been shown to reduce airway eosinophil recruitment in an ovalbumin mouse model of asthma [41]. Further, airway mucosa from asthma patients had a greater number of CXCL12 (aka SDF1) positive cells than that from healthy subjects, and the number of CXCL12+ cells correlated with vascularity, suggesting that CXCL12 may play a role in airway remodeling [42]. Finally, IL8 is a neutrophil chemoattractant whose levels have been reported to be elevated in the sputum of asthma patients, suggesting that it may play a role in severe asthma that is characterized by neutrophilia [43,44]. Our RNA-Seq results showing that CCL2, CCL13, and IL8 mRNA levels were elevated in fatal asthma-derived ASM are consistent with these previous studies showing elevated levels in ASM (for CCL2) and other tissues. Our RNA-Seq results for CXCL12 showing decreased mRNA levels in fatal asthma have an opposite effect from that expected based on increased CXCL12+ cell types in airway mucosa, but a direct comparison of these studies' results is limited because of the different experimental designs.
Our work extends what is known about these four cytokines by suggesting a potential role of vitamin D in regulating their levels in ASM. However, dysregulation of these cytokines is a complex process with high variability across individuals, and further studies are required to better understand the mechanisms by which vitamin D may modulate their effects.
Changes in expression and secretion of individual genes in ASM cells in response to vitamin D treatment have been studied previously. For example, in a study of the inhibitory properties of vitamin D on ASM proliferation, MMP9 and ADAM33 were selected as candidate genes via which vitamin D might affect ASM hyperplasia and both were differentially expressed with vitamin D treatment [14]. In our results, these genes were not differentially expressed in either of the three comparisons made, but other genes representing the extracellular matrix were significantly represented, including other matrix metalloproteinases (e.g., MMP1, MMP12) and ADAMTS proteases (e.g., ADAMTS10, ADAMTS15) [S5 and S6 Tables]. A previous study of the role of vitamin D in ASM inflammation found that vitamin D decreased expression of KPNA4 (a.k.a. importin α3), a mediator of NFκB-induced inflammation, by TNFα and IL1β in the ASM [12], but in our study, KPNA4 mRNA was not differentially expressed in any of the comparison groups. Another study found that TNFα/IFNγ-induced secretion of CX3CL1 was significantly decreased by vitamin D treatment [11]. Consistent with this result, we found that CX3CL1 expression was decreased in non-asthma-derived ASM treated with vitamin D vs. at baseline (q-value 0.011; log2 fold-change 2.3), but the findings did not hold in fatal asthmaderived ASM. TNFα-induced CCL5 (aka RANTES) and CXCL10 (aka IP-10) chemokine secretion were decreased in ASM cells treated with vitamin D [11], but we did not observe any mRNA expression changes for CCL5 or CXCL10. These differences between protein secretion and transcript expression suggest that some changes in cytokine secretion levels are not transcriptionally mediated, or that stimulation with an appropriate cytokine is necessary to observe the changes.
On a transcriptomic level, the effects of vitamin D on ASM gene expression were measured in a previous microarray study that used cells from a single non-asthmatic male donor stimulated with 100nM 1α,25(OH) 2 D 3 for 24 hours [17]. Our results are consistent with the results reported in that paper. For example, CDC42EP3, CEBPB, G6PD, HSD11B1, IL6, and ITPR1 were highlighted by this previous study as asthma-relevant genes that were differentially expressed in response to vitamin D treatment, and each of these is among the 394 genes with lowest q-value of 1.92E-03 in our baseline vs. vitamin D treated non-asthma-derived ASM comparison [S7 Table]. A subset of these genes (i.e., CDC42EP3, ITPR1, G6PD) was differentially expressed in the fatal asthma-vs. non-asthma-derived ASM at baseline comparison, while a different subset (i.e., CEBPB, CDC42EP3, ITPR1, and HSD11B1) was differentially expressed in fatal asthma-derived ASM when stimulated with vitamin D vs. left untreated [S7 Table]. The differences in significant genes among comparison groups may reflect biological differences in vitamin D responsiveness between fatal asthma-and non-asthma-derived ASM cells. Thus, an important strength of our study compared to previous reports is the inclusion of ASM cells from diseased individuals, which may aid in identifying transcriptomic differences that are most relevant to asthma.
The ontological categories represented by differentially expressed genes in S1, S3, and S4 Tables are consistent with biological processes known to be important in severe asthma. The extracellular matrix is important in ASM remodeling [45], steroid resistance and oxidative stress directly contribute to severe asthma [6], and differences in lung development and response to wounding may contribute to intrinsic differences in ASM that alter asthma severity. Two categories that were not over-represented by differentially expressed genes but might be expected are those related to smooth muscle contraction and smooth muscle proliferation since the ASM regulates airway contractility and ASM proliferation contributes to increased ASM mass [45]. Although statistical significance was not achieved, categories related to these processes characteristic of smooth muscle function were represented by genes significantly differentially expressed in some comparisons. Specifically, GO:0048661~positive regulation of smooth muscle cell proliferation included seven differentially expressed genes (i.e. PRKCA, IL6, S1PR1, PTGS2, VEGFA, IGF1, ITGA2) in non-asthma-derived ASM and six differentially expressed genes (i.e. CDH13, VEGFA, TGM2, IGF1, ITGA2, STAT1) in fatal asthma-derived ASM, with Benjamini-Hochberg corrected p-values of 0.066 and 0.14, respectively. This category was absent in the fatal asthma-vs. non-asthma-derived ASM baseline comparison. The category GO:0006939~smooth muscle contraction included four differentially expressed genes (i.e., ACTG2, EDNRB, AGT, MYH11) in the fatal asthma-vs. non-asthma-derived ASM comparison, six differentially expressed genes (i.e., KCNMA1, EDNRB, CACNA1G, BDKRB2, CAC-NA1C, GDNF) in the non-asthma-derived ASM baseline vs. vitamin D treatment comparison, and four differentially expressed genes (i.e., KCNMA1, EDNRA, ACTG2, AGT) in the fatal asthma-derived ASM baseline vs. vitamin D treatment comparison, with Benjamini-Hochberg corrected p-values of 0.80, 0.22 and 0.69, respectively. Thus, while groups of genes of some relevant ontological categories may not have been differentially expressed, transcriptional changes of individual genes that influence ASM contractility and proliferation were present.
A limitation of our study is that we used ASM cells from white donors to study responsiveness to vitamin D. Dose response rates [46] and vitamin D protein binding levels [47] by race, suggesting that gene expression differences in response to vitamin D treatment also differ by race. Because vitamin D deficiency [48] and asthma prevalence, hospitalizations, and death rates [49] are known to be higher among black than white patients, future studies of the effect of vitamin D on the ASM transcriptome that include individuals of diverse racial/ethnic backgrounds may reveal additional gene candidates and/or biological pathways that shed light on asthma disparities. Another limitation is that our RNA-Seq analyses were limited to the hg19 RefSeq annotation files downloaded from Illumina's iGenomes project. Thus, we did not characterize the expression of long-non-coding RNA or mRNA transcript isoforms that were not part of the reference file used. We opted for use of a well-annotated reference file for our investigation of the ASM transcriptome to reduce the number of false-positive results. Future studies with more comprehensive annotation files, a greater number of individuals, and/or greater sequencing depth will yield additional insight into the ASM transcriptome.
In summary, we found that 838 genes were differentially expressed in fatal asthma-vs. nonasthma-derived ASM at baseline, and vitamin D treatment vs. baseline conditions resulted in differential expression of 867 and 711 genes in non-asthma-and fatal asthma-derived ASM, respectively. Functional categories that were represented in all groups include extracellular matrix, response to steroid hormone stimulus, and response to wounding. Genes differentially expressed by vitamin D also represented cytokine and chemokine activity categories. Further experiments conducted for four cytokines (i.e. CCL2, CCL13, CXCL12, IL8) selected from among the top RNA-Seq results found that vitamin D inhibited TNFα-induced IL8 protein secretion levels to a comparable degree in fatal asthma-and non-asthma-derived ASM even though IL8 had significantly higher baseline levels in fatal asthma-derived ASM. Our results provide transcriptome data to further explore the functional differences in the ASM of fatal asthma vs. non-asthma donors, as well as differences in their vitamin D responsiveness.

Ethics Statement
Lung tissue was obtained from the National Disease Resource Interchange (NDRI) and its use approved by the University of Pennsylvania Institutional Review Board; use of the cells does not constitute human subjects research since all donor tissue is harvested anonymously and de-identified.

ASM Cell Culture and Vitamin D Treatment for RNA-Seq Experiment
Primary ASM cells were isolated from white donors, six who died of fatal asthma and ten with no chronic illness or medication use. ASM cell cultivation was described previously [50,51]. ASM cells up to passage 4 maintained in Ham's F12 medium supplemented with 10% FBS, CaCl 2 , buffered with HEPES, penicillin/streptomycin, primocin, and additional L-glutamine were used in all experiments. The F12 medium was used for culture because it provides Ca 2+ levels that are consistent with seeing contractility of muscles in that media. Following 48 hours of serum deprivation, cells from each donor were treated with 100 nM vitamin D (1,25 (OH) 2 D 3 , Cayman Chemical Company, Ann Arbor, MI) or control vehicle for 18 h. This treatment protocol is biologically relevant (approximately 80μg/dL of vitamin D), maximally induces expression of CYP24A1, a well-known metabolizer of vitamin D, and inhibits cytokine-induced cytokine and chemokine expression in ASM [11,13].

RNA-Seq Library Construction and Sequencing
Total RNA was extracted from cells using the miRNAeasy mini kit (Qiagen Sciences, Inc., Germantown, MD). Approximately 1 μg of RNA from each sample was used to generate RNA-Seq cDNA libraries for sequencing using the TruSeq RNA Sample Prep Kit v2 (Illumina, Inc., San Diego, CA). Sample preparation followed the manufacturer's protocol with a workflow that included isolation of poly-adenylated RNA molecules using poly-T oligo-attached magnetic beads, enzymatic RNA fragmentation, cDNA synthesis, ligation of bar-coded adapters, and PCR amplification. Ambion External RNA Controls Consortium (ERCC) RNA Spike-In Control Mix 1 (Life Technologies Corporation, Carlsbad, CA) was added to the samples. The amplified cDNA fragments were analyzed using the 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA) to determine fragment quality and size.

RNA-Seq Data Analysis
Preliminary processing of raw reads was performed using Casava 1.8 (Illumina, Inc., San Diego, CA). Subsequently, Taffeta scripts (https://github.com/blancahimes/taffeta) were used to analyze RNA-Seq data, which included trimming of adapters using trimmomatic (v.0.22) [52] and using FastQC [53] (v.0.10.0) to obtain overall QC metrics. Trimmed reads for each sample were aligned to the reference hg19 genome and known ERCC transcripts using TopHat [54] (v.2.0.8), while constraining mapped reads to be within reference hg19 or ERCC transcripts. Additional QC parameters were obtained to assess whether reads were appropriately mapped. Bamtools (v.1.0.2) [55] was used to count/summarize the number of mapped reads, including junction spanning reads. The Picard Tools (v.1.58; http://picard.sourceforge.net) RnaSeqMetrics function was used to compute the number of bases assigned to various classes of RNA, according to the hg19 refFlat file available as a UCSC Genome Table. For each sample, Cufflinks [18] (v.2.1.1) was used to quantify ERCC Spike-In and hg19 transcripts based on reads that mapped to the provided hg19 and ERCC reference files. ERCC Spike-ins dose response curves (i.e. plots of ERCC transcript FPKM vs. ERCC transcript molecules) were created following the manufacturer's protocol [56]. Raw read plots were created by displaying bigwig files for each sample in the UCSC Genome Browser. The RNA-Seq data is available at the Gene Expression Omnibus Web site (http://www.ncbi.nlm.nih.gov/geo/) under accession GSE58434.

RNA-Seq Quality Control
Initially, 36 RNA samples from 18 vitamin D treated and untreated pairs were obtained. One pair of RNA samples was excluded from sequencing preparation due to failed RIN scores, indicating very low RNA quality. Two samples with evidence of poor quality RNA-Seq library (based on low percentage of junction spanning reads, percentage of mRNA bases and mean insert size [S8 Table]) and the two paired samples matching them, were also excluded, leaving a total of 30 samples selected for differential expression analysis. Such samples corresponded to 5 fatal asthma and 10 non-asthma-derived donors with vitamin D or control vehicle treatment [ Table 1]. For these 30 samples, we obtained an average of 93.2 million raw sequencing reads per sample (range 73.0-121.7 million reads per sample). Of these reads, an average of 89.1% were aligned to hg19 genome reference files downloaded from Illumina's iGenomes project (range 85.9%-91.6%), and an average of 29.1% of the mapped reads spanned junctions (range 23.2%-31.4%). An average of 92.7% of bases in mapped reads corresponded to mRNA (range 73.6%-96.2%). ERCC spike-in dose response plots, which ideally would have slope and R 2 equal to 1.0, had slope and R 2 > = 0.91. Based on these and other quality control (QC) summary metrics, the sequencing and alignment results for each of the 30 samples were deemed of sufficiently high quality to include in differential expression analyses [S8 Table]. Quantification of transcript and gene expression levels, according to hg19 RefSeq annotation files from Illumina's iGenomes Project, was performed using Cufflinks. As positive controls of gene expression, fragments per kilobase of transcript per million mapped reads (FPKM) values for three housekeeping genes (i.e., B2M, GABARAP, RPL19), one well known metabolizer of vitamin D (i.e., the cytochrome P450 gene CYP24A1 [57] Quantitative Real-Time PCR (qRT-PCR) Analysis Following 48 hours of serum deprivation, ASM cells from 4 fatal and 10 non-asthma-derived donors (all used for RNA-Seq) were treated with 100nM vitamin D or vehicle control for 1h, followed by stimulation with 10ng/mL TNFα or vehicle control for 18h. Total RNA was isolated from cells by using QIAshredder and RNeasy kits (Qiagen Sciences, Inc., Germantown, MD). Oligo(dT)-primed cDNA was prepared from 500 ng of total RNA by using SuperScript III First-strand Synthesis System (Invitrogen, Life Technologies, Grand Island, NY). qRT-PCR was set up in the presence of 0.5 μM primers for CCL2, CCL13, CXCL12, and IL8 by using QuantiTect SYBR Green PCR kit (Qiagen Sciences, Inc., Germantown, MD). qRT-PCR was performed on an Eppendorf Mastercycler ep Realplex 2 real time PCR machine (Eppendorf, Hauppauge, NY). GABARAP was used as an internal control for data normalization as it was the housekeeping gene with the lowest amount of variation among six ASM samples (3 fatal asthma-derived, 3 non-asthma-derived) for four housekeeping genes tested (i.e.

Cytokine ELISAs
To test whether some of our RNA-Seq results extended to protein secretion levels, ELISAs were performed for four cytokines. Human ASM cells were grown to confluence and then serum deprived for 48 h. Monolayers were stimulated with 10 ng/ml TNFα (Roche Diagnostics Corporation, Indianapolis, IN), in the presence and absence of 100 nM vitamin D, or vehicle control. After an 18h incubation, supernatants were harvested and cytokine or chemokine levels measured by single analyte ELISAs (R&D Systems, Minneapolis, MN). Comparisons were made of baseline and TNFα-induced chemokine and cytokine levels in fatal asthma-derived ASM and those derived from non-asthma age-and gender-matched donors. Analyses were performed using ASM cells from 3 fatal asthma and 3 non-asthma donors for CCL2 and CXCL12; 5 fatal asthma and 6 non-asthma donors for CCL13; and 5 fatal asthma and 5 nonasthma donors for IL8 based on availability of cells used after RNA-Seq experiments. Each condition was performed in triplicate. Data was analyzed as means ± standard error of the means and statistics determined significance using one-tailed student's t test.
Supporting Information S1 Fig. RNA-Seq profiling of ASM cells. Volcano plots of overall gene-based differential expression results for A) fatal asthma-vs. non-asthma-derived ASM at baseline, B) nonasthma-derived ASM at baseline vs. when treated with vitamin D, C) fatal asthma-derived ASM at baseline vs. when treated with vitamin D. The y-axis corresponds to the negative log (base 10) of P-values while the x-axis corresponds to the negative log (base 2) of the fold change for difference in expression between categories. Differentially expressed genes according to an , and IL8 by qRT-PCR using ACTB as a housekeeping reference gene shows similar differential expression results as those in Fig 2. (TIFF) S1 Table. Functional annotation clusters obtained with the NIH DAVID tool using all differentially expressed genes in fatal asthma-vs. non-asthma-derived ASM at baseline. Clusters with enrichment scores >1.5 are shown. Individual P-values listed correspond to EASE Scores, or modified Fisher Exact P-Values computed by DAVID. (DOCX) S2 Table. Individual annotation categories for top 838 fatal asthma-vs. non-asthmaderived ASM differentially expressed genes obtained from NIH DAVID cluster analysis tool. Categories selected were from clusters with enrichment scores >2.50 and with individual Benjamini-Hochberg corrected p-values <0.05 that correspond to known asthma-related structures and processes, plus two categories that met these criteria in the vitamin D treated cells vs. those at baseline. Genes listed were the differentially expressed ones for the corresponding category. (DOCX) S3 Table. Functional annotation clusters obtained with the NIH DAVID tool using all differentially expressed genes in non-asthma-derived ASM at baseline vs. when treated with vitamin D. Clusters with enrichment scores >1.5 are shown. Individual P-values listed correspond to EASE Scores, or modified Fisher Exact P-Values computed by DAVID. (DOCX) S4 Table. Functional annotation clusters obtained with the NIH DAVID tool using all differentially expressed genes in fatal asthma-derived ASM at baseline vs. when treated with vitamin D. Clusters with enrichment scores >1.5 are shown. Individual P-values listed correspond to EASE Scores, or modified Fisher Exact P-Values computed by DAVID. (DOCX) S5 Table. Individual annotation categories for top 867 non-asthma-derived ASM at baseline vs. with vitamin D treatment differentially expressed genes obtained from NIH DAVID cluster analysis tool. Categories selected were from clusters with enrichment scores >2.50 and with individual Benjamini-Hochberg corrected p-values <0.05 that correspond to known asthma-related structures and processes, plus categories that met these criteria in Table 3. Genes listed were the differentially expressed ones for the corresponding category. (DOCX) S6 Table. Individual annotation categories for top 711 fatal asthma-derived ASM at baseline vs. with vitamin D treatment differentially expressed genes obtained from NIH DAVID cluster analysis tool. Categories selected were from clusters with enrichment scores >2.50 and with individual Benjamini-Hochberg corrected p-values <0.05 that correspond to known asthma-related structures and processes, plus categories that met these criteria in Table 3. Genes listed were the differentially expressed ones for the corresponding category. (DOCX) S7 Table. Differential expression results for all comparison groups of asthma-related genes that were differentially expressed in ASM cells from a single male donor treated with vitamin D highlighted in previous work by Bosse, Y et al [17]. (DOCX) S8 Table. RNA-Seq Metrics Computed for Quality Control. The table contains the number of raw reads for the paired-end samples, the percentage of mapped reads among the total number of raw reads, the percentage of junction spanning reads among the mapped reads, the percentage of mapped bases that mapped to mRNA, the mean insert size of mapped reads, and the slope and R 2 values of ERCC spike-in dose response curves. Samples that were not included in differential expression analyses are highlighted in grey. The two samples marked by an asterisk ( Ã ) were dropped because of low percentage of junction spanning reads and mRNA bases, and mean insert size. The other three samples highlighted in grey were dropped because their pairs were not successfully sequenced. (DOCX)