A Newly Identified lncBCAS1-4_1 Associated With Vitamin D Signaling and EMT in Ovarian Cancer Cells

Long noncoding RNAs (lncRNAs) were identified rapidly due to their important role in many biological processes and human diseases including cancer. 1α,25-dihydroxyvitamin D3 [1α,25(OH)2D3] and its analogues are widely applied as preventative and therapeutic anticancer agents. However, the expression profile of lncRNAs regulated by 1α,25(OH)2D3 in ovarian cancer remains to be clarified. In the present study, we found 606 lncRNAs and 102 mRNAs that showed differential expression (DE) based on microarray data. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis indicated that the DE genes were mainly enriched in TGF-β, MAPK, Ras, PI3K-Akt, and Hippo signaling pathways, as well as the vitamin D-related pathway. We further assessed the potential lncRNAs that linked vitamin D signaling with EMT, and lncBCAS1-4_1 was identified in the first time. Moreover, we found that the most upregulated lncBCAS1-4_1 showed 75% same transcripts with CYP24A1 (metabolic enzyme of 1α,25(OH)2D3). Finally, the lncBCAS1-4_1 gain-of-function cell model was established, which demonstrated that the knockdown of lncBCAS1-4_1 inhibited the proliferation and migration of ovarian cancer cells. Furthermore, lncBCAS1-4_1 could resist the antitumor effect of 1α,25(OH)2D3, which was associated with upregulated ZEB1. These data provide new evidences that lncRNAs served as a target for the antitumor effect of 1α,25(OH)2D3.


INTRODUCTION
Ovarian cancer is the leading cause of death caused by gynecologic malignancies (1). Despite the significant medical advances during the past decades, the 5-year survival rate of ovarian cancer is lower than 50% (2). Long noncoding RNAs (lncRNAs), with transcripts more than 200 nucleotides in length, were suggested to play fundamental roles in the development of tumor (3), due to the fact that lncRNAs may exhibit tumor suppressive or promoting functions through the regulation of transcription, translation, protein modification, and the formation of RNA-protein or protein-protein complexes (4). Therefore, lncRNAs are possible candidates of cancer biomarkers and/or therapeutic targets (5). Recently, 13 published papers investigated the expression of lncRNA in normal ovaries, ovarian cysts, and benign and malignant ovarian cancer (6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18), suggesting the important role of lncRNAs in ovarian cancer development and chemotherapeutic survival outcomes of patients. Thus, it is important to explore the potential of lncRNAs as a therapeutic target of ovarian cancer.
In the present study, the lncRNA and mRNA networks were constructed using microarray data, which were used to explore the profile of lncRNA in 1a,25(OH) 2 D 3 -treated human ovarian cancer SKOV3 cells comprehensively. Moreover, the potential lncRNAs that linked vitamin D signaling with EMT were analyzed, and lncBCAS1-4_1 was identified. Besides, the effect of lncBCAS1-4_1 on the proliferation and migration in 1a,25 (OH) 2 D 3 -treated ovarian cancer cells was investigated.

Microarray Expression Profiling
SKOV3 cells were treated with 1a,25(OH) 2 D 3 (100 nmol/L) or vehicle (the same concentration of ethanol) for 72 h. Total RNA was extracted with a TRIzol reagent (Thermo Fisher Scientific, Scotts Valley, CA, USA) and quantified using NanoDrop ™ ND-2000 (Thermo Fisher Scientific). After RNA integrity was assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA), sample labeling, microarray hybridization, and washing were performed based on the manufacturer's standard protocols (OE Biotech Company, Shanghai, Design ID: 076500). Briefly, total RNA was transcribed to double-stranded cDNA and then synthesized into cRNA, which was labeled with Cyanine-3-CTP. The labeled cRNAs were hybridized onto the microarray. After washing, the arrays were scanned by the Agilent Scanner G2505C (Agilent Technologies).

Differentially Expressed Gene Analysis
Limma (Version 3.8) package in R software was used to identify the differently expressed mRNAs (DE-mRNAs) and -lncRNAs (DE-lncRNAs) with a threshold of |log2 (fold change [FC])| > 2.0 and a false discovery rate [FDR (adjusted p-value)] < 0.05. The heatmap and volcano were constructed by the gplots package in R software.

Functional Enrichment Analysis
To reveal the functions of DE genes, the Enricher database was used to conduct Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses (27). The GO terms comprised of the following three divisions: biological process (BP), cellular component (CC), and molecular function (MF). A significance level of p < 0.05 was set as the cutoff criterion, and the plots were constructed by the gplots package in R software.
A PPI network of DE mRNAs was constructed using STRING 11.0 (http://string-db.org), with a combined score > 0.9 as the cutoff value. Significant modules in the PPI network were identified using MCODE 1.5.1 (a Cytoscape software plug-in).

Construction of the lncRNA-mRNA Co-Expression Modules
The lncRNAs and mRNAs co-expression modules were further selected using Pearson correlation analysis. The lncRNA-mRNA pairs with a correlation coefficient > 0.9 and p < 0.05 were used for bidirectional clustering.

Quantitative Real-Time PCR
Reverse transcription reactions consisted of 1 mg RNA and 2 mL of 5xPrimerScript RT Master Mix (TaKaRa, Japan) with a total volume of 10 mL. The primer sequences of RNA are shown in Table 1. Reactions were performed in a C100 PCR System (Bio-Rad, Hercules, CA, USA) for 15 min at 37°C. GAPDH was used as the internal control. The qPCR was performed using the SYBR Green (Roche, Basel, Switzerland) dye detection method on the ABI 7500 PCR instrument (Applied Biosystem, Foster City, CA, USA) under default conditions: 95°C for 10 s, 40 cycles of 95°C for 5 s, and 55°C for 30 s. The relative gene expression levels were analyzed by the 2 -DDCt method, where DCt = Ct(target) -Ct(GAPDH).

Construction of lncBCAS1-4_1 Loss/Gain Cell Model
Overexpression adenoviruses (OE) as well as control adenoviruses (empty vector, EV) of lncBCAS1-4_1 were purchased from GeneChem Corporation (Shanghai, China). The knockdown lncBCAS1-4_1 was produced by siRNA interference. Scramble control and silncBCAS1-4_1 were purchased from RiboBio Co., Ltd (Guangzhou, China) and transfected using riboFECT ™ CP (Guangzhou, China) according to the manufacturer's instructions. All oligonucleotide sequences are listed in Table 2.

Cell Proliferation Assay
Cell colony formation and the CCK-8 counting were used to assay cell proliferation, respectively. Briefly, 3×10 4 OVCAR8 or 2.5×10 4 SKOV3 cells were seeded onto 60-mm culture plates and transfected by adenoviruses or siRNAs for 72 h. After the cells were treated with 1a,25(OH) 2 D 3 (100 nmol/L) or ethanol for 48 h, they were fixed with 75% alcohol and stained with 0.3% methyl violet for 20 min at room temperature. Then the colonies were dissolved by glacial acetic acid, and the absorbance value (AU) was detected at 585 nm with a microplate reader (Filter Max F5, Molecular Devices, CA, USA). The cell proliferation ratio was calculated as (AU treatment group -AU blank group )/(AU control group -AU blank group ). All experiments were performed in triplicate.
According to the manufacturer's instructions, approximately, 3×10 3 OVCAR8 or 2.5×10 3 SKOV3 cells per well were plated in triplicate into 96-well plates and treated with 1a,25(OH) 2 D 3 for 48 h. The control group was treated with ethanol. At each of the desired time points, 10 mL of the CCK-8 solution was added and incubated for 1 h at 37°C, followed by measurement of absorbance at 450 and 630 nm with a microplate reader for quantifying the relative cell density. Cell viability was calculated as: (AU 450-treatment group -AU 630-treatment group )/ (AU 450-control group -AU 630-blank group ). All experiments were performed in triplicate.

Cell Migration Assay
The cell migration was assessed using a wound healing assay. Cells were plated into a 6-well plate with FBS-free media for 12 h. Afterwards, cells cultured in the bottom of the well were scratched using a pipette tip to create a wound area. After 24 and 48 h, wounds (three images each well) were imaged under a microscope (40, CKX41F, Olympus, Tokyo, Japan) to detect the width of the gaps. Wound healing assay data are displayed as the migration index (%), which is calculated by the formula [(initial width) -(final width)]/(initial width). Values were normalized by the control group. Data points in the figure represent three independent experiments.

Statistical Analysis
All microarray statistical data were analyzed in the R environment (R version: 3.6.3). Wilcoxon/Mann-Whitney test was used to analyze continuous variables, and Fisher's exact test or chi-square test was used to analyze the categorical data. Experimental data were performed using GraphPad Prism 8 (GraphPad Software Inc., La Jolla, CA, USA). Quantitative data were presented as the mean ± standard deviation (SD). Statistical data were analyzed using an unpaired Student's t-test. For all statistical analyses, a p-value less than 0.05 was regarded as statistically significant.
To further validate the findings of the microarray analysis results, five dysregulated lncRNAs were confirmed using quantitative RT-PCR. lnc-BCAS1-4_1 and lnc-RWDD4-5_1 were selected as target lncRNAs with the most upregulated/ downregulated expression. Lnc-ZNF599-3_6 was selected for its potential function of trans-regulating, and the other two (lnc-MBOAT1-4_2 and lnc-KRT7-2_2) were randomly selected. Consistently, their expressions from quantitative RT-PCR results were similar with those of the microarray analysis ( Figure 1E). Similarly, the transcriptional levels of lncRNABCAS1-4_1 and CYP24A1 were indeed dramatically increased after 1a,25(OH) Vitamin D-Regulated lncRNA-mRNA Network in Ovarian Cancer Cells To explore the potential responsible mechanism of cancer cells for 1a,25(OH) 2 D 3 , KEGG pathway analysis was performed on the DE genes. The results indicated that DE mRNAs were mainly enriched in TGF-b, regulating pluripotency of stem cells, and Hippo signaling pathways (Figure 2A). The hub genes with a degree connectivity in PPI network were enriched in insulin-like growth factor 1 (IGF1), which is known to induce cell proliferation (28), TGF-b2 (29), insulin-like growth factor-binding protein 3 (IGFBP3) (28), and COL1A1 (30), which are closely associated with the vitamin D endocrine system ( Figure 2B). Then, we identified a top 5 lncRNAs-mRNAs networks including 5 lncRNAs and 140 mapped mRNAs ( Figure 2C and Supplementary Data Table S1). GO enrichment analysis and subpathway analysis showed that "phagocytosis", "cytoplastic side of plasma membrane", and "growth factor activity" were significantly related to this module ( Figure 2D). KEGG analysis for 140 mRNA from top 5 lncRNA-mRNAs networks revealed that cancer-related pathways were enriched in this network, e.g., Ras, MAPK, TGF-b, Rap1, and PI3K-Akt signaling pathway ( Figure 2E).
Construction of the lncBCAS-1_4-1 as a Core of EMT Signal Pathway in 1a,25(OH) 2

D 3 Treated Ovarian Cancer Cells
Next, we selected the most dysregulated lncRNA, lncBCAS1-4_1, to construct the lncRNA-mRNA network, and 83 mapped mRNAs were involved to explore the function of this module ( Figure 3A and Supplementary Data Table S2). GO analysis showed that "epithelial cell proliferation", "collagen-containing extracellular matrix", and "growth factor activity" were highly enriched in this network ( Figure 3B). KEGG analysis also revealed that these genes mainly enriched in TGF-b, regulating pluripotency stem cells, MAPK, Ras, and Hippo signaling pathways ( Figure 3C). Because the TGF-b signaling pathway repeatedly occurred, the EMT-related genes were applied to identify the significant pathway associated with 1a,25(OH) 2 D 3 ; as shown in Figure 3D. The EMT pathway was significantly activated in this network ( Figure 3D).

The Role lncBCAS-1_4-1 on Proliferation and Migration of Ovarian Cancer Cells
To validate the function of lncBCAS-1_4-1, SKOV3 cells were used to build up lncBCAS1-4_1 gain-of-function cell models ( Figure 4A), while OVCAR8 cells were used to build up lncBCAS1-4_1 loss-of-function cell models ( Figure 4B). The results of CCK8 ( Figure 4C) and platting efficiency ( Figure 4D) assay showed that overexpressed lncBCAS1-4_1 promoted proliferation, while knockdown of lncBCAS1-4_1 inhibited proliferation. Similarly, we found that the gain of lncBCAS1-4_1 increased migration, and the loss of lncBCAS1-4_1 decreased cell migration ( Figure 4E). We then detected the expression of mRNAs associated with the EMT signaling pathway. The result demonstrated that overexpression of lncBCAS1-4_1 significantly upregulated the EMT mesenchymal marker including N-cadherin and Vimentin, as well as the EMT-related transcriptional factor (ZEB1) ( Figure 4F).

The Inhibition of 1a,25(OH) 2 D 3 on Proliferation and Migration of Ovarian Cancer Cells Is Disrupted by lncBCAS-1_4-1
To ascertain the impact of lncBCAS1-4_1 on the antitumor action of vitamin D, ovarian cancer cells treated with 1a,25(OH) 2 D 3 were interfered or overexpressed by the siRNA or adenovirus vector of lncBCAS1-4_1, respectively. As expected, the knockdown of lncBCAS1-4_1 significantly enhanced the 1a,25 (OH) 2 D 3 mediated antitumor effect, while overexpressed lncBCAS1-4_1 resisted the antitumor effect of 1a,25(OH) 2 D 3 in vitro (Figures 5A-C). The results from Figure 5D showed that the expressions of Vimentin, ZEB1, and Twist1 were significantly reduced by 1a,25(OH) 2 D 3 as compared to mock-vehicle negative control. However, the reduced ZEB1 levels in overexpressed lncBCAS1-4_1 SKOV3 cells were increased after treatment with 1a,25(OH) 2 D 3 . Taken together, these data indicated that the overexpression of lncBCAS1-4_1 significantly resisted the antitumor effect of 1a,25(OH) 2 D 3 , which was associated with upregulating ZEB1.
The co-expression/regulatory networks of lncRNAs-mRNAs indicated that "TGF-b signaling pathway", "epithelial cell proliferation", and "Hippo signaling pathway" were significantly involved in 1a,25(OH) 2 D 3 treated cancer cells. Up to date, there are lots of reports about how coding genes or non-coding genes to regulate EMT progress or EMT associated genes and transcription factors (48)(49)(50)(51). Moreover, 1a,25(OH) 2 D 3 was reported to have the effect on inhibiting the progression of EMT (31). Our results also showed that EMT signaling pathway was significantly activated in 1a,25(OH) 2 D 3 treated ovarian cancer cells. It is plausible that these lncRNAs could mediate the EMT process by vitamin D signaling pathway, which supports our hypothesis that 1a,25(OH) 2 D 3 has inhibitory effects on ovarian cancer cells by regulating lncRNA expression patterns.
In the present study, the most upregulated lncRNA was lncBCAS1-4_1, which has the closest relationship with the mRNA transcript of CYP24A1, because their 75% transcripts are the same. CYP24A1 is the gene coding the metabolic enzyme of 1a,25(OH) 2 D 3 , resulting in the loss of physiological activity by 1a,25(OH) 2 D 3 (34). In vitro and in vivo studies also showed that CYP24A1 has been deemed as a candidate oncogene in many cancers, such as ovarian cancer (35), colorectal cancer (36,37), prostate cancer (38), lung cancer (39), breast cancer (40), thyroid cancer (41), and so on. Moreover, a recent study showed that the upregulation of CYP24A1 and PFDN4 as well as nearby lncRNAs may be used as the potential diagnostic biomarker in colorectal cancer (52). Interestingly, it has been reported that mice with CYP24A1 knockout exhibited a fourfold reduction in thyroid tumor growth compared with wild-type CYP24A1 mice. They found that this phenotype was associated with the repression of the MAPK, PI3K/Akt, and TGFb signaling pathways, and a loss of EMT in CYP24A1 knockout cells was also associated with the downregulation of genes involved in EMT, tumor invasion, and metastasis (53). Furthermore, functional analysis revealed that the TGF-b pathway was associated with lncBCAS1-4_1. Based on the 75% similarity with CYP24A1 and the relationship between CYP24A1 and EMT, as well as the key role of TGF-b in the EMT process, we focused on the link of lncBCAS1-4_1 and EMT. For the lncBCAS1-4_1 loss/gain cell model, the oncogenic role of lncBCAS1-4_1 was validated in vitro, and the overexpression of lncBCAS1-4_1 significantly resisted the antitumor effect of 1a,25(OH) 2 D 3 , which was associated with upregulated ZEB1. Thus, it was worthy to reveal the molecular mechanism of EMTrelated lncRNAs in cancer and to demonstrate that lncBCAS1-4_1 can be a potential therapeutic target for patients.
Additionally, we also found that the most downregulated lncRNA (lnc-RWDD4-5_1) and IGFBP3 mRNAs were negatively correlated. After treatment with 1a,25(OH) 2 D 3 , the expression of lnc-RWDD4-5_1 was dramatically decreased, while that of IGFBP3 was increased (2.2-fold change). The most hub gene in the PPI network was IGF1, which can bind to IGFBP3. IGF1 and its binding proteins can promote cellular proliferation and inhibit apoptosis. In vitro studies showed that IGF1 increased ovarian cell growth and invasive potential (42). It is well documented that high IGF1 levels are significantly associated with early-stage cancer, nonserous histology, and optimal cytoreduction in epithelial ovarian cancers (43)(44)(45). Considerably, it is noteworthy that the most downregulated lncRNA (lnc-RWDD4-5_1) has a potential relationship with the hub gene (IGF1). However, the potential molecular mechanisms needed to be further verified.
There are also a couple of limitations to this study. Firstly, although the SKOV3 cell line is a useful model of ovarian cancer cells, it could not be used to predict the performance of 1a,25 (OH) 2 D 3 in actual tumors. And the action of 1a,25(OH) 2 D 3 refers to different sets of genes in different cell lines (54)(55)(56). Secondly, the expressions of lncRNAs and mRNAs were analyzed by ovarian cancer cells, and further testing of these in the tumor tissues of patients is needed. Thirdly, the relationships among noncoding RNAs, mRNAs, and proteins need to be further investigated using bioinformatic prediction to understand their full function. Nevertheless, this is the first study of lncRNA expression patterns regulated by 1a,25(OH)

CONCLUSIONS
In summary, we identified the 606 DE lncRNAs and 102 DE mRNAs in 1a,25(OH) 2 D 3 -treated ovarian cancer cells, which were mainly enriched in the cancer-related and vitamin Drelated pathway. Moreover, by the lncBCAS1-4_1-mRNA core network, the EMT signal was identified, indicating the linkage of lncBCAS1-4_1 between EMT and vitamin D signaling. Furthermore, we established the lncBCAS1-4_1 loss/gain cell model and found that lncBCAS1-4_1 could abolish the antitumor effect of 1a,25(OH) 2 D 3 , which was associated with upregulating ZEB1. These data provide new evidence that lncRNAs can serve as targets for the antitumor effect of 1a,25 (OH) 2 D 3 .

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GEO Database and accession number GSE17363 https://www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc=GSE173633.