Pivotal regulatory network and genes in osteosarcoma.

INTRODUCTION: Understanding the transcriptional regulatory networks that map out the coordinated responses of transcription factors and target genes would represent a significant advance in the analysis of osteosarcoma, a common primary bone malignancy. The objective of our study was to interpret the mechanisms of osteosarcoma through the regulation network construction.
MATERIAL AND METHODS: Using GSE14359 datasets downloaded from Gene Expression Omnibus data, we first screened the differentially expressed genes in osteosarcoma. We explored the regulation relationship between transcription factors and target genes using Cytoscape. The underlying molecular mechanisms of these crucial target genes were investigated by Gene Ontology function and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis.
RESULTS: A total of 1836 differentially expressed were identified and 98 regulatory relationships were constructed between 32 transcription factors and their 60 differentially expressed target genes. Furthermore, BCL2-like 1 (BCL2L1), tumor protein p53 (TP53), v-rel reticuloendotheliosis viral oncogene homolog A (avian) (RELA), interleukin 6 (IL6), retinoic acid receptor, alpha (RARA), nuclear factor I/C (CCAAT-binding transcription factor) (NFIC), and CCAAT/enhancer binding protein, beta (CEBPB) formed a small pivotal network, in which IL-6 could be regulated by TP53, NFIC, RARA, and CEBPB, but BCL2L1 may be only regulated by TP53 and RELA. These genes had been demonstrated to be involved in osteosarcoma progression via various biological processes and pathways, including regulation of cell apoptosis, proliferation, antigen processing and presentation pathway, and phosphatidylinositol signaling system.
CONCLUSIONS: In general, we have obtained a regulatory network and several pathways that may play important roles in osteosarcoma, identified several pivotal genes in osteosarcoma, and predicted several potential key genes for osteosarcoma.

Introduction

Osteosarcoma (OS) is a relatively uncommon cancer but the most common primary bone malignancy in both children and young adults. The overall relapse-free survival rate over 5 years remains approximately 65% despite modern treatment protocols that combine chemotherapy, surgery, and radiotherapy. Hence, identifying molecular targets that are specific for OS will be critical to the development of new treatment strategies to improve patient outcomes [1-3].

The clinical manifestation of a cancer is based on six essential alterations in cell physiology [4]: 1) Self-sufficiency in growth signals. Tumor cells have the ability to proliferate in the absence of growth signals such as growth factors. Many oncogenes can mimic normal growth signaling, such as FBJ murine osteosarcoma viral oncogene homolog (c-fos) and v-myc myelocytomatosis viral oncogene homolog (avian) (c-myc), which both are demonstrated to be over-expressed in OS [5]. 2) Insensitivity to growth inhibitory signals. Nearly all the anti-proliferative signals converge on the retinoblastoma protein (pRb). Once the pRb-pathway is disrupted by tumor growth factor β (TGF-β), the cell becomes insensitive to antigrowth factors. High-grade OS is found to express TGF-β1 in significantly higher amounts than low-grade OS [6]. 3) Apoptosis evasion. Resistance to apoptosis can be acquired by a mutation in the p53 tumor suppressor gene. The PI3K-AKT/PKB pathway is another way to transmit anti-apoptotic survival signals, activated by insulin-like growth factor 1 and insulin-like growth factor 2 [7]. 4) Limitless replicative potential. Telomere dysfunction might have major implications in tumor progression in patients with OS. 5) Sustained angiogenesis. Angiogenesis can be stimulated by positive signals, such as vascular endothelial growth factor (VEGF), which bind to tyrosine kinase receptors on endothelial cells. The mRNA expression of VEGF in OS has been reported [8]. 6) Tissue invasion. Invasion of the surrounding tissues by OS involves degradation of the extracellular matrix. Matrix metalloproteinases are principally involved in the breakdown of the extracellular matrix [9].

Microarray analysis has been used to screen for gene expression alterations and identify potential targets in human OS cell lines [10, 11]. With the use of genome-wide cDNA microarrays, the transcriptome profile of two OS cell lines (namely, Sa OS and U-2 OS) has been investigated. The author identifies 1,098 differentially regulated spots in Sa OS versus U-2 OS cells including 796 functionally characterized genes [12]. Microarray analysis is also performed to identify histological subtype specific differentially expressed genes (DEGs), namely, osteoblastic and non-osteoblastic OS. The results show that 75 genes are up-regulated and 97 genes are down-regulated in osteoblastic versus non-osteoblastic OS samples, respectively [13].

In this study, we used a similar strategy to identify gene expression profiles that distinguish OS patients from healthy controls. Furthermore, relevant transcription factor (TF) genes, target genes, and pathways in the network were analyzed to explain the potential interaction mechanisms between them in the OS.

Material and methods

Affymetrix microarray data

The transcription profile GSE14359 was obtained from the public functional genomics data repository Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/). Samples were obtained from the Max-Delbrueck-Center, Division of Pathology, Berlin-Buch, Germany. In their study, mRNA from 5 frozen conventional OS samples (including 3 males and 2 females, the age ranging from 7 to 74 years old) and 1 non-neoplastic primary osteoblast sample with 2 replicates, were extracted and hybridized to the Affymetrix Human Genome U133A Arrays. Total 12 microarrays were used to identify the DEGs.

Pathway data

Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) [14] connects known information on molecular interaction networks, such as pathways and complexes, information about genes and proteins generated by genome projects and information about biochemical compounds and reactions. As for OS, a total of 130 pathways and 2287 genes were collected from this database.

Regulation data

The TRANSFAC database contains data on TFs, their experimentally proven binding sites, and regulated genes [15]. The Transcriptional Regulatory Element Database (TRED) has been built in response to increasing needs of an integrated repository for both cis- and trans- regulatory elements in mammals [16]. The TRED makes the curation for transcriptional regulation information, including TF binding motifs and experimental evidence. The curation is currently focused on target genes of 36 cancer-related TF families. Seven hundred and seventy-four pairs of regulatory relationship between 219 TFs and 265 target genes were collected from TRANSFAC (http://www.gene-regulation.com/pub/databases.html). Five thousand seven hundred and twenty-two pairs of regulatory relationship between 102 TFs and 2920 target genes were collected from TRED (http://rulai.cshl.edu/TRED/). We combined the two regulation datasets, and then in total 6328 regulatory relationships between 276 TFs and 3002 target genes were collected.

Statistical analysis

For the GSE14359, the limma method [17] was used to identify DEGs. The 5 conventional OS samples were compared to 1 non-neoplastic primary osteoblast sample with 2 replicates each. The original expression datasets from all conditions were processed into expression estimates using the robust multiarray average method with the default settings implemented in Bioconductor, and then to construct the linear model. Only DEGs with the fold change > 2 and p-value < 0.05 were considered significant.

One thousand eight hundred and thirty-six genes were selected as DEGs from GSE14359.

To demonstrate potential regulatory relationships, the Pearson correlation coefficient (PCC) was calculated for all pair-wise comparisons of gene-expression values between TFs and DEGs. Regulatory relationships whose absolute PCC were larger than 0.75 were considered significant.

Regulation network construction

Based on the two TRANSFAC and TRED regulation datasets, we built the regulation networks using Cytoscape [18]. Based on the significant relationships (PCC > 0.75 as the threshold) between TFs and their target genes, 98 putative regulatory relationships were predicted between 32 TFs and 60 target genes.

Gene ontology (GO) enrichment

DAVID [19] provides a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind large lists of genes. For any given gene list, DAVID tools are able to identify over-represented GO categories in a biological process.

Pathway analysis

We adopted an impact analysis that not only includes the statistical significance of the set of pathway genes, but also other crucial factors, such as the magnitude of each gene's expression change, the topology of the signaling pathway, and their interactions [20].

In this model, the impact factor (IF) of a pathway (Pi) was calculated as the sum of two terms:

The first term is a probabilistic term that captures the significance of the given pathway Pi from the perspective of the set of genes contained in it.

The second term is a functional term that depends on the identity of the specific genes that are differentially expressed as well as on the interactions described by the pathway (i.e., its topology).

Results

Regulation network construction in OS

To get DEGs of OS, we obtained publicly available GSE14359 microarray data sets from GEO. After microarray analysis, the genes with the fold change value larger than 2 of GSE14359 and p-value less than 0.05 were selected as DEGs. A total of 1836 genes were selected as DEGs from GSE14359. To get the regulatory relationships, the co-expressed value (PCC > 0.75) was chosen as the threshold. Finally, 98 regulatory relationships consisting of 32 TFs and their 60 differently expressed target genes were selected. By integrating the regulatory relationships and the DEGs, we built a regulation network of OS between TFs and their differently expressed target genes (Figure 1), including 98 regulatory relationships, 32 TFs, and 60 differently expressed target genes. But 8 TFs, including E2F transcription factor 4, p107/p130-binding (E2F4), regulatory factor X, 3 (RFX3), estrogen receptor 1 (ESR1), SMAD family member 2 (SMAD2), spleen focus forming virus (SFFV) proviral integration oncogene (SPI1), SP3, nuclear transcription factor Y, alpha (NFYA), and MYC, regulated their target genes independently. And we also found 3 regulationships between TFs: retinoic acid receptor, alpha (RARA), retinoic acid receptor, beta (RARB), and retinoic acid receptor, gamma (RARG), 3 retinoic acid receptor subtypes [21].

Figure 1

Regulation network of osteosarcoma The triangle denotes the transcription factor and the circle denotes target genes. The black line suggests that the transcription factor could activate their target genes in OS. In contrast, the grey line suggests that the transcription factor could inhibit the expression of their target genes in OS

To present the key regulationship for the regulatory network, we screened out genes in the network with connectivity degree less than 5. Finally, the pivotal network of OS was constructed, which would explain the core function and regulation in the OS (Figure 2). In the pivotal network, interleukin 6 (IL-6) and BCL2-like 1 (BCL2L1) were both regulated by tumor protein p53 (TP53) and v-rel reticuloendotheliosis viral oncogene homolog A (avian) (RELA). In addition, IL-6 could also be regulated by nuclear factor I/C (CCAAT-binding transcription factor) (NFIC), RARA, and CCAAT/enhancer binding protein (C/EBP), beta (CEBPB). BCL2L1 may be regulated by NFIC by indirectly influencing TP53.

Figure 2

Pivotal regulation network of osteosarcoma. The triangle denotes the transcription factor and the circle denotes target genes. The black line suggests that the transcription factor could activate its target gene in OS

GO enrichment of the regulation network

To depict the regulation network, the GO enrichment was applied to the regulation network genes. Several biological process categories were significantly enriched, including regulation of apoptosis, regulation of programmed cell death, regulation of cell death, and others (Table I).

Table I

GO enrichment analysis of biological process

Category	Term	Count	Value of p	FDR
BP	GO:0042981 ~ regulation of apoptosis	24	2.03 × 10−13	3.35 × 10−10
BP	GO:0043067 ~ regulation of programmed cell death	24	2.50 × 10−13	4.13 × 10−10
BP	GO:0010941 ~ regulation of cell death	24	2.70 × 10−13	4.46 × 10−10
BP	GO:0042127 ~ regulation of cell proliferation	22	1.17 × 10−11	1.94 × 10−8
BP	GO:0009628 ~ response to abiotic stimulus	15	6.71 × 10−10	1.11 × 10−6
BP	GO:0048545 ~ response to steroid hormone stimulus	12	7.96 × 10−10	1.32 × 10−6
BP	GO:001003 ~ response to organic substance	19	1.47 × 10−9	2.44 × 10−6
BP	GO:0043065 ~ positive regulation of apoptosis	15	5.01 × 10−9	8.28× 10−6
BP	GO:0043068 ~ positive regulation of programmed cell death	15	5.47 × 10−9	9.05 × 10−6
BP	GO:0010942 ~ positive regulation of cell death	15	5.81 × 10−9	9.60 × 10−6
BP	GO:0042981 ~ regulation of apoptosis	24	2.03 × 10−10	3.35 × 10−10
BP	GO:0043067 ~ regulation of programmed cell death	24	2.50 × 10−13	4.13 × 10−10
BP	GO:0010941 ~ regulation of cell death	24	2.70 × 10−13	4.46 × 10−10
BP	GO:0042127 ~ regulation of cell proliferation	22	1.17 × 10−11	1.94 × 10−8

BP – biological process, FDR – false discovery rate

Significant pathway in OS

To further interpret the regulatory network at a higher regulation level and obtain a better biological interpretation, we adopted a pathway based impact analysis of the network. This method yielded many significant pathways, including Antigen processing and presentation, Phosphatidylinositol signaling system, Focal adhesion, Adherent junction, and others (Table II).

Table II

Significant pathway analysis result

Database name	Pathway name	Impact factor	% Pathway genes in input	Corrected gamma value of p
KEGG	Antigen processing and presentation	125.52	19.101	3.89 × 10−53
KEGG	Phosphatidylinositol signaling system	54.672	10.526	1.00 × 10−22
KEGG	Focal adhesion	27.034	27.094	5.09 × 10−11
KEGG	Adherens junction	23.718	24.359	1.24 × 10−9
KEGG	Systemic lupus erythematosus	17.418	22.222	5.02 × 10−7
KEGG	ECM-receptor interaction	15.867	29.762	2.17 × 10−6
KEGG	Prostate cancer	15.421	26.667	3.30 × 10−6
KEGG	Cell cycle	13.954	24.576	1.30 × 10−5
KEGG	Pathways in cancer	13.084	17.879	2.93 × 10−5
KEGG	Melanoma	12.465	22.535	5.20 × 10−5

Discussion

According to the results, we found that many TFs and pathways closely related to OS were identified by our methods. The model provided a global view of TF regulation and response of target genes in OS. Many genes in the regulatory network, including TP53, RELA, NFIC, MYC, IL-6, BCL2L1, RARA and CEBPB, have been identified related to OS based on previous papers. For example, in the regulation network, MYC exhibited a significantly different regulation mode by up-regulating 5 target genes and down-regulating 4 genes. MYC protein is a multifunctional, nuclear phosphoprotein that plays a role in cell cycle progression, apoptosis and cellular transformation. Over-expression of MYC in bone marrow stromal cells leads to OS development and loss of adipogenesis [22]. Additionally, MYC has been examined as a therapeutic target for OS. Down-regulation of MYC enhances the therapeutic activity of methotrexate against OS cells [23]. The results suggest that the approach we used could reliably identify gene co-expression networks.

Further, among the regulation network (Figure 1), BCL2L1, TP53, RELA, IL6, RARA, NFIC, and CEBPB formed a small pivotal network, in which IL-6 could be regulated by TP53, NFIC, RARA, and CEBPB, but BCL2L1 may be only regulated by TP53 and RELA. We will discuss the role of these genes in OS and their interaction relationships as follows based on previous reports.

IL-6 gene encodes a cytokine as one osteoclast differentiating factor involved in osteoclast formation. Active osteoclasts are frequently present in OS, which requires IL-6 mRNA high expression in OS tissue to stimulate osteoclast activity, facilitate OS further invasion, and cause release of pro-resorptive cytokines [24]. The IL6 mRNA up-regulation in OS may be attributed to TP53 mutation because wild-type (wt) human TP53 preferentially represses the IL-6 promoter in HeLa cells [25].

TP53 is postulated to bind to a p53-binding site and activate expression of downstream genes that inhibit growth and/or invasion, and thus function as a tumor suppressor. Mutants of p53 frequently occurred in human OS cells fails to bind the consensus DNA binding site, and hence causes the loss of tumor suppressor activity [26]. And p53 mutation has been shown to be more common in high-grade conventional OS versus low grade central OS [27].

NFIC is a member of the NFI gene family, which plays wide reaching roles in viral DNA replication, regulation of gene transcription, cell proliferation, and development. NFIC is also found expressed in human OS cell lines mediated by IGFBP5 promoter activity [28]. Although our results indicated that NFIC may regulate IL-6 expression directly, no experimental evidence was supported here. NFIC is thought to be a cofactor to regulate the transcription of p53 [29]. And p53 could regulate IL-6 [25]. Thus, NFIC may indirectly regulate the expression of IL-6.

RELA is often bound to NFκB1 to form the p65 (RELA)/p50 (NFκB1) complex and induce the expression of IL-6 [30]. The regulatory role of the p65/p50 subunit in tumor cells shows great diversity. Cisplatin treatment in the U-2 OS cell line represses RelA activity and inhibits expression of the NF-κB antiapoptotic target gene BCL2L1. In contrast, another chemotherapeutic drug, etoposide, could activate NF-κB and induce BCL2L1 gene expression. These observations suggest that it may be possible to minimize the ability of RelA to inhibit OS therapy by diagnostically predicting the type of chemotherapeutic drug [31].

Agonists of several members of the nuclear receptors have been shown to inhibit proliferation and promote differentiation in OS cells. Among them, RARs [α, β, or γ] attract the most attention [32]. Over-expression of RARα has been demonstrated to effectively inhibit OS cell proliferation [33].

CEBPB is a bZIP transcription factor which can bind as a homodimer to certain DNA regulatory regions. C/EBPB expression increases from the growth to maturation developmental stages of osteoblasts. C/EBPB also could activate osteocalcin gene transcription and synergize with runt-related transcription factor 2 (Runx2) at the CEBP element to regulate bone-specific expression in an OS cell line [34]. In addition, CEBPB is downstream of the mammalian target of rapamycin kinase (mTOR), a target of immunosuppressive and anticancer drugs. Therefore, C/EBPB may represent a novel therapeutic approach in OS [35].

Identically, no experimental evidence was provided here to demonstrate that RARα and CEBPB could directly interact with IL-6.

In brief, in the pivotal network, IL-6 could be regulated by 5 TFs, 3 regulationships identified by previous works. This suggest that IL-6 is a pivotal gene in OS.

BCL2L1 protein, also known as BCLXL, belongs to the BCL-2 protein family. The levels of Bcl-xL mRNA expression are significantly higher in OS tissues than corresponding non-tumor tissues. Bcl-xL down-regulation could significantly enhance in vitro chemo- or radiosensitivity of OS cells. Taken together, over-expression of Bcl-xL may play important roles in OS progression and this molecule will be a potential chemo- or radiotherapeutic molecular target for OS therapy [36].

Decreased p53 expression in HNSCC lines is predicted to enhance activation of BCL-XL. Over-expression of p53 induces greater increase in the BAX/ BCL-XL ratio [37].

T cell antigen receptor-independent CD28 signal leads to selective transcription of the survival gene Bcl-xL mediated by the specific recruitment of RelA and p52 NF-κB subunits to target promoters [38].

Taking the above together, NFIC, TP53, RELA, RARA, and CEBPB were supposed to be the pivotal regulators for BCL2L1 and IL-6, which belong to the signature genes for OS. They may be potential molecular targets for OS therapy.

Besides the gene regulationship, several significantly correlated OS pathways were also found. For example, the phosphatidylinositol 3-kinase (PI3K)/AKT pathway plays an important role in various cellular processes including cell growth, survival, and motility. Recently, accumulating evidence has indicated that PI3K/Akt is involved in anoikis resistance of human OS cells SAOS-2. The pharmacological inhibition of PI3-K activity is correlated with an increase in anoikis among initially resistant SAOSar cells mediated by Src activation [39]. Hence, abnormal function of the PI3K/AKT pathway might be a potential target for OS chemotherapy [40].

There is also evidence that adherens junction is related to OS formation. A study has demonstrated that osteoblasts depend on both pRb and cell-to-cell contacts for their differentiation and function. RB knock-out mice shows loss of adherens junction at osteoblast membranes. Further, as a well-known regulator of adherens junction assembly, merlin is inactivated in osteosarcoma cells deficient in pRb. Therefore, pRb inactivation results in adherens junction defect in osteoblasts and then contributes to OS formation [41].

In addition, the antigen processing and presentation pathway is also suggested to be involved in OS progression. HLA-C, HLA-DOA, HLA-DPB1, HLA-DPA1, and HLA-E, all belonging to the antigen-presentation pathway, show down-regulated expression in OS lesions compared to non-malignant bone. High-mobility group box 1 (HMGB1) expression is also down-regulated in patients. The interaction of HMGB1 with Toll-like receptor 4 on dendritic cells has been shown to be essential for tumor antigen processing and presentation [42].

In conclusion, we used network analysis as a conceptual framework to explore the pathobiology of OS, based on the assumption that OS is a contextual attribute of distinct patterns of interactions between multiple genes. The salient results of our study included many related TFs, target genes, and pathways, which were correlated with OS directly or indirectly. We suggest that NFIC, TP53, RELA RARA, and CEBPB may be the key TFs for OS therapy. However, further experiments are still indispensable to confirm our conclusion.