Chemokine expression in the early response to injury in human airway epithelial cells

Basal airway epithelial cells (AEC) constitute stem/progenitor cells within the central airways and respond to mucosal injury in an ordered sequence of spreading, migration, proliferation, and differentiation to needed cell types. However, dynamic gene transcription in the early events after mucosal injury has not been studied in AEC. We examined gene expression using microarrays following mechanical injury (MI) in primary human AEC grown in submersion culture to generate basal cells and in the air-liquid interface to generate differentiated AEC (dAEC) that include goblet and ciliated cells. A select group of ~150 genes was in differential expression (DE) within 2–24 hr after MI, and enrichment analysis of these genes showed over-representation of functional categories related to inflammatory cytokines and chemokines. Network-based gene prioritization and network reconstruction using the PINTA heat kernel diffusion algorithm demonstrated highly connected networks that were richer in differentiated AEC compared to basal cells. Similar experiments done in basal AEC collected from asthmatic donor lungs demonstrated substantial changes in DE genes and functional categories related to inflammation compared to basal AEC from normal donors. In dAEC, similar but more modest differences were observed. We demonstrate that the AEC transcription signature after MI identifies genes and pathways that are important to the initiation and perpetuation of airway mucosal inflammation. Gene expression occurs quickly after injury and is more profound in differentiated AEC, and is altered in AEC from asthmatic airways. Our data suggest that the early response to injury is substantially different in asthmatic airways, particularly in basal airway epithelial cells.


Introduction
The epithelium lining human lung airways serves as a primary defense interface against particulates, inhaled pollutants, pathogens, and xenobiotics. Basal airway epithelial cells (AEC) comprise approximately 30 percent of the epithelium in the major airways. These cells operate as the progenitors for ciliated and goblet columnar cells within the airway [1][2][3] Human Research Protection (45 CFR46.102 (f)) did not apply, as the donors were deceased and that we did not obtain information about or samples from living individuals, so that the research described in this manuscript did not constitute engagement in human subjects research. No donated tissue and organs came from any vulnerable populations. Lungs were collected from two groups of donors: those without known lung disease, and those with a history of chronic asthma. For the latter group, no specific information generally was available concerning disease severity, duration, or treatment. No other clinical or identifying information provided by the Gift of Hope was used or is disclosed in this study.
To obtain either basal cell or differentiated cell phenotypes, AEC were grown in submersion or ALI culture (21 days) respectively [17,22,23]. Passage 1 cells used in ALI culture were seeded at 1.0 x 10 5 cells onto collagen IV-coated, 12.5 mm diameter Millicell CM culture inserts (pore size 0.4 μm) (Millipore, Inc.). Cells were grown to 80% confluence, the apical medium was removed and cells were fed via the basolateral side every other day to day 21. Passage 1 cells used in submersion culture were seeded similarly onto collagen IV-coated plastic wells, 12.5 mm diameter, and grown to 100% confluence. Cell morphology was confirmed in cultures from each donor using immunostaining antibodies directed against cytokeratin 5 (CK5) to mark basal cells (either goat IgG sc-17090 or mouse IgG1 sc-32721, both used at 1:250 dilution Santa Cruz) [24,25], alpha-acetylated tubulin (AA-tubulin) to mark ciliated cells (mouse IgG2b #32-2700, 1:500 dilution, Zymed) [26,27] or MUC5AC (either goat IgG sc-16910 or rabbit IgG sc20118, both used at 1:250 dilution, Santa Cruz), to mark mucoid cells [28,29]. Cells were first fixed in 100% ethanol at 4˚C for 30 minutes, after which cells were washed, and primary antibodies were incubated overnight at 4˚C. Secondary antibodies directed against each with non-overlapping conjugated fluoroprobes then were added (2 mg/ ml each as required) at 20˚C for 1 hr: donkey-anti rabbit IgG conjugated with Alexa Fluor 594; donkey-anti mouse IgG (H+L) conjugated with Alexa Fluor 488; donkey anti-goat IgG (H+L) conjugated with Alexa Fluor 594; donkey anti-goat IgG (H+L) conjugated with Alexa Fluor 488; and donkey anti-mouse IgG (H+L) conjugated with Alexa Fluor 594 (all Invitrogen). Nuclei were labeled with Hoechst 33258. Cells were imaged using an Olympus DSU disk scanning confocal system on an Olympus IX81 microscope platform equipped with a Photometrics Evolve-back-thinned CCD camera and UPlan Apo 40x or 60x objectives. Images were assembled from z-stacks using ImageJ software (Wayne Rasband, National Institutes of Health, Bethesda, MD).

Mechanical injury experiments
Cells grown to confluence in either culture system were injured using a mechanical rake to ensure a large number of wound edges [30], then incubated for up to 24 hr. Control cells without mechanical injury (MI) were collected at the same time intervals.

Microarray processing
Total RNA was extracted from cells using QIAzol lysis reagent followed by RNeasy mini kits (Qiagen, Inc., Valencia, CA). Gene expression profiling was assessed using HumanHT-12 v4 expression BeadChips (Illumina, San Diego, CA) that included probes for more than 47,000 genome-wide transcripts.

Differential expression analysis
We first grouped the expression data into 12 groups, each group with only one variant of injury versus no-injury. The 12 groups contain each time point (2,8, and 24 hour), two culture conditions (submersion and ALI) and two disease states (asthma and normal). Two major methods of analysis, weighted gene correlation network analysis (WGCNA) and differential expression based analysis, were applied to these groups. For microarrays passing initial quality control, expression levels were extracted using Limma after background subtraction and quantile normalization using the method of Xie et al. [31]. Significant DE was determined following correction for multiple testing and false discovery by the method of Benjamini and Hochberg (23) using Limma [32]. Genes in DE between ALI and submersion conditions (absolute FC ! 1.5, corrected P value < 0.05) for every time point under consideration were extracted.
Gene expression pattern analysis was performed to identify the sets of genes differentially expressed under each culture condition over time [33]. Heat-maps and a linear plot for the significantly differentially expressed genes were generated to indicate the levels of expression and the relationship among the genes at different times and conditions. Enrichment analysis for identification of GO categories, diseases, phenotypes, pathways, and other features over-represented in the sets of differentially expressed genes, was performed using online analytical tools available in the Lynx integrated system [34]. Statistical significance was estimated using Bayes factor and P-value calculated for each category under consideration. Bonferroni correction was applied to the results of analyses to correct for multiple testing. MIAME-compliant raw data has been deposited in the Gene Expression Omnibus (GEO) site (http://www.ncbi.nlm. nih.gov/geo) (accession numbers GSE 59128 and GSE109170).

Network-based gene prioritization
Genes differentially-expressed under each condition were submitted to the Lynx/PINTA server. The PINTA [35] heat kernel diffusion algorithm was used to generate a gene association score that estimates the strength of functional associations between the differentially expressed genes and prediction of additional genes potentially involved in the inflammation/wound healing response under each condition. The original heat kernel diffusion algorithm from the PINTA tool was modified to accommodate weighted data types. The STRING platform was used as the underlying genome wide functional association network for gene prioritization and network reconstruction [36,37]. Two iterations for network diffusion; 0.5 as the diffusion value, and 5,000 network randomizations for estimating p-values, were used as parameters for the heat kernel diffusion algorithm. Bonferroni correction was applied to the results of analyses to correct for multiple testing. Genes with the score ! 0.03 were extracted for each condition. The resulted patterns and network-based predictions were visualized by the STRING webserver using default settings for confidence, evidence and molecular action [36].

Weighted gene co-expression network analysis
Systems level analysis of cellular responses was done by WGCNA [38,39]. This analysis was performed separately for each group as well for the whole experiment to avoid the dilution of the co-expression patterns by mixing the experimental conditions. The initial step included the reconstruction of the co-expression network based on Pearson correlation. The threshold for this network was determined by the connectivity and degree distribution. Module detection was performed on the resulting co-expression network by hierarchical clustering. Each module contained subsets of genes in the input expression files; modules did not overlap with each other. The genes within each module were considered to be co-expressed. The modules then were evaluated by whether the overall expression of genes within the module correlated, positively or negatively, with injury/no-injury as the condition variable. An estimated P-value was generated for each module to estimate the significance of the corresponding correlation. Functional enrichment analysis then extracted significant categories from each module which was compared across the groups.

Epithelial transcriptomes of normal un-injured basal AEC and dAEC are different under standard quiescent culture
Basal AEC grown in submersion culture expressed cytokeratin 5 (KRT5) and did not express significant quantities of either gel-forming mucin 5AC (MUC5AC), characteristic of goblet epithelial cells, or of alpha-acetylated tubulin, characteristic of ciliated epithelial cells, as reported by the previous studies [6,22]. As expected, dAEC culture demonstrated abundant expression of both MUC5AC and alpha-acetylated tubulin in addition to cytokeratin 5, indicating the presence of goblet and ciliated cells, respectively ( Fig 1A). The whole genome transcriptomes of normal dAEC and basal AEC from four donor lungs also demonstrated significant differences as assessed by both PCA and volcano plots of DE ( Fig 1B) as expected [6] for differences in cell type and development. Table A in S1 File provides a list of genes in DE between the two culture conditions. The enrichment analysis of the genes differentially expressed (logFC>1.3, FDR-corrected P-value < 0.05) between the normal uninjured basal AEC and differentiated dAEC was performed using both the Lynx [34,40] and ToppGene [41] bioinformatics platforms. There were substantial differences between these two culture conditions at 8 hr after MI. The genes highly expressed in the basal AEC cells were predominantly involved in keratinization (GO:0031424) and cell cycle (GO:0007049). However, the genes highly expressed in the dAEC cells at the same time point were enriched in the categories related to cytokine activity (GO:0005125), extracellular space (GO:0005615), and interferon signaling (GO:0060333) ( Table B in S1 File).

Epithelial transcriptomes in normal basal and dAEC differ after mechanical injury
Mechanical injury to normal AEC and dAEC was done using a rake that generated a substantial number of cells at wound edges [30]. The analysis of differential expression of genes in both normal differentiated and basal AEC at 2, 8 and 24 hours after MI demonstrated (Fig 2A) a significant number of DE genes. Table C in S1 File lists 87 genes of interest to inflammation; as shown in Fig 2B, a number of inflammatory cytokine and chemokine ligands, and key signaling intermediates for these ligands, were up-regulated in both differentiated and basal AEC 2 to 24 hr hours after mechanical injury.
Tables D-G in S1 File list genes related to cell cycle, differentiation, proliferation, and cell migration. Using one-tailed paired t-tests, we examined gene expression after injury in differentiated AEC as well as basal AEC for genes associated with cell cycle and differentiation. In both culture conditions at 8 hr after MI, there were genes in significant DE for both cell cycle (P < 0.01) and differentiation (P < 0.02) related genes. Therefore, genes involved in these two biological process have significant (P < 0.05) difference between dAEC and basal AEC as expected.
Enrichment analysis was performed for identification of the functional categories overrepresented in the sets of the genes differentially expressed in normal vs injured AEC or dAEC cells. Table 1 presents the results of the enrichment analysis of the DE genes in normal dAEC and AEC 8 hr after injury. Both basal AEC and dAEC demonstrated significant activation of molecular pathways involved in cytokine activity (GO:0001817, REACT_75790), inflammation (GO:0006954), and to defense responses against bacteria (GO:0002684). Activation of the response to injury was observed in both cell culture conditions but different sets of genes were activated in the two different cell types (Table 1).
Systems level analysis of these responses was characterized using weighted gene co-expression network analysis (WGCNA) [38,39]. The analysis was done for each time point and cell type using only significantly differentially expressed genes (FC > 1.5, adjusted P value < 0.05) to identify gene functional modules specific for the condition under observation. An absolute correlation coefficient > 0.8 and adjusted P value < 0.05 was required for significance. Table H in S1 File presents the statistically significant modules characteristic for each condition. The significant modules detected at 8 hr after MI were enriched for cytokine activity (GO:0005125), interleukin-1 receptor binding (GO:0005149), and inflammatory responses (GO:0006954) in both differentiated and basal AEC, while differentiated AEC also showed enrichment for toll-like receptor pathways (GO:0034142, GO:0002224, GO:0034138) ( Table I in S1 File). A greater number of GO functional categories were identified as significant for dAEC as compared to basal AEC.
Reconstruction of the molecular networks using DE genes (P < 0.05 and FC > 1.5) as seed genes for dAEC and AEC at each time point are presented in Fig 3. There was an increase in membership and complexity of the anti-inflammatory response at 8 and 24 hr, as compared to 2 hr, after MI for both basal and differentiated AEC. However, at 8 hr after MI, the DE genes network for differentiated AEC (86 nodes) are more complex than their basal AEC counterparts (47 edges) (Fig 3 and Table J in S1 File). Especially at 8 hr after MI, the differentiated AEC network contains twice as many edges as basal AEC network and the differentiated AEC network has a higher average degree (8.047) than basal AEC network (6.809). These data indicated that in the immediate hours after injury, both differentiated and basal AEC elicited gene expression that is related to a significant inflammatory and pro-defense response that in differentiated AEC is more intense.

Epithelial transcriptomes in basal and dAEC differ after mechanical injury specific to asthmatic cells
As with cells collected from normal donor lungs, transcriptomes of dAEC and basal AEC collected from asthmatic donor lungs demonstrated appreciable numbers of ciliated and goblet cells when grown in air-liquid interface and a significant number of up-regulated and downregulated genes illustrated by volcano plot (FC ! 1.5 and corrected P value [FDR] < 0.05) (Table A in S1 File).
After the mechanical injury, both differentiated and basal AEC from asthmatic subjects expressed a limited, discrete set of up-regulated genes, and fewer down-regulated genes ( Fig  4A). Furthermore, the expression of cytokine and chemokine ligands was similar in differentiated and basal AEC (Fig 4B). These patterns were consistent to the ones observed in normal AEC and dAEC.
Reconstruction of the molecular networks for the significantly DE genes (p-value<0.05, fold change > 1.5) using the STRING10 server is presented in Fig 5. The connections in dAEC were more complex and denser at 8 hr (with 771 edges, 9.071 average degree) versus 2 hr (with 76 edges 4.75 average degree). There were decidedly few connections observed at 24 hr with only 6 edges, which fit with the clear change in cytokine/chemokine expression changes and the far fewer functional pathways identified at that time point.  Epithelial cell chemokine expression after injury WGCNA analysis for each condition only detected 1 or 2 significant modules (Table K in S1 File). In 8 hr, the significant module detected for dAEC culture shows enrichment for cytokine activity, inflammatory response, interferon signaling pathway, cell death, TNF pathway and toll-like receptors while in AEC culture the detected module is enriched for regulation of phosphorylation, cell death, and differentiation (Table L in S1 File).

Epithelial transcriptomes in normal and asthmatic basal and dAEC differ after mechanical injury
The analysis of the transcription in both differentiated dAEC and basal AEC after mechanical injury demonstrated significant differences for the cells collected from normal and asthmatic subjects at all time points after injury. These included variations in expression patterns of  ADAM8 CCL2 CD83 CSF2 CYP1B1 HERC5 IFIH1 IKBKE IL1A IL1B IL23A IL36G IL6 ISL1  KLF4 MAP2K3 NFKB2 PDE4B PELI1 PTGS2 STAT5A TICAM1 TNF TNFAIP3 ZBP1  ZC3H12A   Reconstruction of molecular networks for normal basal and differentiated AEC after mechanical injury. String networks generated from WGCNA modules using the genes in significant DE at 2-24 hr after injury for each cell type to demonstrate functional relationships in molecular networks. Highly connected network activation was richer and more substantial in differentiated (dAEC) versus basal AEC, particularly at 2 and 8 hr after MI, but less so at 24 hr. Probe labels are as shown.
cytokines, chemokines and genes involved in related signaling pathways both in asthmatic basal and differentiated AEC. Fig 6 shows the gene DE between asthmatic and normal AEC, and dAEC at 8 hours after the injury. The expression of the genes involved in cytokine response to injury generally was lower in asthmatic epithelium both in AEC and dAEC. Epithelial cell chemokine expression after injury

Fig 5. Reconstruction of molecular networks for asthmatic basal and differentiated AEC after mechanical injury.
String networks generated from WGCNA modules using the genes in significant DE at 2-24 hr after injury for each cell type to demonstrate functional relationships in molecular networks. Highly connected network activation was similar at 2 hr, richer and more substantial in differentiated (dAEC) versus basal AEC at 8 hr, and less significant at 24 hr after MI in dAEC compared to basal AEC. Probe labels are as shown.
However, the expression of IL8, IL1A, IL1RL1, as well as CCL20 genes was higher in asthmatic basal AEC in comparison to normal basal AEC. The expression of IL-19, TNFSF14, as well as the NFKB-independent IL1RL1 interleukin-like receptor also was increased in asthmatic dAEC. These cells also demonstrated higher expression levels of the genes that encode the interferon-induced IFIT proteins (IFIT1, IFIT2, and IFIT3). The DE of key genes at each time point after MI are presented in Table M in S1 File. These data suggest a substantially different pattern of gene expression in asthmatic epithelial cells as compared to normal AEC.

Enrichment analysis
The enrichment analysis of genes differentially expressed (p < 0.05, FC >1.5) at each time point has revealed the significant activation of functional categories related to transcription  Table 2, basal AEC; and Table 3, differentiated AEC). Asthmatic AEC transcriptomes were enriched for the genes involved in the pathways related to transcription and cell death 2 to 8 hr after MI, and to cell cycle 24 hr after MI in comparison to functional pathways generated in normal AEC (Table 1). Asthmatic DE genes also contained fewer cytokine and chemokine genes. Examination of functional pathways based on the number of DE genes, confirmed the initial analysis of gene DE enrichment analysis. At each time point in many of the functional pathways, AEC from normal donors had a greater number of DE genes and higher Bayes factors compared to AEC from asthmatic donors, and differentiated cells had a greater number of DE genes compared to basal AEC (Fig 7).

Discussion
An important issue in tissue repair after airway mucosal injury is to define the mechanisms by which inflammatory factors, growth factors, and other effectors influence the renewal and differentiation of resident cells. Especially important are early events after injury that initiate or perpetuate inflammatory responses in the immediate and neighboring regions. Because multiple factors are produced in the first hours after injury, it is critical to understand the initial responses of cells so as to determine whether these signals, alone and in combination, serve as pro-or anti-reparative signals or as autocrine or paracrine signals for repair and inflammation in neighboring cells and tissue. While many studies to date have evaluated the role of an individual protein or a handful of factors in elegant 'reductionist' models to illustrate their key roles in these processes, few have examined the orchestrated response of epithelial cells, particularly in airways, after discrete injury. These orchestrated responses, from transcriptional to post-translational, serve to coordinate reparative and inflammatory signals that are appropriate in normal responses to injury, and that may become aberrant in disease states depending on context.  -08  51  1952  MYC,SOX9,GADD45B,SPRY2,ABL2,TP53INP1,TNFSF10,GAS1,ETS1,IRAK2,  DAPK1,PPP1R15A,ITGA5,PHLDA2,ADAM8,F3,TNFRSF12A,SEMA3A,BMP2,  PPARGC1A,MALT1,TSC22D1,BTG2,SPHK1,TAOK1,VAV3,PMAIP1,CYP26B1,  RARB,PHLDA1,TSC22D3,WNT7A,WNT7B,G0S2,DUSP6,ID1,ID3,EDN1,LGALS3,  FOSL1,TGFA,PPP2CB,LATS2,CYR61,ANGPTL4,SERPINB2,PAK2,ATF3,TNFAIP3,  CKAP2,IL1A GO:0042981 regulation of apoptotic process 5.4E-11 3. 9E-08  44  1519  MYC,SOX9,GADD45B,SPRY2,ABL2,TP53INP1,TNFSF10,GAS1,ETS1,IRAK2,  DAPK1,PPP1R15A,ITGA5,ADAM8,F3,TNFRSF12A,BMP2,PPARGC1A,MALT1,  TSC22D1,BTG2,SPHK1,VAV3,PMAIP1,RARB,PHLDA1,TSC22D3,WNT7A,G0S2 We have examined the AEC transcriptome in the first hours after injury with an emphasis on the expression of cytokines and chemokines that can influence inflammation in the airway mucosa and submucosa, particularly that with identified importance in asthma. We identified a distinct gene expression profile of such mediators in both differentiated and basal normal AEC that left unchecked would contribute to inflammatory cell influx into the mucosa. Differentiated AEC demonstrated a richer gene expression pathways of pro-inflammatory chemokines than basal AEC, suggesting that a normal, homeostatic airway epithelium is more capable of responding to injury with pro-inflammatory mediator expression. Basal AEC have a central role in repair after injury in airways, from the initial steps of spreading and migration to later proliferation and phenotype shifting into needed new cell subtypes [1,2,4,42,43]. In contrast to the response seen in differentiated AEC, these pluripotent basal cells were less capable of expressing key cytokines and chemokines. Gene expression networks relevant to the cytokines response differed significantly in the first 24 hr after injury. In an airway with chronic damage to or loss of ciliated and goblet cells, as is frequently encountered in severe asthma, basal AEC are more dominant, and thus the pro-and anti-inflammatory responses signaled by them may differ from that seen with an intact, differentiated epithelial layer. These observations have clear implications for the responses of the neighboring structural and inflammatory cells and may become magnified over time with repeated cycles of injury and (aberrant) repair.

Go term Name P value
Adjusted Pvalue
Of interest is the observation that more differentiated cells have a richer DE early (2 and 8 hr) after the injury reflected by the variety of functional categories in enrichment analysis, than basal cells. Dissecting this difference is hampered due to the fact that the differentiated ciliated and goblet columnar cells require a mixed culture in which basal cells provide for anchoring and survival and cannot be isolated and grown in pure culture. The relative contribution of each cell type then is not clear from our findings, but may become important in models of chronic airway inflammation in which goblet cell metaplasia is prominent [60][61][62][63].

Conclusions
Asthmatic airway inflammation frequently is initiated by injury to the epithelial cells lining central airways. Gene expression is essential for orchestrating the early response to epithelial cell injury. We report that basal and differentiated airway epithelial cells collected from asthmatic lungs have diminished immediate cytokine and chemokine response to injury in comparison to non-asthmatic, normal AEC and dAEC. Gene expression in asthmatic cells instead is centered on reparative processes including proliferation, cell cycle, and differentiation. While airway inflammation is the sine qua non of asthma, intrinsic dysfunction of the asthmatic airway epithelium as demonstrated by decreased expression of inflammatory pathways after injury suggests that paracrine factors are required to augment epithelial-derived inflammatory responses after injury.
Supporting information S1 Letter. Letter from the University of Chicago institutional review board. (PDF) S1 File. Single file that contains supporting information Tables A-M