Upregulation of 15 Antisense Long Non-Coding RNAs in Osteosarcoma

The human genome encodes thousands of natural antisense long noncoding RNAs (lncRNAs); they play the essential role in regulation of gene expression at multiple levels, including replication, transcription and translation. Dysregulation of antisense lncRNAs plays indispensable roles in numerous biological progress, such as tumour progression, metastasis and resistance to therapeutic agents. To date, there have been several studies analysing antisense lncRNAs expression profiles in cancer, but not enough to highlight the complexity of the disease. In this study, we investigated the expression patterns of antisense lncRNAs from osteosarcoma and healthy bone samples (24 tumour-16 bone samples) using RNA sequencing. We identified 15 antisense lncRNAs (RUSC1-AS1, TBX2-AS1, PTOV1-AS1, UBE2D3-AS1, ERCC8-AS1, ZMIZ1-AS1, RNF144A-AS1, RDH10-AS1, TRG-AS1, GSN-AS1, HMGA2-AS1, ZNF528-AS1, OTUD6B-AS1, COX10-AS1 and SLC16A1-AS1) that were upregulated in tumour samples compared to bone sample controls. Further, we performed real-time polymerase chain reaction (RT-qPCR) to validate the expressions of the antisense lncRNAs in 8 different osteosarcoma cell lines (SaOS-2, G-292, HOS, U2-OS, 143B, SJSA-1, MG-63, and MNNG/HOS) compared to hFOB (human osteoblast cell line). These differentially expressed IncRNAs can be considered biomarkers and potential therapeutic targets for osteosarcoma.


Introduction
Osteosarcoma (OS), also known as osteogenic sarcoma, is the most common primary malignant solid tumour of bone [1]. The peak incidence is in children and adolescents with a smaller second peak in incidence after the age of 65 years associated with Paget's disease of bone [2]. OS commonly develops in the extremities of long bones such as the distal femur, proximal tibia, proximal humerus, and proximal femur [3]. It is an aggressive-invasion sarcoma type that frequently metastasizes to the lung and other bones in the body [4]. OS usually presents with pain, tenderness and swelling around the affected bone, and diagnosis is achieved by a combination of imaging and histology with the characteristic appearance of malignant cells forming osteoid [5]. Cytotoxic chemotherapy was introduced by Rosen in the 1970s and improved the prognosis from 20% to a 70% five-year survival rate with no further significant improvements in outcome since then [6]. Current treatments of OS include neo-adjuvant chemotherapy with drugs such as doxorubicin, methotrexate, and cisplatin with the aim of reducing tumour size as well as eradicating micro-metastases. Ablative surgery is then followed by further chemotherapy determined by the cell death rate observed in the surgical specimens [7,8]. Current OS therapeutic agents are limited to cytotoxic drugs interfering with transcription and DNA replication [6]. This is a reflection of our knowledge of the pathways involved in OS initiation and progression, which are insufficient to understand the underlying molecular mechanisms of the disease.
The sense strand of DNA provides the template for production of messenger RNA (mRNA) to be translated into proteins [9], but the Human Genome Project highlighted that only 1.5% of the human genome contains protein-coding genes. In addition, the Encyclopedia of DNA elements (ENCODE) and the Functional Annotation of the Mammalian Genome (FANTOM) have suggested that the majority of the genome is transcribed and produces a various amount of non-coding RNA species (ncRNAs) [10,11]. The ncRNAs, fittingly, are RNA molecules that do not encode proteins. They are generally classified based on their length, with an artificial cut off of 200 base pair (bp), small ncRNAs (sncRNAs) less than 200 bps, whereas long ncRNAs (lncRNAs) greater than 200 bps [12]. The lncRNAs frequently regulate epigenetic silencing through chromatin remodeling [13]. They also play critical role in regulating splicing, recruiting transcription factors, and controlling mRNA stability [14]. The natural antisense RNAs belong to the lncRNAs family and are transcribed from the opposite strand of a protein-coding gene [12]. They can stimulate, reduce or completely silence gene expression of the sense transcripts at multiple levels, play a functional role in physiological, and pathological processes, and may eventually lead to diseases [10][11][12][13][14][15].
Nuclear RNA duplex formation may occur locally immediately upon transcription, consequently resulting in inhibition of sense RNA processing ( Figure 1A) [16]. The processing of RNA includes capping, polyadenylation, nuclear localization and mRNA transport, all these events may be affected by nuclear sense-antisense RNA duplex formation. Natural antisense transcripts (NATs) can cover donor (5 ) and acceptor (3 ) splice sites in the sense precursor mRNA transcript to modify alternative splicing patterns and develop mature RNA in different isoforms ( Figure 1B) [16,17]. Another possible consequence of nuclear RNA duplex formation is RNA editing through the adenosine deaminases that act on the RNA (ADAR) enzyme responsible for binding to double stranded RNA and converting adenosine (A) to inosine (I) by hydrolytic deamination (an adenosine loses an amine group) ( Figure 1C) [18,19].
NATs can be subdivided into two categories: cis-NATs and trans-NATs. cis-NATs are antisense RNAs transcribed from the same genomic locus. Consequently, in the section of the overlap, sense and antisense transcripts share complete complementarities such as head-to-head overlap in cis (Figure 1Da), embedded overlap in cis (Figure 1Db), and tailto-tail overlap in cis (Figure 1Dc) [20]. On the other hand, trans-NATs are antisense RNAs transcribed from a different genomic region of their paired sense transcript, displaying partially complementarities such as overlap in trans (Figure 1Dd) [20][21][22].
Considering the growing evidence of the antisense lncRNAs in cellular process and their involvement in various disease types including cancer we investigated their expression pattern in OS using RNA sequencing (RNA-seq).

Sample Description
The study has been approved by the Ethics Review Committee of the University of Western Australia and Sir Charles Gairdner Hospital (2019/RA/4/20/5211). The patient informed consent forms were signed and dated by the participants and patient representatives before the limb sparing or amputation surgery.
Twenty-one Australian patients underwent the surgery to remove malignant tissue and the diagnosis of OS was confirmed by a specialist sarcoma pathologist. Cancerous and normal bone formalin-fixed paraffin-embedded (FFPE) tissue samples were collected from PathWest (QEII Medical, Nedlands, WA, Australia).

Total RNA Extraction and Sequencing
Total RNA was extracted from recently cut 5 sections of ≤ 20 µm thick FFPE samples using the FFPE RNA purification kit (Norgen Biotek, Thorold, ON, Canada), following to the manufacturer's protocol. Extracted RNA was completely dissolved in RNase-free water and stored at -80 • C. The quality of total RNA was measured using Agilent 2100 Bioanalyzer and the RNA 6000 Nano Kit (Agilent Technologies Inc., Santa Clara, CA, USA).
The RNA samples were processed using the Takara SMARTer V2 Total RNA Mammalian Pico Input protocol using 2 ng of Total RNA input as per manufacturer's instructions (Takara Bio Inc., Mountain View, CA, USA). The sequencing was completed using an Illumina NovaSeq 6000 and an S4-300 cycle lane (150PE) with v1.5 sequencing chemistry [23]. The FastQC (version 0.11.9) was used to determine the quality score distribution of the sequencing reads. The low-quality reads, Phred score ≤ 20, were trimmed out using Trimmomatic (version 0.39) [24]. Trimmed reads were mapped to the human genome hg38 (GRCH38) reference using STAR (version 2.7.7a) [25].

Data and Statistical Analyses
Differential gene expression and statistical data analyses were performed through DESeq2 package for R [26]. DESeq2 is a Bioconductor package specifically designed to detect differential expressed genes between individual samples. Differential gene expression levels between tumour and normal samples were obtained by the package through investigating the logarithmic-2-fold changes of the genes (logFC, the cut-off value of 0.5). The DESeq2 package also provides the Benjamini-Hochberg method to adjust p-value (padj). The significant level was set at padj < 0.05.

FANTOM-CAT Analysis
FANTOM-CAGE-Associated Transcriptome (FANTOM-CAT) is an efficient software to analyze lncRNAs structure and function (http://fantom.gsc.riken.jp/cat/ (accessed on 7 April 2021)). Zenbu tool was used to determine the type of sense-antisense overlap of the genes through FANTOM-CAT [27].

Transcript
Forward Primer Reverse Primer Quantitative PCR was conducted using a Viia 7 Real-Time PCR machine (Applied Biosystems, Foster City, CA, USA). The thermal cycler protocol used for the RT-qPCR as following; 50 • C for 2 min, 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, and 60 • C for 60 s. GAPDH was used as a housekeeping gene to normalize gene expressions. Relative gene expressions were calculated using the 2 −∆∆CT method. Data analysis was performed using GraphPad Prism software, version 8 (GraphPad software, San Diego, CA, USA). The data are presented as the mean ± standard error (SE) of values from 3 independent experiments. Statistical significance was determined using one-way ANOVA, with p < 0.05 considered statistically significant.

Participants' Characteristics
The samples were collected after surgical removal of the affected bone from 21 Australian OS patients. Only 16 samples contained paired tumour and healthy bone biopsies. Total samples number were noted as 24 tumour and 16 healthy bone biopsies ( Table 2).

Differential Gene Expression Analysis between Tumour and Normal Samples
Total RNA was collected from 24 tumour and 16 normal FFPE samples. Statistical data and differential gene expression analyses were performed through DESeq2 package for R. The 3D principal component analysis (PCA) plot has been provided insights into the similarities between tumour and normal samples and indicated the quality of the gene expression data (Figure 2).  Table 3). The data were visualized using heatmap ( Figure 3) and circos plot ( Figure 4). Further, the type of overlap of the antisense lncRNA transcripts was determined using FANTOM-CAT analysis and listed in Table 4.

Discussion
NATs are a growing focus of cancer genomics studies. They have been dysregulated in various cancer types, are implicated in several malignant phenotypes [28][29][30], and are emerging as pre-dominant players in carcinogenesis through their involvement in gene expression regulation, epigenetic modification, evasion of growth suppressors and reprogramming energy metabolism [28].
Dysregulation of RUSC1-AS1 (also known as C1orf104) has been associated with several cancer types. According to some studies RUSC1-AS1 is highly expressed in laryngeal squamous cell carcinoma, cervical cancer, and breast cancer cells [31][32][33]. Another study supports that the transcript promotes cell proliferation in hepatocellular carcinoma through modulating NOTCH signaling [34]. Our RNA-seq data also highlighted that the transcript was upregulated in OS samples compared to normal. According to our RT-qPCR result RUSC1-AS1 was significantly upregulated in U2-OS cell line by more than a 2-fold ( Figure 5B). Whereas, the transcript was dramatically downregulated in SaOS-2, HOS and 143B cell lines. The cell lines have different characteristics, morphology, and metastatic properties, consequently the differences can affect transcript expression. Further, passage number affects a cell line's characteristics over time such as cell lines with high passage numbers can experience alterations in morphology, response to stimuli, cell growth rates, gene and protein expression and transfection efficiency, compared to lower passage cells [35][36][37][38][39].
A study highlighted a potential regulatory connection of TBX2-AS1 and TBX2. The same study also suggested that TBX2-AS1 tightly co-expressed with TBX2 suggesting cisregulation and their association with neuroblastoma [40]. Our RT-qPCR also has validated alongside with RNA-seq data that TBX2-AS1 was upregulated in G-292, SJSA-1, 143B, U2-OS, and MNNG/HOS by more than 15-fold compared to hFOB ( Figure 5C).
The knockdown of heterogeneous nuclear ribonucleoprotein K (hnRNPK) eventually reduced PTOV1-AS1 expression in HeLa cervical carcinoma cells. The study also investigated reduced expression of hnRNPK or PTOV-AS-1 suppressed heme oxygenase-1 (HO-1) expression by increasing the enrichment of HO-1 mRNA in miR-1207-5p-mediated miRISC. The knockdown or decreased expression of either hnRNPK or PTOV-AS-1 resulted in inhibition of the proliferation and clonogenic ability of HeLa cells [41]. We also observed upregulation of PTOV-AS-1 in tumour samples compared to normal bone tissue. Interestingly, we also observed hnRNPK upregulation in tumour compared to normal samples ( Figure S1).
The single-nucleotide polymorphisms (SNPs) of ZMIZ1-AS1, located at 10q22.3, has been associated with colorectal cancer and patients' survival among Korean population [42,43]. The transcript also interferes with ZMIZ1 gene regulating several tumour suppressors such as SMAD4 and p53 [44,45]. According to some studies RNF144A-AS1, also known as GRASLND, is highly expressed in bladder cancer, and overexpression of the transcript is correlated with poor prognosis. The same study also observed that the knockdown of RNF144A-AS1 eventually inhibited cell proliferation, migration, and invasion in J82 and 5637 cell lines, in addition xenograft growth of cells was reduced compared to negative control in nude mice [46,47]. Another study highlighted RNF144A-AS1 as an important regulator of mesenchymal stem cell chondrogenesis [48].
Over expression of TRG-AS1 has been observed in several cancer types including squamous cell carcinoma of the tongue, hepatocellular carcinoma, and glioblastoma. Furthermore, TRG-AS1 has been associated with poor prognosis [49][50][51].
Interestingly, GSN-AS1 was downregulated in breast cancer patients [52]. The transcript was upregulated in OS tumour compared to normal samples in our study.
Interestingly, overexpression of OTUD6B-AS1 inhibits cell proliferation, migration, invasion, and promotes cell apoptosis in colorectal cancer by sponging miR-21-5p and regulating PNRC2 [61]. Another study has also suggested the similar findings, overexpression of the transcript inhibits cell proliferation, migration, and invasion by downregulation of microRNA-3171 [62]. Gang Wang et al. suggested that OTUD6B-AS1 expression was downregulated in renal cell carcinoma via the Wnt/β-catenin signaling pathway and low expression of the transcript has been correlated with shorter overall survival than patients with high OTUD6B-AS1 expression [63]. The overexpression of the transcript also reduced cell migration and invasion in thyroid carcinoma cells [64]. Whereas high expression of OTUD6B-AS1 indicates poor prognosis in ovarian cancer [65]. The transcript is also upregulated in OS samples compared to normal bone tissues. The RT-qPCR result has also validated that the transcript was upregulated in SaOS-2, and HOS cell lines compared with hFOB ( Figure 5A).
Chaoyang Zhou et al. found that COX10-AS1 acts as a competing endogenous RNA to positively regulate ACTG1 expression by sponging miR-361-5p and promotes glioblastoma cell proliferation and inhibits apoptosis [66]. The transcript is also upregulated in SaOS-2, and G-292 cell lines compared to hFOB by more than 2-fold in our RT-qPCR result ( Figure 5D).
The upregulation of SLC16A1-AS1 was observed in hepatocellular carcinoma, and glioblastoma [67,68]. Hong Yue Liu et al. suggested that the transcript was dramatically downregulated in non-small cell lung cancer and over expression of SLC16A1-AS1 inhibits the cell viability and proliferation of lung cancer cell [69]. This is the first time of ERCC8-AS1, UBE2D3-AS1, RDH10-AS1, and ZNF528-AS1 expressions in cancer have been reported. According to GeneCards (genecards.org (accessed on 10 April 2021)) ERCC8-AS1 has been associated with Cockayne Syndrome A and Cockayne Syndrome which one of main clinical features is cachectic dwarfism. It is a photosensitive, DNA repair disorder which has been associated with progeria that is caused by a defect in the transcription-coupled repair sub-pathway of nucleotide excision repair [70,71]. Another study has suggested that individuals with Cockayne Syndrome have mutations in ERCC8 and ERRC6, resulting in defective transcription-coupled nucleotide excision repair. In addition, ERCC1 or ERCC4 mutation also have been reported in Cockayne Syndrome [72]. ERCC family also widely involved with Fanconi anemia which leads to bone marrow failure, several moderate skeletal abnormalities and a predisposition to leukemia and solid tumours [73,74]. Interestingly, downregulation of ERCC4 has displayed worse survival outcome in OS patients [75]. Our RT-qPCR result also highlighted that ERCC8-AS1 was upregulated in G-292, SJSA-1, HOS, 143B, US-OS and MNNG/HOS cell lines compared to hFOB ( Figure 5E). UBE2D3-AS1 is also upregulated in SaOS-2, and G-292 cell lines compared to hFOB ( Figure 5G).
The limitation of this study is that some patients received chemotherapy may lead to gene expression alterations of the results. Another limitation of the study is the cell lines have had different passage numbers which may affect RT-qPCR results.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/genes12081132/s1, Figure S1. Table of   Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.