RNA-Seq and Network Analysis Revealed Interacting Pathways in TGF-β-Treated Lung Cancer Cell Lines.

Whole transcriptome shotgun sequencing (RNA-Seq) is a useful tool for analyzing the transcriptome of a biological sample. With appropriate statistical and bioinformatic processing, this platform is capable of identifying significant differences in gene expression within the transcriptome and permits pathway and network analyses to determine how these genes interact biologically. In this study, we examined gene expression in two lung adenocarcinoma cell lines (H358 and A459) that were treated with transforming growth factor-β (TGF-β) as a model for induction of the epithelial-to-mesenchymal transition (EMT), commonly associated with disease progression. We performed this study in order to illustrate a workflow for identifying interesting genes and processes that are regulated early in EMT and to determine their gene pathway/network relationships and regulation. With this, we identified 137 upregulated and 32 downregulated genes common to both cell lines after TGF-β treatment that represent components of multiple canonical pathways and biological networks associated with the induction of EMT. These findings were also verified against reposited Affymetrix U133a expression profiles from multiple trials examining metastatic progression in patient cohorts (n = 731 total) to further establish the clinical relevance and translational significance of the model system. Together, these findings help validate the relevance of the TGF-β model for the study of EMT and provide new insights into early events in EMT.

Introduction

In the study of cancer mechanisms, the use of next-generation sequencing (NGS) to gain valuable biological insights via pathway/network analyses is still limited, because of the relative infancy of these high-throughput, computationally intensive methods. Novel bioinformatics methods offer a means to reveal previously unappreciated mechanisms that have potential of producing new candidate biomarkers or targets for pharmaceutical intervention and, thereby, impact clinical practice once fully developed.

Whole transcriptome shotgun sequencing (RNA-seq) is a biochemical sequencing method that interrogates the entire transcriptome of a sample at a particular time using NGS.1 In contrast to studies performed using gene expression microarrays, which target specific sequences for measurement with probes, RNA-seq takes into account greater base coverage of DNA sequences as well as the different types of RNA (microRNA, transfer RNA, etc.). RNA-Seq accomplishes this by amplifying and sequencing the RNA content of a sample using specific, bar-coded RNA sequences of interest, and aligning the sequences by a reference genome to determine the gene biomarkers associated with those bar-coded sequences. This high-throughput technique uses extensive bioinformatics computation to perform the necessary sequence alignments. The resulting data can be used to analyze gene expression, single nucleotide polymorphism, alternative splicing, and other genomic studies of interest to cancer biology.

Pathway and network analysis via modeling biological systems as computable networks is another bioinformatics procedure that converts existing biological information on canonical pathways and biological networks into searchable databases, then uses computing algorithms to find linkages between related pathways and networks.1,2 This permits researchers to filter the existing literature for useful, non-obvious relationships between biological systems that can be used to build a biological model for a system that can be further investigated with experimentation.

Lung cancer is the leading cause of cancer-related deaths in the United States and is projected to account for an estimated 159,260 deaths in 2014.3 Dissemination of primary tumor cells to distant structures is the primary event that leads to mortality in lung cancer and is a major focus of disease staging for prognosis and treatment plan development. Biological processes commonly associated with metastasis include: increased tumor cell motility, altered cellular adhesion, enhanced capacity for remodeling extracellular matrices for extravasation and intravasation, resistance to apoptosis (cell death), induction of angiogenesis, and the ability to thrive in the secondary site.4–8 Many of these phenotypic features are consistent with the popular concept that metastatic progression in epithelial cancers is driven by a type III epithelial-to-mesenchymal transition (EMT).9,10 Furthermore, this apparent phenotypic “transdifferentiation” is also associated with the adoption of specific “cancer stem cell” characteristics, including capacity for self-renewal, altered surface antigens, modulation of signaling pathways, and transcriptional regulation mechanisms.11–14 Perhaps the most archetypal example of EMT features is the loss of E-cadherin expression and upregulation of alternate adhesion molecules, such as N-cadherin or fibronectin.11,15 Several cell-based models for EMT have emerged that are thought to be generally consistent with the biology observed in early metastatic events clinically, with the model using transforming growth factor-beta (TGF-β) induction among the most well studied.16–19 In this study, we leverage the strengths of RNA-Seq and pathway/network analysis to further refine our understanding of events immediately resulting from TGF-β induction as a model of early EMT events. Our goal is to identify features that further validate this model, identify potentially underappreciated processes and/or effector molecules intimate to EMT for future study, and reveal novel candidate biomarkers that may have translational significance.

Materials and Methods

Sample collection for the TGF-β model of EMT induction and phenotypic characterization

H358 and A549 lung adenocarcinoma cell lines were purchased from American Type Culture Collection (ATCC) and maintained in RPMI 1640 supplemented with 2.5% fetal bovine serum (FBS) at 37 °C in a humidified 5% CO2 atmosphere. The cell lines were cultured for 3 days in the presence or absence of 10 ng/mL TGF-β in RPMI 1640 supplemented with 2.5% FBS. After treatments, cells were rinsed in cold PBS, scraped, and total RNA isolated using an RNeasy kit (Qiagen GmbH) according to manufacturer-suggested protocols. A total of four conditions (H358 ± TGF-β; A549 ± TGF-β) in duplicate provided eight samples for analysis.

A parallel set of samples were cultured and treated as defined above for western blot analysis of apoptosis-related targets. Briefly, cellular lysates were prepared for western blot analysis as previously described,20,21 with proteins (30 μg/lane) resolved on 10–20% Criterion TGX tris-glycine gels (Bio-Rad) and transferred to nitrocellulose overnight. Blots were probed with the Apoptosis Antibody Sampler Kit (Cell Signaling Technologies) for caspase-3, cleaved caspase-3, PARP, cleaved PARP, as recommended, and normalized to β-actin. Also, cellular proliferation was evaluated in both the A549 and H358 cell lines using the MTT Cell Proliferation Assay kit (ATCC) using recommended protocols, with cell treatments consistent with those defined above (adapted for the 96-well format). Each condition was evaluated a minimum of two times with triplicate sampling of cellular numbers.

Whole transcriptome shotgun sequencing

For each of the eight samples, total RNA (1 μg/sample) was processed into an NGS-compatible library with the TruSeq RNA sample preparation kit (Illumina). Poly-A containing mRNA was purified by bead selection and the collected mRNA sized by heat fragmentation. After conversion to cDNA, bar-coded adaptors were added by A/T supported ligation, and the resulting library was amplified by PCR. Equal amounts of each library were pooled and submitted to paired-end sequencing on a HiSeq2000 sequencer (Illumina). Clusters were sequenced for 50 bases from each end. The Tophat (version 2.0.0) software package was used to align the raw RNA-seq data to the human genome and identify individual read sequences via barcodes followed by the use of Cufflinks (version 2.0.0) to count and convert the read results into genes with annotation via Ensembl.

Statistical analysis

After the gene expression levels were derived, expression values were averaged across four lanes, log2 transformed, and subjected to statistical analyses using Agilent Genespring GX 12.6. Within each cell line, the gene lists were processed using an independent Student’s t-test and fold-change analysis to identify upregulated and downregulated genes in TGF-β-treated specimens, relative to control. A significance threshold of P ≤ 0.05 and absolute fold change ≥2.0 were used to define modulated genes of interest. Venn diagrams were used to illustrate commonality in genes being up/downregulated in both cell lines.

Pathway and network analysis

A pathway analysis was performed to identify associations of the genes identified in our cell line studies with canonical pathways and gene networks. The pathways and networks were generated through the use of Ingenuity Pathway Analyses (IPA; Ingenuity® Systems, www.ingenuity.com). Thresholds of two-fold or greater in changes in expression and a P-value of 0.05 or less for significance were used to filter the findings from the analysis with IPA software.

Correlation of gene expression to metastasis in clinical lung cancer samples

We identified several studies in the literature that used Affymetrix U133a expression arrays to assess global transcription levels of genes in primary lung cancer tumor tissue associated with metastatic progression to the locoregional lymph nodes. In total, 731 raw CEL files each representing separate and individual tumors were combined and preprocessed using the RMA preprocessing in the SimpleAffy Bioconductor package in the statistical analysis package R.22 Once expression estimates were generated, individual genes identified as differentially expressed by TGF-β treatment in cell lines in the above studies were correlated to surrogates of metastatic progression in the clinical specimens (ie, using disease stage and/or disease-free survival annotation as available) using a Cox regression function in the R survival package. Briefly, profiles originated from studies conducted by Hoang et al.23 (n = 22), Landi et al.24 (n = 107), Raponi et al.25 (n = 130), Shedden et al.26 (n = 338), Shah et al.27 (n = 30), and Spira et al.28,29 (n = 104).

Results

Integrity of the adenocarcinoma cell lines with TGF-β induction

Prior to expression profiling, it was crucial to demonstrate that culturing conditions and/or treatment with TGF-β does not result in widespread apoptosis or unexpected changes in cellular proliferation. With this, we evaluated via western blot analysis specific indicators of apoptosis induction, including caspase-3 and PARP activation (cleavage), to find no significant differences in any of the conditions tested (Supplementary Fig. 1, panels A and B). We did identify a change in cellular proliferation rate upon chronic treatment with TGF-β that resulted in a 25.2% and 44.1% lower cellular density at day 3 for the H358 and A549 cells (Supplementary Fig. 2, panels A and B); however, these observations are not unexpected and have been reported by others.30,31 Together with the adoption of more fibroblastic morphological features upon TGF-β treatment and lack of vacuoles (data not shown), we determined the general condition of the cells to be suitable for evaluation via RNA-sequencing.

Comparison of RNA-seq findings of the TGF-β-induced adenocarcinoma cell lines

Paired-end whole transcriptome sequencing for two lung adenocarcinoma cell lines, H358 and A549, was accomplished as a means to investigate modulations in gene expression and pathway/network recruitments resulting from a TGF-β model of EMT. Our objective was to identify features that further validate this model, identify potentially underappreciated processes and/or effector molecules intimate to EMT, and reveal novel candidate biomarkers that may have translational significance.

After expression values were calculated from the RNA-seq results, a total of 54,607 candidates with unique Ensembl gene IDs were identified. When subjected to filtering with statistical (P ≤ 0.05) and fold change (≥2) thresholds, the H358 cell line was found to have a total of 1177 upregulated genes and 1108 downregulated genes, whereas the A549 cell line had 525 upregulated and 167 downregulated genes, all relative to the cell line controls (eg, no treatment). These general findings are illustrated in Figure 1 as “box and whisker” plots. In general, there were fewer downregulated genes that met our statistical and fold-change thresholds than those upregulated. In comparing the findings of unique genes modulated in each cell line, a total of 137 upregulated and 32 downregulated genes were shared between H358 and A549 cell lines with the TGF-β induction, as shown in Figure 2 as a Venn diagram. Table 1 shows the top 10 common upregulated and downregulated genes in the two cell lines after TGF-β treatment. The largest increases in gene expression in this model were SERPINE1 in the H358 cells, with a 344.39-fold change, and DOCK2 in the A549 cells, with a 103.29-fold change (Table 1). Conversely, we observed the most dramatic decreases in gene expression in this model to be MUC5B in the H358 cells, with a 64.65-fold change, and ST6GALNAC1 in the A549 cells, with a 118.58-fold change (Table 1). Supplementary Tables 1–4 provide additional examples of genes up- and downregulated in response to the TGF-β treatment, separated by cell line and direction of change. Many of these targets are typically associated with or considered hallmark biomarkers of the EMT thought to mechanistically underlie tumor metastasis in vivo.9,32

Figure 1

Comparison of expression value distributions between samples. Average expression value was calculated in each cell line with or without TGF-β treatment (control), respectively.

Figure 2

Intersection of upregulated and downregulated genes across cell lines.

Table 1

Top 10 upregulated and downregulated overlapping genes.

GENES	P-VALUE-358	FOLD CHANGE-358	P-VALUE-549	FOLD CHANGE-549
Upregulated
CLRN3	<1.0E-005	4.00	<1.0E-005	9.00
ISM2	2.51E-005	29.50	4.62E-003	22.25
LPAR5	2.98E-004	12.43	1.52E-002	33.54
IQSEC3	6.19E-004	5.74	<1.0E-005	4.00
SERPINE1	6.20E-004	344.39	6.93E-003	87.83
MAOB	7.06E-004	6.30	1.99E-002	5.05
DOCK2	7.86E-004	59.46	8.69E-003	103.29
CLDN14	9.94E-004	23.36	3.64E-002	15.02
JUN	1.02E-003	6.98	2.00E-002	15.46
CGB	1.22E-003	244.82	4.19E-002	19.37
Downregulated
NRAP	<1.0E-005	2.00	3.97E-002	42.67
ABCC12	9.88E-004	3.57	4.86E-002	8.06
INHBB	1.70E-003	62.98	4.00E-002	10.02
SSTR5	4.46E-003	23.47	4.29E-003	56.03
APOD	5.42E-003	3.29	3.83E-002	9.40
MUC5B	6.46E-003	64.65	3.10E-002	48.93
FGFBP1	8.58E-003	4.91	2.86E-002	22.37
ST6GALNAC1	1.27E-002	54.91	3.59E-003	118.58
TSPAN8	1.58E-002	6.32	4.14E-002	18.44
UGT2B15	1.66E-002	6.86	2.86E-002	33.85

Note: P-value and fold change were obtained by comparing the gene expression level in each cell line with and without TGF-treatment.

Biological pathways modulation in the TGF-β-induced adenocarcinoma cell lines

We further classified these genes with varying expression into different pathways according to their canonical pathways using the IPA suite. In order to distinguish the relative impact of TGF-β on each cell line, we provide a selection of the pathways found to be modulated separated based on being observed in either the A549 cells (Table 2) or the H358 cells (Table 3). The main point of dissimilarity between the cell lines is that the A549 cells appear to be very focused on the activation of pathways associated with hepatic fibrosis/hepatic stellate cell activation, human embryonic stem cell pluripotency, and granulocyte adhesion and diapedesis, which are consistent with active transdifferentiation events that are part of the EMT program.9,32 Although, it should be noted that the activation of pathways such as axonal guidance signaling and leukocyte extravasation signaling is also observed in the A549 cells, but these pathways are observed to be much more dominate in the H358 cells, suggesting that both cell lines are going through similar changes but potentially to different degrees. The Supplementary results further illustrate these commonalities by listing the top 10 upregulated and top 10 downregulated pathways in the two cell lines after TGF-β induction (Supplementary Table 5). The most significantly up- and downregulated pathways common between the cell lines were Gα12/13 signaling and serotonin degradation, respectively. We also observed that several canonical pathways with involvement of the 137 upregulated genes shared between the two cell lines were associated with cancer-related elements, such as p53 signaling, HER-2 signaling, and bladder cancer signaling (Supplementary Table 5). The downregulated genes did not seem to be associated with cancer pathways in particular, but this is expected because of the relatively low number of downregulated genes found.

Table 2

A selection of modulated ‘canonical pathways’ from the A549 lung adenocarcinoma cell line.

PATHWAYS	P-VALUE (−LOG)	R ATIO	GENES INVOLVED
Upregulated
Hepatic Fibrosis/Hepatic Stellate Cell Activation	6.76E00	7.61E-02	COL4A1,IGFBP5,MMP2,COL4A2,PDGFB,MYL9,COL5A1,COL1A1,IGF2,IGF1,COL6A3,COL22A1,SERPINE1,MMP1,TNFRSF11B
Human Embryonic Stem Cell Pluripotency	3.15E00	5.97E-02	NTF4,S1PR5,WNT9A,BMP2,LEF1,LEFTY2,PDGFB,INHBA
Granulocyte Adhesion and Diapedesis	2.97E00	5.08E-02	CLDN4,CXCL12,CCL14,MMP2,CLDN14,CLDN9,MMP1,HSPB1,TNFRSF11B
Axonal Guidance Signaling	2.4E00	3.24E-02	EFNA2,GNG4,ADAMTS8,NTF4,WNT9A,BMP2,CXCL12,MMP2,PDGFB,MYL9,IGF1,ABLIM3,ADAM19,SEMA7A
Role of Osteoblasts, Osteoclasts and Chondrocytes in Rheumatoid Arthritis	2.35E00	4.11E-02	COL1A1,JUN,IGF1,WNT9A,BMP2,LEF1,MMP1,IL11,TNFRSF11B
Gα12/13 Signaling	2.16E00	5.13E-02	MYL9,JUN,CDH4,TBXA2R,LPAR5,CDH19
Chondroitin Sulfate Biosynthesis	2.13E00	7.41E-02	CHST1,XYLT1,HS3ST6,SULT1B1
Leukocyte Extravasation Signaling	2.1E00	4.04E-02	TIMP3,CLDN4,RASGRP1,CXCL12,MMP2,CLDN14,CLDN9,MMP1
TGF-β Signaling	2.08E00	5.75E-02	JUN,BMP2,SERPINE1,INHBA,PMEPA1
Heparan Sulfate Biosynthesis	2.04E00	7.02E-02	CHST1,XYLT1,HS3ST6,SULT1B1
Downregulated
Guanosine Nucleotides Degradation III	2.71E00	1.54E-01	GDA,ACPP
FXR/RXR Activation	2.39E00	3.15E-02	NR0B2,NR1H4,SLC51B,APOD
Atherosclerosis Signaling	1.6E00	2.44E-02	CCR3,RARRES3,APOD
Phospholipases	1.46E00	3.51E-02	RARRES3,PLA1A
Serotonin Degradation	1.41E00	3.28E-02	ALDH3A1,UGT2B15
IL-17 A Signaling in Airway Cells	1.37E00	3.12E-02	MUC5AC,MUC5B
Histamine Degradation	1.19E00	7.69E-02	ALDH3A1
Extrinsic Prothrombin Activation Pathway	1.1E00	6.25E-02	F7
Fatty Acid α-oxidation	1.1E00	6.25E-02	ALDH3A1
Glutathione Redox Reactions I	1.05E00	5.56E-02	GPX2

Note: P-value (-log) was obtained by comparing the expression level of all genes with and without TGF-β treatment.

Table 3

A selection of modulated ‘canonical pathways’ from the H358 lung adenocarcinoma cell line.

PATHWAYS	P-VALUE(-LOG)	RATIO	GENES INVOLVED
Upregulated
Axonal Guidance Signaling	1.01E01	9.49E-02	ADAMTS7,PDGFA,EPHB2,FZD1,NGF,NTNG1,EFNB2,WNT7A,RHOD,ABLIM3,MRAS,ADAM19,WNT4,FZD2,RTN4R,ITGA4,EFNA2,ITGB1,TUBB3,ADAMTS1,RRAS,WNT9A,ITGA2,TUBA4A,ITGA5,VEGFC,L1CAM,MMP2,PDGFB,MYL1,MYL9,TUBA1A,ADAM12,GLIS2,TUBB6,NFATC2,SEMA3C,GLI1,MMP9,WNT11,SEMA7A
Hepatic Fibrosis/Hepatic Stellate Cell Activation	6.92E00	1.12E-01	CTGF,FN1,PDGFA,FGFR1,KLF6,SMAD7,VEGFC,MMP2,COL17A1,COL15A1,MYL1,PDGFB,FGF1,MYL9,COL16A1,COL1A1,COL6A3,TNFSF9,SERPINE1,COL9A2,MMP9,TIMP2
Regulation of the Epithelial-Mesenchymal Transition Pathway	4.96E00	9.78E-02	LOX,SNAI2,RRAS,mir-8,WNT9A,FGFR1,SNAI1,DVL1,MMP2,FZD1,FGF1,CDH2,WNT7A,MRAS,WNT4,FZD2,MMP9,WNT11
Human Embryonic Stem Cell Pluripotency	4.92E00	1.12E-01	PDGFA,FGFR1,WNT9A,DVL1,SMAD7,BMPR2,FZD1,NGF,PDGFB,WNT7A,MRAS,WNT4,FZD2,WNT11
ILK Signaling	4.89E00	9.68E-02	ITGB1,SNAI2,FBLIM1,FN1,LIMS2,SNAI1,VEGFC,VIM,MYL1,MYL9,TGFB1I1,JUN,RHOD,RPS6KA4,ITGB4,ITGB6,MMP9,ACTN1
Leukocyte Extravasation Signaling	4.52E00	9.09E-02	ITGB1,MMP28,ITGA2,THY1,ITGA5,MMP2,CLDN6,MMP23B,RASGRP1,NCF2,CD44,CLDN14,MMP9,ACTN1,ITGA4,TIMP2,MSN,ITK
Epithelial Adherens Junction Signaling	4.48E00	1.03E-01	TUBB3,SNAI2,RRAS,FGFR1,SNAI1,TUBA4A,BMPR2,MYL1,FGF1,MYL9,CDH2,TUBA1A,TUBB6,MRAS,ACTN1
p53 Signaling	4.44E00	1.22E-01	CDKN2A,JUN,GADD45B,SNAI2,GADD45G,TP73,THBS1,CDKN1A,SERPINB5,SFN,SERPINE2,TP53I3
Integrin Signaling	4.41E00	8.91E-02	ITGB1,RRAS,MYLK2,ITGA2,TSPAN2,ITGA5,MYLK,TNK2,PDGFB,MYL9,TLN2,RHOD,MRAS,ITGB4,ITGB6,NEDD9,ACTN1,ITGA4
Downregulated
Cell Cycle Control of Chromosomal Replication	7.31E00	3.33E-01	MCM5,MCM3,MCM6,MCM2,CDC6,CHEK2,MCM4,DBF4,RPA2
Mitotic Roles of Polo-Like Kinase	5.46E00	1.67E-01	PLK4,ESPL1,CDC20,PTTG1,PKMYT1,FBXO5,CDK1,CHEK2,KIF11,CDC25A,CCNB1
Role of CHK Proteins in Cell Cycle Checkpoint Control	4.5E00	1.64E-01	PCNA,RFC2,RFC5,CDK1,CHEK2,E2F2,CDC25A,CHEK1,RFC3
GADD45 Signaling	3.73E00	2.63E-01	PCNA,CCNE2,CCNE1,CDK1,CCNB1
Cyclins and Cell Cycle Regulation	2.66E00	1.03E-01	CCNA2,CCNE2,CCNE1,CDKN2C,CDK1,E2F2,CDC25A,CCNB1
Cell Cycle: G2/M DNA Damage Checkpoint Regulation	2.5E00	1.22E-01	TOP2A,PKMYT1,CDK1,CHEK2,CHEK1,CCNB1
Granulocyte Adhesion and Diapedesis	2.17E00	6.78E-02	CXCL8,CXCL3,SELL,SELE,MMP7,MMP20,ITGAM,NGFR,CCL22,CLDN2,CX3CL1,CLDN3
Nicotine Degradation III	1.76E00	9.8E-02	CYP2F1,UGT2B17,CYP4B1,CYP2S1,UGT2B15
Wnt/β-catenin Signaling	8.83E-01	4.73E-02	SOX2,MMP7,SFRP2,NR5A2,FZD5,WNT8B,SOX5

Note: P-value (-log) was obtained by comparing the expression level of all genes with and without TGF-β treatment.

Network analysis for cancer-related pathways

The up- and downregulated genes used for pathway analysis were also used for network analysis to find biological networks associated with those genes. Network analysis showed strong association of modulated gene expressions with known cancer networks, represented by high network scores. These scores are negative log-transformed P-values of the probability that a particular network is associated with the focus genes by chance. Table 4 lists a selection of five cancer-related gene networks modulated in these two cell lines, with three upregulated and the other two downregulated. The network of pathways involved in cellular movement, organismal injury and abnormalities, and cancer was observed to be significantly upregulated with TGF-β treatment (P-value = 5.82 × 10−11), as shown in Figure 4. A total of 36 unique genes were involved in this network (Table 4), with most of them upregulated in both cell lines after TGF-β treatment. Similarly, the network of pathways related to cancer, tissue development, and hematological disease was also upregulated (P = 9.54 × 10−7), with 35 genes involved (Table 4 and Fig. 5). There were 36 genes involved in the downregulated network of dermatological diseases and conditions, cancer, and neurological disease (P-value = 2.84E-14), and 35 genes in the network of cellular growth and proliferation, tissue development, and organ morphology, which was also observed to be downregulated (P-value = 3.81 × 10−6) in the two NSCLC cell lines after TGF-β treatment (Table 4, Figs. 6 and 7). As with the presentation of the findings from the canonical pathway analysis, we also separated the findings of the network analysis, listing the findings for the A549 cell line separately from the H358 cell line in order to permit an appreciation of the differences in events ongoing in the two cell lines. For the A549 cells, the networks commonly known to promote adoption of an embryonic phenotype and angiogenesis were noted in addition to those typically associated with an EMT; however, the inactivation of systems controlling cellular morphology and cellular adhesion was also observed (Supplementary Table 6). The H358 cells appeared to have networks activated that similarly promoted cellular locomotion and angiogenesis, and also display the inhibition of cellular proliferation, which is commonly associated with EMT (Supplementary Table 7).

Table 4

Selection of gene networks modulated by TGF-β in both lung adenocarcinoma cell lines.

TOP DISEASES AND FUNCTIONS	SCORE	FOCUS MOLECULES	GENES INVOLVED*
Upregulated
Cellular Movement, Organismal Injury and Abnormalities, Developmental Disorder	36	18	Alp,Alpha catenin,CDH4,CGb (includes others),COL1A1, COL6A3,collagen,Collagen Alpha1,Collagen type I,Collagen type III,Collagen type IV, Collagen(s),ERK1/2,FOXS1,ITGB6, JINK1/2,Laminin, LTBP2,mir-27,MMP2,NES,Pdgf (complex),PDGF BB, PDGFB, PMEPA1,RASGRP1,SEMA7A,SERPINE1,Smad,Smad2/3, Smad2/3-Smad4,SNAI1,Tgf beta,TGFBI,THBS1
Cellular Assembly and Organization, Cellular Function and Maintenance, Neurological Disease	33	17	ACKR3,Calmodulin,Cg,CGB7,CGB1/CGB2, CRLF1,EFNA2,estrogen receptor,FLRT2,Focal adhesion kinase, GALNT9,Gpcr,GPR87,GPR115,Jnk,LAMC2, Lh,LPAR5, Mapk,Metalloprotease,Mmp,NPBWR1,P2RY2,P38 MAPK,Pka,Pkc(s),PLEK2,Rac,RAMP1,Sos,STAT,TBXA2R, Tnf (family),TSPAN2,Vegf
Cellular Movement, Cellular Growth and Proliferation, Hematological System Development and Function	31	16	ABLIM3,ADAM19,Akt,Cadherin,caspase,CD3,Cdc2,CDK14, CDKN1A,Cofilin,Cyclin A,Cyclin B,Cyclin E, DOCK2,E2f, F Actin,Hsp27,Ige,IL11,Immunoglobulin, JUNB,MAF,Mek, mir-147,MYL9,PIK3IP1,PPP1R13L, Rb,RBP1,Rock,SEC14L2, TPM1,TSH,Ubiquitin,WNT9A
Downregulated
Cancer, Dermatological Diseases and Conditions, Endocrine System Disorders	45	17	Akt,ALDH3A1,Ap1,APOD,AQP3,Creb,CSF2RA,ERK, ERK1/2, FGFBP1,FSH,GDF15,HGF,Ifn,IL1,INHBB,Insulin,Jnk,LCN2, Mapk,Mek,MUC2,MUC5AC,MUC5B,NFkB (complex),P38 MAPK,Pkc(s), Pro-inflammatory Cytokine, RARRES3,SLPI, TCF,TSPAN8, UGT2B15, Vegf,VTCN1
Cellular Development, Tissue Development, Tissue Morphology	10	5	ACVR1,ACVR2A,AMPH,CHGA,COL4A1,EXOSC2,FAM3C, FAM50A,FNTA,FSHR,FST,FSTL3,GLDN,GOPC,GPRC5A, IL33,JUN,LHCGR,MAN2B2,MAP2K4,MAP2K7,NOG,NRAP, S1PR3,SLC52A2,SLC7A1,SOX2,SOX21, SSTR5, ST6GALNAC1,TAF4B,TGFB1,TNC,UBC,ZDHHC5

Note:

Bold gene names indicate those genes observed modulated in both cell lines in response to TGF-β treatment.

Figure 4

Upregulated network in cellular movement, organismal injury and abnormalities, and cancer. IPA network legend is on the right side.

Figure 5

Upregulated network in cancer, tissue development, and hematological disease. IPA network legend is on the right side.

Figure 6

Downregulated network in dermatological diseases and conditions, cancer, and neurological disease. IPA network legend is on the right side.

Figure 7

Downregulated network in cellular growth and proliferation, tissue development, and organ morphology. IPA network legend is on the right side.

Verification of study findings using reposited microarray data from clinical trials

The significant findings from the analyses of the cell lines were then contrasted against genes, pathways, and networks found to be modulated in a composite of U133a microarray data originating from seven studies and representing 731 patient profiles. These findings are presented primarily in the Supplementary information, given the findings are presented only to corroborate the general findings of the TGF-β induction model. Overall, we found a general agreement of the correlations between the cell line data and the clinical dataset, with 24/137 upregulated and 7/32 downregulated genes being shared between H358 and A549 cell lines with the TGF-β induction. More impressive, however, is the agreement between significant gene expression changes between the clinical dataset (Supplementary Tables 8 and 9) and either H358 or A549 findings (Supplementary Tables 1–4). Further, the agreement at the level of canonical pathways between the datasets was very striking, with “Hepatic Fibrosis/Hepatic Stellate Cell Activation,” “Axonal Guidance Signaling,” and “ILK Signaling” being among the top hits in both the TGF-β model cell line data (Supplementary Table 5) and the clinical datasets (Supplementary Table 10). With further examination of these findings, other examples of pathway modulations are extremely supportive of parallels between the TGF-β model for EMT and signatures of metastatic progression in clinical specimens, such as “Regulation of the Epithelial–Mesenchymal Transition Pathway” and “Integrin Signaling” in the clinical datasets and “TGF-β signaling” and “HIF1α Signaling” in the cell line-based analyses.

Discussion

In this study, we report the modulations of multiple genes and biological pathways that were up- or downregulated in two lung adenocarcinoma cell lines after induction by TGF-β treatment. This investigation was accomplished as a means to identify features that further validate this model of tumor progression, identify potentially underappreciated processes and/or effector molecules intimate to EMT for future study, and reveal novel candidate biomarkers that may have translational significance. Our results revealed a reasonable overlap in upregulated and downregulated significant genes in the two cell lines 358 and 549, indicating that the two cell lines share similar genetic characteristics and similar response to TGF treatment, and these genes may play similar roles in the pathological process. However, a large portion of the up- and downregulated genes did not overlap in the two cell lines (Fig. 2), which could be explained by the high heterogeneity of tumorigenesis in the two cell lines.

Although the “process of metastasis” is an idea well implanted in the mindsets of most cancer researchers, metastasis is far from being well understood on a molecular level.33 Features commonly associated with metastasis for many scientists include the following: increased tumor cell motility, altered cellular adhesion molecules, cell cycle arrest, the ability to remodel extracellular matrices (esp. basement membranes) for extravasation and intravasation, resistance to apoptosis, and the ability to thrive in the secondary site.4–8 One of the more important aspects of this study was to validate that the TGF-β induction model for EMT was consistent with the hallmarks of this process. Indeed, at the individual gene level, we observed upregulation of many of the molecules associated with EMT, including vimentin, N-cadherin, fibronectin-1, Snai1 (“snail”), Snai2 (“slug”), and MMP-2 and -9 just to list a few, and observed the downregulation of cadherin 13, MUC2, MUC5AC, cytokeratin-4, -20, and -40, and collagen type IX (as listed in Supplementary Tables 1–4). At the level of the canonical pathways and network analysis, we further confirm the validity of this model with functional pathways associated with cellular locomotion, ability to remodel basement membranes, and altered cellular signaling (including P53) pathways all activated, whereas pathways associated with epithelial cell signaling and cell cycle control (leading to decreased proliferation) were silenced. These finding are all consistent with our current understanding of EMT. Of interest would be future studies examining the changes in gene expression and pathway/network activation with treatment time. The findings reported here were after 3 days of incubation with TGF-β; it is unclear if these changes would continue in a manner that further supports our understanding of EMT or if other factors, such as IGF-1, need to be recruited to advance to the next phase of this process. Along these lines, we observe an over 29-fold upregulation of IGF-1 in the A549 cells as well as other components of the “IGF-axis” (eg IGFBP-5 and -7) in both cell lines, supporting the idea that IGF-1 is a key mediator in EMT in vitro34–37 and tumor progression in vivo.38–41

Multiple molecular signaling pathways have been extensively investigated in lung cancer,12 and the main pathways observed here are believed to be able to provide roadmaps for therapy of this disease, which include: EGFR/Ras/PI3K growth promoting pathways,42–45 p53/Rb/P14ARF growth inhibitory pathways,46 and Bcl-2/Bax/Fas/FasL apoptotic pathways.47,48 Our observations in this study further provided valuable clues to understand these pathways in lung cancer cells under the circumstance of TGF transformation. For example, we found that the p53-signaling pathway was upregulated via TGF stimuli in these two cell lines (Table 2). Our data also suggested some novel pathways related to tumor cell’s response to TGF treatment, which have not been well investigated in lung cancer yet.

In this study, we also provided results from a network analysis on multiple pathways interaction in the two lung cancer cell lines after TGF-β treatment. We demonstrated four cancer- related networks significantly up- or -downregulated in the two lung cancer cell lines after TGF-β treatment. These networks include multiple interacting genes and reciprocal pathways, which function in various biological processes including cellular movement, growth and proliferation, organismal injury and abnormalities, tissue development, hematologi-cal and dermatological disorders, and cancer development (Figs. 4–7). These observations shed light on the understanding of the complex mechanism of lung cancer development as well as the potential application of these pathways/networks in the treatment of this disease, while further studies are warranted in the future. Further, the parallels with data we present from the analysis of 731 U133a expression microarray profiles are substantial and further legitimize the use of the TGF-β model of EMT induction in the A549 and H358 cell lines for future efforts.

Conclusion

RNA-seq is a useful tool along with statistical, pathway, and network analyses to identify cancer mechanisms. The findings of this study will help direct cell-based studies in our lab modeling early events in lung cancer progression and have also exposed a range of new candidate biomarkers to be explored for potential prognostic value. Our overall objective is to improve treatment options and selection methods as a means to one day improve clinical outcomes for this dreaded disease.

Supplementary Table 1. Selection of 50 genes upregulated in the A549 lung adenocarcinoma cell line.

Supplementary Table 2. Selection of 50 genes downregulated in the A549 lung adenocarcinoma cells.

Supplementary Table 3. Selection of 50 genes upregulated in the H358 lung adenocarcinoma cell line.

Supplementary Table 4. Selection of 50 genes downregulated in the H358 lung adenocarcinoma cells.

Supplementary Table 5. Top 10 canonical pathway modulations for overlapping genes.

Supplementary Table 6. Selection of gene networks modulated by TGF-β in the A549 adenocarcinoma cell lines.

Supplementary Table 7. Selection of gene networks modulated by TGF-β in the H358 adenocarcinoma cell lines.

Supplementary Table 8. A selection of upregulated genes from clinical data set analysis.

Supplementary Table 9. A selection of downregulated genes from clinical data set analysis.

Supplementary Table 10. A selection of modulated canonical pathways from clinical data sets.

Supplementary Figure 1. Impact of TGF-β induction on apoptosis markers. Evaluation of markers for the induction of apoptosis by western blot analysis in both control and TGF-β induced cultures over a 5-day interval, as indicated. All immunoblots were performed according to manufacturer-recommended protocols and optimized against positive control cultures (not shown).

Supplementary Figure 2. Impact of TGF-β induction on the lung adenocarcinoma cell lines. Cellular proliferation of cell lines was monitored using the MTT Cellular Proliferation Assay Kit (ATCC) using 6,000 cells per well in a 96-well plate over 4 days. All cells were adjusted to 2.5% serum for 24 hours prior to initiation of the experiment. Upon each time point absorbance was measured at 570 nm with a Bio-Tek Powerwave XS with normalization to an empty well (Note: *P < 0.05; **P < 0.001).