RNA sequencing identifies gene expression profile changes associated with β-estradiol treatment in U2OS osteosarcoma cells

Abstract

This study was conducted to identify gene expression profile changes associated with β-estradiol (E2) treatment in U2OS osteosarcoma cells by high-throughput RNA sequencing (RNA-seq). Two U2OS cell samples treated with E2 (15 μmol/L) and two untreated control U2OS cell samples were subjected to RNA-seq. Differentially expressed genes (DEGs) between the groups were identified, and main biological process enrichment was performed using gene ontology (GO) analysis. A protein–protein interaction (PPI) network was constructed using Cytoscape based on the Human Protein Reference Database. Finally, NFKB1 expression was confirmed by quantitative real-time polymerase chain reaction (qRT-PCR). The map ratios of the four sequenced samples were >65%. In total, 128 upregulated and 92 downregulated DEGs were identified in E2 samples. After GO enrichment, the downregulated DEGs, such as AKT1, were found to be mainly enriched in cell cycle processes, whereas the upregulated DEGs, such as NFKB1, were involved in the regulation of gene expression. Moreover, AKT1 (degree =117) and NFKB1 (degree =72) were key nodes with the highest degrees in the PPI network. Similarly, the results of qRT-PCR confirmed that E2 upregulated NFKB1 expression. The results suggest that E2 upregulates the expression of NFKB1, ATF7IP, and HDAC5, all of which are involved in the regulation of gene expression and transcription, but downregulates that of TCF7L2, ALCAM, and AKT, which are involved in Wnt receptor signaling through β-catenin and morphogenesis in U2OS osteosarcoma cells.

Keywords:

• differentially expressed genes
• Wnt receptor signaling
• β-catenin
• protein-protein interaction network

Introduction

17β-Estradiol (E2) is a primary sex hormone in human beings that is essential for the development and maintenance of female reproductive organs.¹ However, it also has important effects on many other tissues, such as bone,² liver,³ and brain.⁴ E2 is mainly produced by the granulosa cells of ovaries in women,⁵ but it can also be produced by the testes in men.^6,7 E2 has been used for treating menopausal syndrome⁸ and preventing osteoporosis in postmenopausal women.⁹ In addition, it has been reported that older men with total E2 deficiency are more likely to be osteoporotic.¹⁰ E2 has also been implicated in cancer progression.^11,12 Recently, Tchafa et al¹³ have found that E2 promotes the cellular invasion and proliferation of breast cancer cells. Gunter et al¹² have reported that endogenous E2 levels are positively associated with the risk of colorectal cancer.

Osteosarcoma is the most common primary malignancy of bone and exhibits a high risk of metastasis and poor prognosis.^14,15 Although E2 is known to play a critical role in osteosarcoma, its effects in this disease are controversial. E2 can inhibit purine metabolic and biosynthetic pathways in human osteosarcoma cells to achieve an antagonistic effect on cell proliferation.¹⁶ Previous studies have also shown that E2 protects against cell death in estrogen receptor-α and -β-expressing human U2OS osteosarcoma cells.^17,18 Furthermore, 2-methoxyestradiol, a mammalian metabolite of E2, has been reported to induce cell cycle arrest and osteosarcoma cell apoptosis.¹⁹ Therefore, defining the molecular mechanism(s) of E2 actions in osteosarcoma cells is necessary. In the present study, high-throughput RNA sequencing (RNA-seq) and bioinformatics methods were used to identify changes in the gene expression profile that are associated with E2 treatment of U2OS osteosarcoma cells.

Materials and methods

Cell lines and culture conditions

Human U2OS osteosarcoma cells were purchased from Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, People’s Republic of China). Cells were maintained in phenol red-free Roswell Park Memorial Institute (RPMI)-1640 medium supplemented with 10% fetal bovine serum and 1% penicillin–streptomycin (Thermo Fisher Scientific, Waltham, MA, USA) in a humidified atmosphere of 95% air and 5% CO₂ at 37°C. Before beginning treatments, cells were washed twice with phosphate-buffered saline (PBS) to remove residual serum and were grown in serum-free RPMI-1640 medium for 24 h. Subsequently, E2 (15 μmol/L, dissolved in dimethyl sulfoxide [DMSO]) was added to the medium, and an equal volume of DMSO was added to the control U2OS cell medium. Two U2OS cell samples (14710C-3 and 14710C-4) treated with E2 and two untreated control U2OS samples (14710C-1 and 14710C-2) were subjected to RNA-seq.

RNA extraction and sequencing

After an incubation period of 48 h, cells were washed twice with PBS and harvested. Total RNA was isolated from cultured cells using TRIzol^® Reagent (Thermo Fisher Scientific) according to the manufacturer’s instruction. RNA quality and quantity were assessed by an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). cDNA libraries were constructed using an NEBNext^® Ultra™ RNA Library Prep Kit (Illumina, Shanghai, People’s Republic of China) following the manufacturer’s instruction. Subsequently, libraries were sequenced on Illumina HiSeq 2000 at Beijing Berry Genomics Co., Ltd. Sequenced reads were generated by base calling using the Illumina standard pipeline. Paired-end RNA-seq data were generated with a read length of 100 bp. The raw sequencing data have been uploaded to the National Center for Biotechnology Information database under the BioProject accession no SRP101761.

Alignment of sequenced reads

Raw reads were first filtered to obtain clean reads using FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/).²⁰ High-quality reads were mapped to the human reference genome hg19 using TopHat (version 2.0.12) software²¹ with default parameters. Alignment was independently performed for reads from each sample, and reads mapping to more than three genomic sites were discarded.

Differentially expressed genes (DEGs) screening

Based on the value of reads per kilobase per million mapped reads, gene expression levels were determined using Cufflinks software (version 2.21).²² Subsequently, the Cuffdiff program²³ of Cufflinks was used to identify DEGs. Only the genes with |log(fold change)| >1 and P<0.01 were considered as DEGs.

Gene ontology (GO) enrichment analysis

Database for Annotation, Visualization and Integrated Discovery (DAVID)²⁴ is an online tool used for functional annotation of genes. A GO functional enrichment analysis of DEGs was performed using the DAVID. P<0.05 was chosen as the cutoff criterion.

Construction of protein–protein interaction (PPI) network

The Human Protein Reference Database (http://www.hprd.org/)²⁵ is a database for experimentally derived information about the human proteome, including that on PPIs, post-translational modifications, and tissue expression. DEG-encoding proteins were mapped to the Human Protein Reference Database to search for interaction relationships. The PPI network was visualized using Cytoscape.²⁶ Finally, the hub nodes with a high degree of connectivity²⁷ in the PPI network were also identified.

NFKB1 expression levels using quantitative real-time polymerase chain reaction (qRT-PCR)

To confirm RNA-seq results, NFKB1 and AKT1 expression levels were detected using qRT-PCR. Total RNA was extracted using TRIzol Reagent following the manufacturer’s instructions (TaKaRa, Dalian, People’s Republic of China). Subsequently, the first-strand cDNA was prepared from total lens RNA using a TaKaRa PrimeScript II First Strand cDNA Synthesis Kit (RR036A-1; TaKaRa) according to the manufacturer’s instructions. Glyceraldehyde-3-phosphate dehydrogenase was used as a control. The primers used for NFKB1, AKT1, and glyceraldehyde-3-phosphate dehydrogenase were based on the rat sequences: 5′-AACAGCAGATGGCCCATA CC-3′ (forward), 5′-AACCTTTGCTGGTCCCACAT-3′ (reverse); 5′-GCCTGTCAGCTGGTGCAT-3′ (forward), 5′-CCGCCAGGTCTTGATGTACT-3′ (reverse); and 5′-CA GTGCCAGCCTCGTCTCAT-3′ (forward), 5′-AGGGGCCATCCACAGTCTTC-3′ (reverse), respectively.

Statistical analysis

Differences between the two treatment groups were analyzed using unpaired Student’s t-test. Data analysis was completed using SPSS 22.0 (IBM Corporation, Armonk, NY, USA). P<0.05 was considered to indicate statistically significant difference.

Results

Sequence alignment

The results of TopHat alignment of clean reads with the human reference genome are shown in Table 1. In total, 9,485,360 (68.63%) and 7,843,445 (65.76%) of clean reads were mapped to the human reference genome for the two E2-treated U2OS cell samples; 8,129,145 (69.71%) and 8,531,536 (69.07%) cleans reads were mapped for the two control samples.

Table 1 Summary of clean reads alignment to the reference genome

Sample title	Treatment	Clean reads	Mapped reads	Map ratio (%)
14710C-1	Control	11660638	8129145	69.71
14710C-2	Control	12351701	8531536	69.07
14710C-3	E2	13820514	9485360	68.63
14710C-4	E2	11927581	7843445	65.76

Abbreviation: E2, estradiol.

DEGs and GO enrichment analysis

In total, 220 genes, including 128 upregulated and 92 downregulated genes, were identified as being significantly differently expressed between E2-treated U2OS osteosarcoma cells and the controls. According to the GO enrichment analysis, the top five GO terms of upregulated DEGs NFKB1, ATF7IP, HDAC5, MEN1, and EPC1 were significantly related to the regulation of gene expression and transcription (Table 2). Meanwhile, the top five GO terms of downregulated DEGs included Wnt receptor signaling through β-catenin involving RARG, TBL1X, and TCF7L2; axonogenesis involving ALCAM, NRP1, and SLC26A6; and cell cycle processes involving AKT1, DSN1, and POLD1 (Table 3).

Table 2 The top five GO terms of upregulated DEGs

Category	Term	Description	P-value	Genes
BP	GO:0010628	Positive regulation of gene expression	0.002272	NFKB1, ATF7IP, HDAC5, etc
BP	GO:0010629	Negative regulation of gene expression	0.002581	NFKB1, ATF7IP, HDAC5, etc
BP	GO:0016568	Chromatin modification	0.002959	HDAC5, MEN1, EPC1, etc
BP	GO:0045944	Positive regulation of transcription from RNA polymerase II promoter	0.004248	NFKB1, HDAC5, MEN1, etc
BP	GO:0016481	Negative regulation of transcription	0.004572	NFKB1, ATF7IP, HDAC5, etc

Abbreviations: BP, biological process; GO, gene ontology; DEGs, differentially expressed genes.

Table 3 The top five GO terms of downregulated DEGs

Category	Term	Description	P-value	Genes
BP	GO:0060070	Wnt receptor signaling pathway through β-catenin	0.002951	RARG, TBL1X, TCF7L2
BP	GO:0007409	Axonogenesis	0.01382	ALCAM, NRP1, SLC26A6, etc
BP	GO:0048667	Cell morphogenesis involved in neuron differentiation	0.01801	ALCAM, NRP1, SLC26A6, etc
BP	GO:0022402	Cell cycle process	0.01832	AKT1, DSN1, POLD1, etc
BP	GO:0048812	Neuron projection morphogenesis	0.01916	ALCAM, NRP1, SLC26A6, etc

Abbreviations: BP, biological process; GO, gene ontology; DEGs, differentially expressed genes.

Construction of the PPI network

In total, 1,185 nodes, including 91 upregulated DEGs, 55 downregulated DEGs, and 1,064 non-DEGs, were present in the PPI network (Figure 1). The top five DEGs with the highest degree of connectivity in the network were AKT1 (117), NFKB1 (72), ATF7IP (64), NCOA3 (45), and HDAC5 (36).

Figure 1 PPI network of DEGs.

Note: Red nodes represent upregulated DEGs, green nodes represent downregulated DEGs, and blue nodes represent non-DEGs.
Abbreviations: PPI, protein–protein interaction; DEGs, differentially expressed genes.

Differences in NFKB1 and AKT1 expression levels

As shown in Figure 2A, NFKB1 expression levels significantly increased when U2OS osteosarcoma cells were treated with E2 (P=0.002), which confirmed the reliability of the bioinformatics method. Although AKT1 expression levels were increased when U2OS osteosarcoma cells were treated with E2, they were not significantly different (P>0.05; Figure 2B).

Figure 2 Related NFKB1 (A) and AKT1 (B) expression levels in U2OS osteosarcoma cells treated with E2.

Note: **P<0.01.
Abbreviations: E2, estradiol; DMSO, dimethyl sulfoxide.

Discussion

In the current study, RNA-seq was used to explore changes in the gene expression profile that are associated with E2 treatment of U2OS osteosarcoma cells. We found that E2 treatment induced the upregulation of genes related to the regulation of gene expression and transcription (eg, NFKB1, ATF7IP, and HDAC5) and downregulation of those involved in Wnt receptor signaling through β-catenin and morphogenesis (eg, TCF7L2, ALCAM, NRP1, SLC26A6, and AKT).

Our results demonstrated that NFKB1 was mainly enriched in the regulation of gene expression and transcription. NFKB1 belongs to the NF-κB family, which contains a group of proteins involved in carcinogenesis, immune response, cell adhesion, proliferation, angiogenesis, and apoptosis.²⁸ NF-κB is a transcription factor that participates in the regulation of viral and cellular genes.²⁹ Constitutive NF-κB activation has been observed in 67% of colorectal cancer cell lines and promoted tumor growth.³⁰ NFKB2, another member of the NF-κB family, can stimulate cell proliferation in U2OS osteosarcoma cells.³¹ Furthermore, several studies have shown the relationship between NFKB1 and tumors. For example, a functional insertion/deletion polymorphism in the promoter region of NFKB1 increases the risk of nasopharyngeal carcinoma.³² Riemann et al³³ discovered that the NFKB1 promoter polymorphism was a useful molecular marker for the risk of recurrence in superficial bladder cancer. However, although most studies have analyzed the correlation of the NFKB1 promoter polymorphism with tumors,^31,34 those analyzing the effects of NFKB1 expression in osteosarcoma are limited. Thus, further experiments are needed to explore whether NFKB1 expression has any impact on osteosarcoma progression.

We also found that E2 treatment may repress the expression of genes, such as RARG, TBL1X, and TCF7L2, involved in Wnt receptor signaling through β-catenin in U2OS osteosarcoma cells. The activation of Wnt signaling and the accumulation of β-catenin have been reported in many carcinomas,^35,36 including osteosarcoma.³⁷ TCF7L2 encodes the transcription factor TCF-4, which can be activated by dephosphorylated β-catenin via binding to a conserved N-terminal region in the nucleus, thereby initiating the expression of target genes, including the proto-oncogenes c-jun and fra-1.³⁸ Thus, E2 treatment may have an unfavorable effect on U2OS osteosarcoma cells.

Among the other downregulated DEGs, AKT1 was observed to have the highest degree in the PPI network. AKT, also known as protein kinase B (PKB; a serine/threonine kinase), is one of the most critical and versatile protein kinases involved in the mechanism of human physiology and disease.^39,40 The activation of AKT pathways plays a central role in tumor metastasis.⁴¹ Furthermore, Fukaya et al⁴² have demonstrated the important role of AKT signaling in the pulmonary metastasis of osteosarcoma. The AKT family has three members: AKT1/PKBα, AKT2/PKBβ, and AKT3/PKBγ. AKT1 and AKT2 have been reported to be ubiquitously and similarly expressed in various tissues.⁴³ Recently, Zhu et al⁴⁴ have reported that elevated AKT2 expression is associated with poor outcomes in patients with osteosarcomas.

Here, AKT1 was enriched in the cell cycle process, which was closely related to tumor progression. A study by Ju et al⁴⁵ has revealed that AKT1 governed breast cancer progression in mice, whereas another study has indicated that AKT1 amplification regulates cisplatin (a chemotherapeutic agent) resistance in human lung cancer.⁴⁶ Collectively, these reports have demonstrated that AKT1 might be closely involved in osteosarcoma metastasis. In the current study, we predicted that AKT1 expression was downregulated. According to qRT-PCR results, although AKT1 expression levels were increased when U2OS osteosarcoma cells were treated with E2, they were not significantly different. However, the effect of AKT expression on U2OS osteosarcoma cells is controversial. For example, Nielsen-Preiss et al⁴⁷ have reported that the downregulation of AKT expression enhances osteosarcoma cell proliferation, whereas Díaz-Montero et al⁴⁸ have found that AKT expression is upregulated in anoikis-resistant human osteosarcoma SAOSar cells. Therefore, we speculated that E2 is involved in osteosarcoma metastasis, but the modulating mechanism is still unclear.

ALCAM was another gene downregulated by E2 treatment in U2OS osteosarcoma cells, which was speculated to function in axonogenesis, cell morphogenesis involved in neuron differentiation, cell cycle processes, and neuron projection morphogenesis. ALCAM encodes the CD166 antigen, which is a 100–105 kDa type-I transmembrane glycoprotein of the immunoglobulin protein superfamily. Similar to CD29, CD44, CD73, CD90, CD105, and CD106, it is known as a marker of the mesenchymal stem cell phenotype.^49,50 Its overexpression has also been reported in colorectal carcinoma⁵¹ and has been demonstrated to be associated with a poor prognosis of several tumors.^51–53 Federman et al have reported that CD166 is highly expressed in the osteosarcoma cell lines HOS, KHOS, KHOS240s, and SJSA, but its expression status is not observed in U2OS osteosarcoma cells.⁵⁴ They further proposed that this gene is a potential candidate for the targeted therapy of osteosarcoma. Our findings suggested that E2 treatment decreases ALCAM expression in U2OS osteosarcoma cells, thereby playing an inhibitory role against osteosarcoma.

Taken together, we found that E2 treatment may mainly upregulate the expression of genes, such as NFKB1, ATF7IP, and HDAC5, related to the regulation of gene expression and transcription and downregulate that of genes involved in Wnt receptor signaling through β-catenin and morphogenesis (eg, TCF7L2, ALCAM, NRP1, SLC26A6, and AKT) in U2OS osteosarcoma cells. Thus, we proposed that E2 has an unfavorable effect against U2OS osteosarcoma cells. However, given that our findings were partly obtained using bioinformatics tools, they need to be further validated.

Acknowledgments

This study was funded by a special fund for the medical service of the Jilin Finance Department (Grant No SCZSY201507).

Disclosure

The authors report no conflicts of interest in this work.