Transcriptome analysis reveals the potential biological function of FSCN1 in HeLa cervical cancer cells

Cell culture and transfection

The human cervical cancer cell line HeLa was purchased from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China). Cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) containing 10% fetal bovine serum (FBS). HeLa cells were transfected with the pGFP-B-RS vector with shRNA sequence (CCCTTGCCTTTCAAACTGGAA) inserted to knock down FSCN1. A vector without the shRNA was transfected into HeLa cells as the control. In addition, FSCN1 siRNA (CCCTTGCCTTTCAAACTGGAA) and nontargeting control siRNA (UUCUCCGAACGUGUCACGUTT) were used to transfect HeLa cells to knock down FSCN1. Both the FSCN1 shRNA and siRNA had the same sequence, which targeted exon 5 in the 3′ UTR of FSCN1. Transfection was performed via Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. After 48 h, transfected HeLa cells were collected for qPCR and Western blot analysis of FSCN1 expression.

Assessment of knockdown efficiency

Total RNA was isolated from HeLa cells transfected with different vectors using TRIzol reagent (Ambion, USA). Then, cDNA was synthesized. qPCR was performed using Bestar SYBR Green RT-PCR Master Mix (DBI Bioscience, Shanghai, China) on a Bio-Rad S1000 thermal cycler. GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was used as the reference gene to assess the efficiency of FSCN1 knockdown. The primers used for qPCR are listed in Table S1. The expression of FSCN1 was then normalized to GAPDH using the 2^−ΔΔCT method.

Western blotting

HeLa cells were lysed using RIPA buffer (Beyotime) on ice for 30 min. Then, centrifugation of the sample was performed at 12,000 rpm at 4 °C for 10 min. The supernatant sample was separated by SDS–PAGE and transferred onto polyvinylidene fluoride membranes. The membranes were blocked using nonfat milk and incubated overnight with primary antibodies against FSCN1 (1:1,000, ABclonal), ANGPTL4 (1:1,000, ABclonal), and GAPDH (1:5,000, ATA). Subsequently, the membranes were incubated with secondary antibodies (Proteintech, Wuhan). Bound secondary antibody was detected using enhanced chemiluminescence (ECL) reagent (Bio-Rad, 170506).

Library preparation and sequencing

Total RNA was extracted from HeLa cells with TRIzol (Ambion) and purified by two phenol–chloroform extractions. Then, the purified total RNA was treated with RQ1 DNase (Promega, Madison, WI, USA) to remove DNA. The quality and quantity of further purified RNA was assessed by measurement of the absorbance ratio at 260 nm/280 nm (A260/A280) using a Smartspec Plus spectrophotometer (BioRad, USA). Agarose gel electrophoresis (1.5%) was performed to assess the integrity of the purified RNA.

For RNA-seq library preparation, 1 µg RNA was used for each sample. Oligo(dT)-conjugated magnetic beads (Invitrogen, Carlsbad, CA, USA) were used to purify and concentrate polyA mRNAs from the total RNA sample. RNA in the purified and concentrated sample was fragmented. End-repair and 5′ adaptor ligation of fragmented RNA were performed. Then, the RNA was reverse transcribed using primers harboring a 3′ adaptor sequence and randomized hexamers for cDNA. The cDNAs were further amplified and purified. The purified cDNAs were stored at −80 °C until sequencing.

High-throughput sequencing libraries were prepared according to the manufacturer’s instructions. Paired-end sequencing (151-bp) was conducted with the Illumina HiSeq4000 system at a commercial company (ABLife Inc., Wuhan, China).

Clean and alignment of raw sequencing data

First, raw sequencing reads with more than 2-N bases were discarded. Then, FASTX-Toolkit (Version 0.0.13) was used to trim the adaptors and low-quality bases from the raw reads. The resulting reads less than 16 nt were further discarded. TopHat2 was used to align the clean reads to the GRch38 genome allowing four mismatches (Li et al., 2013). Then, the uniquely mapped reads for each gene were used to calculate the FPKM (paired-end fragments per kilobase of exon per million fragments mapped) values (Trapnell et al., 2010).

Differentially Expressed Genes (DEG) analysis

FPKM was used to evaluate the expression level of genes. edgeR software was used to identify the differentially expressed genes (DEGs) (Robinson, McCarthy & Smyth, 2010) with fold change threshold (fold change ≥2 or ≤0.5) and a false discovery rate threshold (FDR < 0.05).

Functional enrichment analysis of DEGs

Gene Ontology (GO) analysis was used to predict the gene function and calculate the functional category distribution frequency, which annotates an input set of genes with putative pathways and disease relationships based on mapping to genes with known annotations (Xie et al., 2011). A hypergeometric test with Benjamini–Hochberg FDR correction were used to define the enrichment of each pathway (corrected p value < 0.05).

Validation of DEG expression by qPCR

Quantitative real-time PCR (qPCR) of selected DEGs was conducted to elucidate the validity of the RNA-seq data. GAPDH was used as a reference gene. Primers for qPCR are presented in Table S1. qPCR was conducted on the same RNA samples that were used for RNA-seq. The thermal cycling conditions for PCR were as follows: denaturation for 10 min at 95 °C, followed by 40 cycles of including denaturation for 15 s at 95 °C and annealing and extension for 1 min at 60 °C. Three technical replicates were performed for PCR amplification of each sample.

Analysis of TCGA (The Cancer Genome Atlas) data

The RNA-seq data of cervical cancer in the TCGA database were downloaded from UCSC Xena (https://xenabrowser.net/datapages/). The expression levels (FPKM values) of FSCN1 and FSCN1-regulated DEGs in cervical cancer samples were analyzed. Kaplan–Meier (KM) analysis was used to estimate the survival rate, and the survival curve was plotted with the log-rank test to calculate the P value.

Statistical analysis

Statistical analysis was performed using R (v3.1.3). qPCR data of the two groups were compared by an unpaired two-tailed t test at P values < 0.05. The qPCR data are presented as the mean  ± standard deviation (SD). Pearson correlation analysis was used to assess the correlation between the expression of FSCN1 and FSCN1-regulated DEGs in cervical cancer tissues from TCGA at P values < 0.05.

FSCN1 was overexpressed in cervical cancer tissues and associated with poor prognosis

To verify whether FSCN1 was overexpressed in cervical cancer tissues, we detected the expression of FSCN1 in cervical cancer. We first analyzed the expression of FSCN1 in cervical cancer samples from the TCGA database, which includes the gene expression profile (FPKM) of 306 cervical cancer tumor and 3 normal tissues. The results showed that the expression of FSCN1 was obviously upregulated in cervical cancer tumor tissues compared to normal tissues (Fig. 1A). Moreover, FSCN1 showed relatively high expression in samples of cervical cancer at different stages of development compared to normal tissues (Fig. 1B). Additionally, the expression of FSCN1 was correlated with the prognosis of patients with cervical cancer. Compared to patients with low expression of FSCN1, patients with high expression had shorter overall survival times (Fig. 1C). Thus, FSCN1 might be a potential prognostic factor in cervical cancer.

FSCN1 knockdown broadly affects the gene expression profile of HeLa cells

To explore the targets regulated by FSCN1 in cervical cancer, the expression of FSCN1 was knocked down in HeLa cells by transfection of FSCN1 shRNA (shFSCN1). Compared with empty vector (shCtrl), FSCN1-shRNA effectively reduced the mRNA expression of FSCN1 (Fig. 2A). To further assess the transfection efficiency, we also knocked down FSCN1 with a siRNA having the same target sequence and a nontargeting siRNA. The results showed that this FSCN1 siRNA decreased the expression of FSCN1 at the mRNA and protein levels (Figs. 2B and 2C).

Then, RNA-seq was used to detect the gene expression profiles of FSCN1 knockdown and control HeLa cells. Both shFSCN1 and control HeLa cells had two biological replicates. There were four RNA-seq samples (shFSCN1_1st, shFSCN1_2nd, shCtrl_1st, and shCtrl_2nd). After cleaning the raw sequencing data, at least 73.4 million paired-end reads were obtained for each sample. After mapping of those clean reads to the human genome, more than 57.0 million uniquely mapped reads of each sample were used for gene expression determination and further gene expression analysis (Table S2).

The FPKM value was calculated to represent the expression level of each gene. The results showed that there were 25, 296 genes with FPKM >0 and 13,415 genes with FPKM > 1 in at least one sample (Tables S3 and Table S4). FSCN1 knockdown was further verified by the FPKM values of this gene in HeLa cells (Fig. 2D). Based on the FPKM values of each expressed gene, a correlation matrix was constructed for the four samples that were used for unsupervised hierarchical clustering. As shown, there was a clear separation of the shFSCN1 and control samples, with the two biological replicates clustered together (Fig. 2E). This result demonstrated that FSCN1 knockdown obviously changed the gene expression profile of HeLa cells.

Then, DEGs between the shFSCN1 and control cells were identified by edgeR to further compare the gene expression profile. A total of 3,870 upregulated and 1,173 downregulated DEGs were identified between shFSCN1 and control cells (Fig. 2F). All the DEGs are listed with detailed information including FPKM and fold change values (Table S5). In addition, there was consistent expression of these DEGs in both replicates of the shFSCN1 and control samples (Fig. 2G). These results showed that FSCN1 extensively regulates gene expression in HeLa cells.

Functional analysis of genes upregulated and downregulated by FSCN1 knockdown in HeLa cells

To reveal the potential roles of the DEGs, GO analysis was performed to annotate the 3,870 upregulated DEGs. As an input set of genes, these DEGs were annotated with putative pathways and disease relationships based on mapping to genes with known annotations. The GO analysis was divided into three categories: biological process (GO-P), molecular function (GO-F) and cellular component (GO-C). The results revealed 1,859, 1,859, and 2,124 upregulated genes annotated with 641 GO-P, 237 GO-F, and 186 GO-C terms, respectively (Table S6). The top 10 GO-P terms included DNA-dependent transcription and regulation of transcription (Fig. 3A). The top 10 GO-F terms included DNA binding and nucleic acid binding (Fig. 3B). Moreover, the top 10 GO-C terms included nucleus, chromosome and nucleolus (Fig. 3C). These results indicated that FSCN1 knockdown increased the expression of many transcription factors, which are distributed in the nucleus and can bind DNA.

To further reveal the roles of downregulated DEGs, GO analysis was performed to annotate those DEGs. The results revealed 561, 552, and 707 downregulated genes annotated with 245 GO-P, 74 GO-F, and 77 GO-C terms, respectively (Table S7). The top 10 GO-P terms included angiogenesis and cell adhesion (Fig. 3D). The top 10 GO-F terms included integrin binding and actin binding (Fig. 3E). Moreover, the top 10 GO-C terms included extracellular space and focal adhesion (Fig. 3F). These results indicated that FSCN1 positively regulates the expression of many genes associated with angiogenesis and cell adhesion.

Validation of FSCN1-regulated genes in HeLa cells and cervical cancer tissues

To verify the effect of FSCN1 knockdown on the expression of these DEGs, qPCR was conducted to quantify the changes in the mRNA levels of these genes after FSCN1 knockdown in HeLa cells. Fifteen DEGs were selected for qPCR analysis, namely ten upregulated transcription factors—HBP1 (HMG-box transcription factor 1), HIF1A (hypoxia inducible factor 1, alpha subunit), HLTF (helicase-like transcription factor), LRPPRC (leucine-rich pentatricopeptide repeat containing), TGFBR1 (transforming growth factor, beta receptor 1), DDX17 (DEAD box helicase 17), CNOT1 (CCR4-NOT transcription complex, subunit 10), ZNF664 (zinc finger protein 664), WAC (WW domain containing adaptor with coiled-coil), and PTPN14 (protein tyrosine phosphatase, non-receptor type 14)—and five downregulated genes—KRT18 (keratin 18, type I), BCL2L1 (BCL2-like 1), ANGPTL4, CTGF (connective tissue growth factor), EPHA2 (EPH receptor A2). All selected DEGs had an FPKM value > 1 in at least one sample. The results showed that 13 of the 15 selected DEGs showed a significant increase or decrease after FSCN1 knockdown in HeLa cells, which was in agreement with the RNA-seq analysis (Figs. 4A and 4C). Although there was no significant difference in the qPCR results for ZNF664 and CNOT1, the trend was the same between the qPCR and RNA-seq data (Fig. 4A). These results indicated that FSCN1 negatively regulates the expression of HBP1, HIF1A, HLTF, LRPPRC, TGFBR1, DDX17, CNOT1, ZNF664, WAC, and PTPN14 (Fig. 4A), but positively regulates the expression of KRT18, BCL2L1, ANGPTL4, CTGF, and EPHA2 (Fig. 4C) in HeLa cells.

Then, we analyzed the correlations between the expression levels of FSCN1 and these genes in cervical cancer tissues from TCGA. The results showed that the expression of FSCN1 had a significant negative correlation (P < 0.05) with that of HLTF, HBP1, ZNF664, CTGF, DDX17, and KRT18 in cervical cancer tissues (Fig. 4B). There were significantly positive correlations (P < 0.05) between FSCN1 expression and the expression of PTPN14, TGFBR1, EPHA2, HIF1A and ANGPTL4 (Fig. 4D).

Thus, FSCN1 negatively regulates four transcription factors (HLTF, HBP1, ZNF664, DDX17) and positively regulates ANGPTL4 (an angiogenic factor) and EPHA2 in both HeLa cells and cervical cancer tissues (Fig. 4). Then, we further detected the effect of FSCN1 knockdown on the protein expression level of one selected target in HeLa cells. Western blotting showed that FSCN1 knockdown decreased the protein level of ANGPTL4 in HeLa cells (Fig. 4E). Altogether, these results indicated that FSCN1 is not only an actin binding protein, but also a transcriptional regulator and an angiogenic factor in cervical cancer.