AP003352.1/miR-141-3p axis enhances the proliferation of osteosarcoma by LPAR3

Data collection

The mRNA and lncRNA sequencing data of 88 patients with osteosarcoma and their corresponding clinical information were downloaded from the TCGA database as cluster analysis data and COX analysis data. The miRNA sequencing data of 65 patients with osteosarcoma and their corresponding clinical information were downloaded from the GSE39040 GEO database as cluster analysis data and COX analysis data. The TARGET-OS data downloaded from the TCGA database was used as the training cohort, while the TCGA-sarc data downloaded from TCGA database was used as the validation cohort. The official websites of TCGA database and GEO database are as follows:

TCGA: https://www.cancer.gov/ccg/research/genome-sequencing/tcga.

GEO: https://www.ncbi.nlm.nih.gov/geo/.

Univariate COX analysis

Univariate COX analysis was applied to evaluate the overall survival rate and status of osteosarcoma patients in the TCGA database, which was then applied to assess the prognostic value of the mRNAs, miRNAs, and lncRNAs. Univariate COX analysis was used to assess the prognostic value of risk scores, deriving hazard ratios (HR) and 95% confidence intervals (CI) for each variable.

Differential expression analysis

Before differentially expressed gene (DEG) identification, the data were normalized using the limma package (3.52.4). Meanwhile, the data type used for the DEG identification was the new TPM value after log2 transformation. According to the new standard of log FC ≥ 1 and p < 0.05 for screening the condition, the DEGs between osteosarcoma tumor samples and precancerous samples were analyzed, and the obtained results could be used for unicox analysis and LASSO screening.

Gene enrichment analysis

Gene enrichment analysis is an analytical method used to analyze the expression and enrichment of DEGs. The R package clusterProfiler (4.4.4) was used to analyze the enrichment of intersecting genes on different signaling pathways. P-value < 0.05 indicated that the enrichment genes in corresponding pathways were significant.

Construction of the lncRNA-miRNA-mRNA ceRNA network

The ceRNA network of lncRNA-miRNA-mRNA was created by predicting the correlation between lncRNA, miRNA, and mRNA. miRNA is a key component of the ceRNA network, which can connect with lncRNA and mRNA. The R package multiMiR (1.18.0) was used to analyze the ceRNA network. In this study, the network was visualized by Cytoscape software.

Cell lines and cell culture

HEK293T, U2OS and MG63 cell lines were purchased from ATCC. They were cultured in DMEM with 10% FBS, 100 μg/mL penicillin, and 100 μg/mL streptomycin. Meanwhile, HEK293T, U2OS and MG63 cells were incubated in a humidified incubator with 5% CO₂ at 37 °C.

CCK8 assay

CCK8 assay was performed to test the cell viability. A total treated of 5 × 10³ U2OS cells were inoculated into 96-well culture plates. After 24 h of cell culture, CCK-8 was added into the plates with 10 μL per well 3 h before testing. Afterwards, the absorbance value at was measured 450 nm with an enzyme marker. Each group of experiments was repeated for three times.

Clone formation assay

To clarify the effect of ceRNA axis on clone formation, single cells of U2OS were cultured in six-well plates with 1 × 10³ cells per well. Withing 10 days of incubation, the colonies were clearly visible, and polyformaldehyde was used to fix cells, while crystal violet was used for cell staining.

Western blot

Cells were lysed using RIPA lysis buffer, and the BCA kit was used for cell quantification. The protein was separated by SDS-PAGE and then transferred to the NC membrane. After blocking with 5% skimmed milk, the NC membrane was incubated with primary antibody at 4 °C for 14 h. TBST was used to wash the NC membrane three times, followed by sealing with the secondary antibody for 2 h. Finally, the protein expression was observed under electrophoresis gel imaging. The antibodies used include: LPAR3 (Invitrogen, Waltham, MA, USA; PA5-27074; 1:2,000) and α-Tubulin (CST, Tarzana, CA, USA; 2148S; 1: 1,000).

RT-qPCR assay

RNAiso Plus was used to extract total RNA from the U2OS cells. HiScript III RT SuperMix qPCR kit was used to synthesize the cDNA. The expression of lncRNA relative to mRNA was determined using the ratio of sample transcripts. The primers used are shown below:

LPAR3 forward: 5′-TCGCTTACGTGTTCCTGATG-3′

LPAR3 reverse: 5′-TTCCACAGCAATAACCAGCA-3′

has-miR-141-3p: 5′-TAACACTGTCTGGTAAAGATGG-3′

AP003352.1 forward: 5′-TGCCTCAGCCTCTCAAGTAG-3′

AP003352.1 reverse: 5′-CGTGGCTCACACCTGTAATC-3′

Dual-luciferase reporter assay

PGL4.15 was performed to construct plasmids containing the AP003352.1-wild type and mutant. U2OS cells were inoculated onto the 24-well plate. After cell adhesion, wild type and mutant AP003352.1 were transfected into U2OS with lipofectamine 3000 (Invitrogen, Waltham, MA, USA; L3000001). After transfection, luciferase activity was detected according to the specification of the Dual-Luciferase Reporter Assay System Kit (Promega, Madison, WI, USA; E1910).

Transfection

The inhibitors and mimics of has-miR-141-3p were bought from Biomics (Nantong, China). The ASOs (Antisense Oligonucleotide) of lncRNAs AL358332.1, AC124798.1, LINC02280, AP003352.1, AC073487.1 and AC084855.2 were purchased from RiboBio (Guangzhou, China). The inhibitors and mimics of miR-141-3p and the ASOs of lncRNAs were transfected into U2OS cells using Lipofectamine 3000 (Invitrogen, Waltham, MA, USA; L3000001).

Construction of shRNAs

The shRNAs of LPAR3 were constructed on the PLKO.1 plasmid. PSPAX2, PMD2G and shLPAR3 plasmids were transfected into HEK293T cells. After 48 h, the supernatant was collected and filtered through a 0.45 um filter membrane. Next, U2OS and MG63 cells were infected with virus solution, and 48 h later, fresh medium was replaced and 5 mg/mL puromycin was added to screen stably transfected cell lines. The sequences of shRNAs are shown below:

shLPAR3-1

5′-CCGGAATACATAGGCAATTCCAGCGCTCGAGCGCTGGAATTGCCTATGTATTTTTTTG-3′

shLPAR3-2

5′-CCGGCAGTACATAGAGGATAGTATTCTCGAGAATACTATCCTCTATGTACTGTTTTTG-3′

Pearson’s correlation coefficient

Pearson’s correlation coefficient is used to determine the correlation between two variables, and the value of the Pearson correlation coefficient varies from −1 to 1. The larger the absolute value of the Pearson correlation coefficient, the stronger the correlation between the two variables. A correlation coefficient between close to 1 or −1 indicates that the correlation between the two variables is strong, whereas a correlation coefficient close to 0 indicates that the two variables have no correlation.

Construction and validation of a prognostic model

The expression data of patients with OS from TCGA were treated as the training set, and the expression data of patients with sarcoma from TCGA were treated as the test set. The training set was used to build a prediction model through LASSO COX regression analysis. The risk score model obtained was as follows: risk score = (LPAR3 × 0.3996556 + AP003352.1 × 1.4405687). The test set was divided into low-risk and high-risk score groups based on the median value of the risk score. Receiver operating characteristic (ROC) curves were established to validate this model. The R package glmnet (4.1.6) was used for LASSO analysis.

Statistical analysis

Statistical analyses were performed using R v4.2.1. The predictive ability of the model was determined by ROC curve analysis. KM curves were used to compare the difference in survival between the two risk groups. GraphPad Prism software was used for the analysis of the assay results.

Screening of lncRNAs, microRNAs, and mRNAs related to osteosarcoma progression

During OS tumorigenesis and development, the expression level of tumor markers changes significantly. A majority of tumor markers will be used as oncogenes and antioncogenes to regulate the malignant progression of tumors. However, the exact role of the ceRNA axis, comprising of these genes, in osteosarcoma has not yet been clarified. To explore the internal mechanism, we performed differential expression analysis and univariate Cox analysis on mRNAs, lncRNAs, and miRNAs of osteosarcoma samples and precancerous samples, respectively. The results showed that there were 4,462 differentially expressed genes in mRNAs (up regulated genes: 2,333; down regulated genes: 2,129) and 1,242 survival related genes (oncogenes: 518; tumor suppressor genes: 724). A total of 2,348 differentially expressed lncRNAs (up regulated genes: 1,208; down regulated genes: 1,140) and 562 survival related genes (oncogenes: 319; tumor suppressor genes: 243) were detected. A total of 68 differentially expressed genes were found in microRNAs (upregulated genes: 37; downregulated genes: 31). The PCA diagrams showed the dispersion of differential analysis (Fig. S1A). The upregulated oncogenes and downregulated suppressor genes in osteosarcoma were extracted after obtaining the overlapping differentially expressed genes and survival related genes, including 186 protein-coding mRNAs and 45 lncRNAs (Figs. 1A–1C). The Venn diagrams showed the overlapped genes (Fig. S1B).

Enrichment of tumor progression related genes in osteosarcoma

To further explore the importance of genes in osteosarcoma, we performed GSEA enrichment analysis of mRNAs. The C5-GO enrichment analysis results demonstrated that most of the genes were involved in the malignant progression of osteosarcoma, including the proliferation function (Figs. 2A and 2B). Osteosarcoma is a highly malignant tumor originating from the bone marrow, which is characterized by the proliferation of osteoblastic precursor cells and the production of osteoid or immature bone. Traditional chemotherapy is the first choice for drug treatment of osteosarcoma. Chemotherapy drugs can significantly inhibit the proliferation of osteosarcoma cells, which also indicates that the inhibition of proliferation is the primary treatment for osteosarcoma. At the same time, the proliferation related pathways mTOR and MAPK signaling pathways play an pivotal role in the occurrence and development of osteosarcoma. In conclusion, the proliferation phenotype plays an important role in the development of osteosarcoma, so our research focuses on the proliferation pathway in the GO term analysis.

The construction of ceRNA network in osteosarcoma

According to the above results, the mRNAs, lncRNAs, and microRNAs obtained from survival analysis were used to construct ceRNA networks. The networks regulated the malignant progression of osteosarcoma, including 41 lncRNAs, 8 miRNAs, 36 mRNAs, and 85 groups of interaction axes (Fig. 3A). Cytoscape software was used to draw a network diagram of the genes regulation (Fig. S2). Among mRNAs, LPAR3 is important in the regulation of osteosarcoma. As a regulatory factor, LPAR3 participated in the regulation of the cell signal pathway, which could affect the proliferation and apoptosis of osteosarcoma cells. Therefore, the focus of our research was related to the ceRNA axis of LPAR3. The results showed that the AL358332.1/miR-141-3p, AC124798.1/miR-141-3p, LINC02280/miR-141-3p, AL135960.1/miR-141-3p, AP003352.1/miR-141-3p, AC073487.1/miR-141-3p, and AC084855.2/miR-141-3p axes could regulate LPAR3 (Fig. 3B).

AP003352.1 may regulate LPAR3 at the mRNA and protein levels

To explore the main lncRNA that regulated LPAR3, differential expression analysis and survival analysis were conducted on AL358332.1, AC124798.1, LINC02280, AP003352.1, AC073487.1, and AC084855.2, respectively (The count of AL135960.1 is too low). The results of differential expression analysis showed that four types of lncRNAs were significantly overexpressed in osteosarcoma (Fig. 4A). Therefore, these lncRNAs may be responsible for the overexpression of LPAR3 in osteosarcoma. In addition, the results of the univariate Cox regression analysis showed that the six lncRNAs were associated with a worse survival (Fig. 4B). Next, we conducted Pearson correlation analysis on six genes. And the results indicated that only AP003352.1 and AC084855.2 had significant correlations with LPAR3 (Fig. 4C). The AUC curve showed that the survival sensitivity of lncRNA AP003352.1 was significant (Fig. S3). To further screen the results, ASO was used to knock down AL358332.1, AC124798.1, LINC02280, AP003352.1, AC073487.1, and AC084855.2, respectively. Subsequently, the changes in LPAR3 caused by knocking down the relevant lncRNA were detected at the mRNA and protein levels (Figs. 4D, and 4E). The results revealed that the expression of LPAR3 was particularly regulated by AP003352.1, while no statistical difference was observed when other lncRNAs were knocked down. Thus, AP003352.1 may be the main lncRNA regulating LPAR3.

AP003352.1, miR-141-3p, and LPAR3 regulate the malignant progression of osteosarcoma

According to the prediction results reported above, AP003352.1/miR-141-3p/LPAR3 may be the ceRNA axis that regulated the osteosarcoma development. However, the exact role of the ceRNA axis in the malignant progression of osteosarcoma has not yet been clarified. Therefore, osteosarcoma samples were divided into two groups for differential expression analysis: LPAR3 high expression group and LPAR3 low expression group (AP003352.1 was for the same analysis). The results of differential expression analysis showed that there were many genes related to osteosarcoma proliferation (Figs. 5A and S4). To further evaluate the effect of the three genes on the proliferation of osteosarcoma, CCK8 assay and clone formation were conducted after respectively inhibiting the three genes. When the expression of AP003352.1, miR-141-3p and LPAR3 genes were changed, the value of OD450 detected by CCK8 and the number of clones were regulated accordingly. Therefore, the three genes can regulate the proliferation of osteosarcoma cells U2OS and MG63. The experimental results were consistent with the differential expression analysis results (Figs. 5B and 5C). At the end of the cloning experiment, we detected the mRNA levels of AP003352.1, miR-131-3p and LPAR3 by qPCR and the knockdown effects still existed (Fig. 5D). The above results indicated that all three genes could regulate the malignant progression of osteosarcoma.

AP003352.1 and miR-141-3p can regulate LPAR3 in osteosarcoma

To further explore the relationship between AP003352.1, miR-141-3p, and LPAR3, ASO was used to knock down AP003352.1, and the changes in miR-141-3p and LPAR3 were detected. The results showed that LPAR3 was typically downregulated, while miR-141-3p was upregulated (Fig. 6A). Then, after altering the expression of miR-141-3p by using overexpression mimics and inhibition inhibitors, we detected the changes in the expression of AP003352.1 and LPAR3. The results indicated that AP003352.1 and LPAR3 were negatively regulated by miR-141-3p (Figs. 6B and 6C). All the above results indicated that AP003352.1 and miR-141-3p can regulate the expression of LPAR3 in osteosarcoma, which was consistent with the trend of the ceRNA axis.

As the ceRNA of miR-141-3p, AP003352.1 regulates LPAR3 to affect the malignant progression of osteosarcoma

To determine whether AP003352.1 plays the same role as the ceRNA of miR-141-3p, the predicted binding sites of AP003352.1 and miR-141-3p were mutated complementarily (Fig. 7A). Luciferase reporter gene assay was utilized to identify the influence of miR-141-3p mutation on luciferase activity. The results showed that the luciferase activity of wild-type AP003352.1 was regulated by miR-141-3p, whereas that of the mutant-type AP003352.1 was not statistically significant (Fig. 7B). Then, we jointly predicted the enrichment of genes regulated by AP003352.1 and LPAR3 in the malignant progression pathway of osteosarcoma. The results showed that there was significant gene enrichment in the proliferation pathway (Fig. S5). Although, the genes were related to proliferation pathways in other cells, which suggested that the genes were indirectly related to osteosarcoma proliferation and had certain significance. The osteosarcoma cells and other types of cells may have some of the same proliferation regulation mechanism. To verify the results, we evaluated the cell proliferation in the control group and the osteosarcoma cells using miR-141-3p mimic after knocking down AP003352.1. The results of CCK8 assays showed that in the control group, AP003352.1 knockdown inhibited the proliferation of osteosarcoma cells U2OS, while this inhibition was not observed in the mimic group (Fig. 7C). Besides, the results of clone formation assays indicated that AP003352.1 knockdown also inhibited the colony formation ability of osteosarcoma cells U2OS in the control group and in the mimic group, the inhibition of colony formation ability was disappeared (Fig. 7D). Besides, we also detected the cell proliferation in the control group and the osteosarcoma cells after knocking down LPAR3. The results showed that AP003352.1 knockdown could inhibit the proliferation of osteosarcoma cells U2OS in the control group, but the inhibitory effect was not observed in LPAR3 knockdown group (Figs. 7E and 7F). All above results showed that AP003352.1 regulated LPAR3 and affected the malignant progression of osteosarcoma, similar to the ceRNA of miR-141-3p.

The AP003352.1/miR-141-3p/LPAR3 axis can be a better biomarker for predicting osteosarcoma prognosis

Our previous research results indicated that the three genes were correlated with the malignant progression of osteosarcoma (Figs. S6A and S6B). However, the malignant progression of osteosarcoma was co-regulated by a multi gene network. Therefore, there were some limitations in the diagnosis of osteosarcoma by using a single gene as the biomarker. A prognosis model based on AP003352.1/LPAR3 expression was established through LASSO regression analysis (Figs. 8A–8D, hazard system = LPAR3 × 0.3996556 + AP003352.1 × 1.4405687). In addition, higher risk coefficients were associated with poorer survival, and the AUC values in both the training set and validation set were significantly higher than those in the single gene model (Figs. 8E–8G and Figs. S6C–S6F). Thus, we constructed a novel prognosis model by AP003352.1, miR-141-3p and LPAR3 genes in osteosarcoma.