The circRNA expression profile of GSE101586 was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which included five pairs of matched LAC tissues (GSM2706428, GSM2706429, GSM2706430, GSM2706431 and GSM2706432) and adjacent non-cancerous tissues (GSM2706423, GSM2706424, GSM2706425, GSM2706426 and GSM2706427). The tissues were obtained from five female patients with LAC with no smoking history, and the physiology of the tissues were confirmed by pathologists. The mRNA data of the LAC tissues was downloaded from The Cancer Genome Atlas (https://portal.gdc.cancer.gov).

In total, 60 cases of LAC tissues in the tissue bank were obtained from the patients who underwent surgery in the Fourth hospital of Hebei Medical University (Shijiazhuang, China) between January 2016 and June 2016. A total of 20 pairs of LAC and adjacent normal tissue (5-cm from the tumor tissue) were collected from the above 60 cases. The clinicopathological characteristics and survival status of the patients with LAC were obtained from the follow-up data. (The duration of follow-up was between January 2016 and June 2021). The inclusion criteria were: i) All patients were diagnosed LAC through pathology and ii) All patients did not receive any therapy before surgery. The clinicopathological data for the patients are presented in Table SI. The human tissues were obtained with written informed consent, and the present study was approved by The Clinical Research Ethics Committee of The Fourth hospital of Hebei Medical University (approval no. 2021KY157).

The limma package in R (v3.4.1; http://bioconductor.org/packages/release/bioc/html/limma.html) was used to identify the differentially expressed circRNAs between the LAC tissues and the adjacent non-cancerous tissues from the GEO101586 dataset using the criteria of P<0.05 and |Log₂ fold change (FC)|>1. The ‘edger’ package (v3.52; http://bioconductor.org/packages/edgeR/) was used to screen for the differentially-expressed mRNAs using thresholds of |Log₂FC|>3 and P<0.01. The miRNAs-targeted mRNAs were predicted based on the miRWalk software (v2.0; http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/).

LAC cell lines H1299 and H1975 were obtained from Chinese Academy of Sciences Cell Bank were cultured in RPMI-1640 media (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.), 100 U/ml penicillin and 100 µg/ml streptomycin. Cells were cultured at 37°C under 5% CO₂.

Total RNA from 80% confluent cells and tissues were isolated by TRIzol^® Reagent (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocols. The cDNA was prepared from total RNA according to the GoScript Reverse Transcription System (Promega Corporation) protocol. qPCR analysis was performed using the Promega GoTaq qPCR Master Mix (Promega Corporation). GAPDH was used as the mRNA internal reference. The primer sequences for all qPCR reactions are shown in Table SII. The PCR cycling conditions was as follows: Denaturation at 95°C for 15 sec, annealing at 58°C for 30 sec and extension at 72°C for 30 sec. Three repeat wells were set. The relative expression levels were calculated using the 2^−ΔΔCq method (17). The sequencing of PCR products were obtained from Sangon Biotech Co., Ltd. RNase R treatment was carried out for 15 min at 37°C using RNase R 3 U/mg.

siRNA for negative control and siRNAs targeting the junction sites of circ-0006220, circ-0072088, and circ-0001666 were designed and synthesized by Guangzhou RiboBio Co., Ltd. The sequences of these siRNAs are listed in Table SIII. H1299 cells were transfected with above siRNAs when reached 70–80% confluence by using HiPerFect Transfection Reagent (Qiagen GmbH) based on the manufacturer's instruction. The final siRNA concentrations were 50 nmol/l. Following transfection, cells were cultured at 37°C for 72 h.

H1299 cells were first seeded at a cell density of 5×10³ cells/well into 96-well plates. After attachment overnight, CCK-8 assays were performed after 24, 48 and 72 h. In brief, 10 µl CCK-8 solution (cat. no. HY-K030; MedChemExpress) was added into each well. After 1 h incubation at 37°C, the absorbance readings for each well were performed at 450 nm using the microplate reader (Tecan Group, Ltd.).

The migration and Matrigel invasion assays were performed in 6.5-mm Transwell chambers for migration assays or Matrigel pre-coated (at 37°C for 24 h) chambers for invasion assays based on the manufacturer's protocol (BD Biosciences). The cell suspensions of different groups were added to the upper chambers at a density of 5×10⁴ cells/well and incubated for 24 h. The migratory and invasive cell numbers were then quantified in five random fields per chamber under the inverted microscope at ×100 magnification.

Prediction of interactions between circRNAs and miRNAs was performed using the CircInteractome database (https://circinteractome.irp.nia.nih.gov/). miRNAs with a context score percentile ≥98 were eventually selected.

The ceRNA regulatory network was established according the possible interactions among the four differentially-expressed circRNAs, eight miRNAs predicted to be targeted by these circRNAs and the 232 overlapped mRNAs from the list of predicted target genes and the upregulated genes in LAC. Cytoscape 3.7.1 software (v.3.7.1; http://cytoscape.org/) was used to visualize the established ceRNA regulatory network.

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) software (version 11.5; http://string-db.org/) and visualized using the Cytoscape 3.7.1 software (https://cytoscape.org/). Additionally, suspected hub genes with high predicted degrees of interactions from the PPI network were screened using the CytoHubba plugin in Cytoscape.

The mRNA expression levels of the hub genes and their association with survival were assessed using the web-based GEPIA database (http://gepia.cancer-pku.cn/) with the settings of P≤0.05 and |Log₂FC|≥1.

The 5-µm, paraffin-embedded slides were deparaffinized in xylene, rehydrated using a decreasing alcohol gradient and washed with 1X PBS (pH 7.2) three times at 5 min each. The sections were then heated in a microwave oven for 5 min in 10 mmol/l Na-citrate buffer (pH 6.0) for antigen retrieval and washed with 1X PBS. The sections were immersed in 0.3% hydrogen peroxide in methanol for 20 min to suppress endogenous peroxidase activity. After further washing with 1X PBS, the sections were incubated in 10% normal goat serum (cat. no. SP-9001 or SP-9002; OriGene Technologies, Inc.) at room temperature in a humidified chamber for 30 min to prevent nonspecific immunoglobulin binding. The sections were then treated with the 1:100-diluted KIF2C antibody (cat. no. TA503320; OriGene Technologies, Inc.), KIF18B (cat. no. ab121798; Abcam), CKAP2L (cat. no. ab122617; Abcam), PLK1 antibody (cat. no. TA500383; OriGene Technologies, Inc.), MELK antibody (cat. no. ab129373; Abcam), BIRC5 antibody (cat. no. TA301427; OriGene Technologies, Inc.) at 4°C overnight. Normal IgG in place of the primary antibody served as the negative control. A streptavidin-biotinylated HRP-based detection system was used to reveal specific binding. The sections were counterstained for 2 min at room temperature with hematoxylin for light microscopic review and evaluation. The expression was ranked on the sum of intensity and area from 0 to 7: 0–2, negative expression; 3–7, positive staining (of those, 3–4, weak positive expression; and 5–7, strong positive expression). Staining intensity was graded as follows: 0 for no staining; 1 for mild staining; 2 for moderate staining; and 3 for intense staining. The staining area was scored as follows: 0 for no staining; 1 for 1–25% area; 2 for 26–50% area; 3 for 51–75% area; and 4 for 76–100% area.

Statistical analysis was performed using SPSS 22.0 software (IBM Corp). The data were presented as the mean ± standard deviation and measured using the Student's t-test (paired t-tests were used in Fig. 1C, and unpaired t-tests were used in Fig. 2C-F). ANOVA analysis followed by Dunnett's post hoc test was used for 2G-I. Chi-square test was used to analyze the association of circRNA expression with the hub protein expression (Fig. 7C). Kaplan-Meier analysis was used to evaluate the overall survival and log-rank test was performed to estimate the differences among the various groups. P<0.05 was considered to indicate a statistically significant difference.

The GSE101586 dataset from GEO was analyzed using the limma package in R (P<0.05 and |Log₂FC|>1). In total, 47 upregulated and 21 downregulated circRNAs were identified in LAC tissues compared with adjacent non-cancerous tissues (Table I). After tightening the parameters to P<0.01 and |Log₂FC|>2, four circRNAs, hsa_circ_0043278, hsa_circ_0006220, hsa_circ_0072088 and hsa_circ_0001666, were identified to be significantly upregulated in LAC tissues compared with adjacent non-cancerous tissues (Fig. 1A). The basic structural profiles of these four circRNAs are shown in Fig. 1B. All four of these circRNAs are derived from exons of their parent genes. Among them, three circRNAs (hsa_circ_0006220, hsa_circ_0072088, hsa_circ_0001666) were confirmed to be upregulated in LAC tissues compared with those in the corresponding normal tissues by RT-qPCR (Fig. 1C).

To further explore the characteristics of three differentially-expressed circRNAs, two sets of primers for each circRNA were first designed. Divergent primers were expected to amplify the circular forms of the RNAs, whereas convergent primers were expected to amplify the linear forms of RNAs. Using cDNA and genomic DNA (gDNA) from LAC tissues as templates, all three of the differentially-expressed circRNAs were amplified by the divergent primers in cDNA, but no amplification products were observed in the gDNA samples (Fig. 2A). All three linear forms of the circRNAs were amplified in both cDNA and gDNA samples. The sequences of the PCR products from all three circRNAs were verified by sequencing (Fig. 2B). By using RT-qPCR, it was verified further that all three of these circRNAs were resistant to RNase R, whilst the levels of their corresponding linear forms were significantly reduced after RNase R treatment (Fig. 2C).

Subsequently, the expression of these three circRNAs were observed in the 60 LAC tissue samples. In addition, higher expression levels of hsa_circ_0006220, hsa_circ_0072088 and hsa_circ_0001666 indicated poor prognosis in patients with LAC (Fig. 2D). To further explore the biological functions of these differentially-expressed circRNAs, their expression was then measured in the LAC cell lines (Fig. 2E). siRNA-mediated knockdown of these circRNAs was found to suppress proliferation, migration and invasion in LAC cells (Fig. 2F-I). Taken together, these results suggest that hsa_circ_0006220, hsa_circ_0072088 and hsa_circ_0001666 can serve oncogenic roles in LAC progression and have potential as biomarkers and therapeutic targets for LAC.

Accumulating evidence has demonstrated that several circRNAs derived from exons can exert biological functions by acting as ‘sponges’ or ‘decoys’ to sequester miRNAs. To assess if these three aforementioned circRNAs performed similar functions in LAC, miRNAs that can potentially interact with these three circRNAs were screened using the CircInteractome database. In total, eight miRNAs (miR-520f, miR-1261, miR-1270, miR-620, miR-188-3p, miR-516b, miR-940 and miR-661) with context score percentiles ≥98 were selected (Fig. 3A). The predicted binding sites between the circRNAs and miRNAs of interest are shown in Fig. 3B.

A total of 7,829 target genes of the eight aforementioned miRNAs were obtained from the miRWalk software. Additionally, 795 upregulated genes in LAC were obtained from TCGA, according to the thresholds of false discovery rate value <0.01 and Log₂|FC|>3 (Fig. 4A). By overlapping the predicted target genes and the upregulated genes in LAC, 232 target genes that may serve key roles in LAC were identified (Fig. 4B). From this, a circRNA-miRNA-mRNA network was constructed by integrating the circRNA-miRNA interactions and miRNA-mRNA interactions (Fig. 4C), which provided a preliminary insight into the links among the three differentially-expressed circRNAs (hsa_circ_0006220, hsa_circ_0072088 and hsa_circ_0001666), their eight suspected target miRNAs (miR-520f, miR-1261, miR-1270, miR-620, miR-188-3p, miR-516b, miR-940 and miR-661) and the 232 mRNAs.

After removing the unconnected nodes, a PPI network consisting of 158 nodes and 287 edges was constructed to view the interaction among the 232 target genes (Fig. 5A). Subsequently, a pivotal module of six hub genes, kinesin family member (KIF) 2C, KIF18B, maternal embryonic leucine zipper kinase (MELK), baculoviral IAP repeat-containing 5 (BIRC5), polo-like kinase 1 (PLK1) and cytoskeleton-associated protein 2-like (CKAP2L), were identified from the PPI network using the CytoHubba plugin in Cytoscape (Fig. 5B).

The six hub genes in the PPI network were next evaluated for their expression levels and prognostic value using the GEPIA database. As shown in Fig. 6A-F, the expression of these six hub genes in LAC tissues was significantly higher compared with that in the normal tissues. All of these six hub genes exhibited their potential in the prediction of overall survival based on their expression. High expression levels of each hub gene were associated with poorer overall survival in patients with LAC (Fig. 6A-F).

A circRNA-miRNA-hub gene sub-network was then built to delineate the links among the differentially-expressed circRNAs, target miRNAs and hub genes (Fig. 6G). Of note, the hsa_circ_0006220/hsa-miR-520f/KIF2C, has_circ_0001666/hsa-miR-940/KIF18B, has_circ_0001666/hsamiR-940/CKAP2L, has_circ_0001666/hsa-miR-661/CKAP2L, has_circ_0001666/hsa-miR-188-3p/PLK1, has_circ_0072088/hsamiR-1270/MELK and has_circ_0072088-hsa-miR-1261-BIRC5 regulatory axes were found from this network. Taken together, these results suggest that hsa_circ_0006220, hsa_circ_0072088 and hsa_circ_0001666 may serve oncogenic roles by regulating the expression of these six hub genes.

To further validate the findings, immunohistochemistry was performed in the 60 LAC samples. The statistical results showed that high expression levels of KIF2C, KIF18B, MELK, BIRC5, PLK1 and CKAP2L were associated with the poorer survival in patients with LAC (Fig. 7A and B), which was consistent with the results of GEPIA database analysis. In addition, the expression of hsa_circ_0006220 was positively correlated with the expression of KIF2C, whereas the expression of has_circ_0001666 was positively correlated with the expression of KIF18B, MELK and BIRC5. The expression of has_circ_0072088 was positively correlated with the expression of PLK1 and CKAP2L (Fig. 7C).