Multi-omics analysis identifies potential mechanisms of AURKB in mediating poor outcome of lung adenocarcinoma

Overexpression of AURKB in LUAD patients

TCGA-LUAD cohort was used to compare the expression of AURKB in lung adenocarcinoma patients. As shown in Figure 1A, AURKB was overexpressed in LUAD tumor tissues (513 cases) when compared with 59 normal tissues (rank-sum test, p < 0.001). We further examined the correlation between AURKB expression and the clinical features of the 513 lung adenocarcinoma patients. For TNM system, we found that AURKB was significantly low-expressed in tumors with smaller size (T1 stage) when compared with T2 (p<0.001) or T3-4 (p=0.028) (Figure 1B). Patients without lymph node invasion (N0) showed low expression of AURKB than those with N1 (p=0.011) or N2-3 (p=0.001) (Figure 1C). No significant differences among patients with or without metastasis (p=0.123, Figure 1D) were detected, which might be due to the limited number of metastatic patients in the LUAD cohort. However, AURKB showed obvious variations in the patients’ staging: lower expression in stage I against stage II (p=0.003) or III/IV (p<0.001), as shown in Figure 1E. We also found that AURKB was noticeably high-expressed in patients with a history of smoking (p=0.002, Figure 1F), and that elder or female patients tend to have a relatively lower expression of AURKB (p=0.041, 0.0028, respectively, as shown in Table 1). We additionally compared 7 microarrays and sequencing data sets using the lung cancer explorer (LCE) database, and the results confirmed that AURKB expression was up-regulated in lung adenocarcinoma (p=3.1e-10, Figure 1G).

AURKB indicates a worse outcome in LUAD

As AURKB is high-expressed in lung adenocarcinoma tissues and correlates with more advanced tumor stage, we hypothesized that AURKB might be a prognostic factor in LUAD. Thus, the expression distribution of AURKB and prognosis of each sample were examined. We observed that the expression of AURKB of each individual was (Figure 2A) concentrated, with a median value of 14.33. Based on this, lung adenocarcinoma patients were then divided into AURKB^high and AURKB^low groups, with 245 samples each group. As shown in Figure 2B, 2C, higher expression of AURKB contributed to significantly worse overall survival (OS, p=0.002) and progression-free survival (PFS, p=0.013), indicating that AURKB is associated with poor prognosis in patients with lung adenocarcinoma. The data of multivariate cox analysis also showed that AURKB and N staging could be two independent prognostic factors of LUAD (HR, 95%CI: 1.672, 1.045-2.674 for AURKB, 1.897, 1.13-3.184 for N staging, Figure 2D), irrespective of T stage, M stage, smoking history and EGFR mutation status. The outcome data of 21 data sets involving 2912 patients with LUAD further confirmed that AURKB could be a prognostic factor in LUAD (HR, 95%CI: 1.22, 1.13-1.32, Figure 2E).

To evaluate the potential of AURKB as a prognostic marker in blood samples, we obtained the data of 80 LUAD tumor samples and 80 normal control samples from the GEO database (GSE20189). The original data of the chip were downloaded, and the R software package Affy was used to perform data conversion and RMA standardization for obtaining a gene expression profile data set. The expression of AURKB in tumor samples was found to be higher than that of control samples with significant margins (p=0.072, Supplementary Figure 1A). Further analysis of the changes in the expression of AURKB at the early and advanced stages of LUAD showed that AURKB was significantly high-expressed in tumors more advanced after Stage II (p=0.043, Supplementary Figure 1B). Angiogenic factors are key markers of tumor progression. The ssGSEA method was used to evaluate the scores of angiogenic factors in the samples, and the data indicated that the scores of angiogenic factors in tumor samples were significantly higher than in normal samples (p=0.00012, Supplementary Figure 1C). In LUAD patients, the expression of AURKB was greatly positively correlated with the angiogenic factor score (r=0.52, p<0.001, Supplementary Figure 1D). The above results indicated that the expression of AURKB was of a certain consistency in blood samples and tumor tissue samples.

Genome-wide mRNA profiles associated with AURKB expression

To identify the key genes involved in the AURKB-mediated poor survival of LUAD and related biological functions or pathways, we compared the expressions of differentially expressed genes (DEGs) between AURKB^high group and AURKB^low group and examined their correlations with AURKB using the T-test and the rank-sum test, respectively (mRNAs were excluded if over 50% values are empty or zero). Only genes with correlation coefficients larger than 0.3, fold change >2, and the p-value < 0.01 either in the T-test or the rank-sum test were considered as DEGs. As visualized by volcano plot and heatmap (Figure 3A, 3B), 348 DEGs incorporating 270 upregulated genes and 78 downregulated genes were found. The top 10 upregulated and down-regulated genes identified by fold changes were listed in Table 2. The most correlated mRNA with AURKB was FAM64A, which is a marker of cell proliferation, and has been found to be associated with the aggressive growth of pancreatic cancer and breast cancer [16, 17]. Moreover, study indicated that FAM64A may also regulate Th17 cells and promote inflammation-associated cancers [18]. PRAME, an antigen high-expressed in various types of cancers, is involved in cell apoptosis, differentiation, proliferation, and metastasis [19], and it had the greatest fold change in the two groups. CDCA3 is the most significantly differentiated mRNA, and is already proven to be a prognostic factor in NSCLC patients [20]. Additionally, we further used ClusterProfiler [21] R-package to run KEGG enrichment analysis, and found that DEGs were enriched in 19 pathways (p<0.05), especially in cancer-related pathways including cell cycle, DNA replication, homologous recombination, and p53 signaling pathway (Figure 3C). These pathways were additionally validated by KEGG pathway enrichment analysis using GSEA. Meanwhile, we found some other pathways (24 pathways in total), such as purine metabolism, pyrimidine metabolism, spliceosome, and proteasome, were related to high expression of AURKB (Figure 3D). Since AURKB could mediate EGFR TKI resistance as above-mentioned, we also analyzed the relationship between 19 pathways and EGFR pathway by calculating the score of each patient in these pathways using ssGSEA. The results demonstrated that most of the 19 pathways were significant correlated with the EGFR pathway (89.5%) (Supplementary Figure 2A). However, for the AURKB^low group, no pathways met our criteria that both p-value and FDR q-value should be less than 0.05. Aberrant DEGs and enriched pathways identified in this study were closely associated with cancer formation and therapy resistance, which was consistent with our previous finding that AURKB mediated poor outcome for LUAD patients and might explain for the internal mechanism of action.

Connections between microRNAs and AURKB

We also found certain DEGs were enriched in the miRNA pathways of LUAD, thus, genome-wide miRNA profile analysis was performed. In a further Pearson correlation analysis with AURKB, 926 miRNAs with at least 30% sample expression values were included. A total of 54 miRNAs met the threshold, with a FDR of less than 0.05 and a correlation coefficient of above 0.3. Specifically, 50 of these miRNAs were positively correlated with AURKB expression, while 4 were negatively correlated with AURKB expression. The distribution relationship of correlation coefficients with FDR is shown in Figure 4A. As displayed by the clustered heatmap (Figure 4B), only a small number of miRNAs were negatively correlated with AURKB, while most of the miRNAs were positively correlated with AURKB. Among them, the top miRNA with the highest positive correlation coefficient was miR-130b-3p. Previous research has been reported that miR-130b-3p present in exosomes contributes to CXCL12/CXCR4-induced colorectal cancer metastasis into liver [22] and downregulate PIEZO2, thereby activating the Hedgehog signaling pathway and leading to worsening breast cancer prognosis [23]. Moreover, miR-421 is previously identified as a prognostic biomarker indicating promoted tumor progression of non-small cell lung cancer (NSCLC) [24]. miR-29b-2-5p and miR-101-3p were negatively correlated with AURKB. From reviewing previous studies, miR-29b-2-5p is reported to target CBI-b to inhibit the proliferation of pancreatic cancer cells [25], and miR-101-3p could hinder the growth and metastasis of NSCLC through attenuating MALAT-1 mediated PI3K/AKT signal pathway [26].

Correlated miRNAs were analyzed by the T-test and rank-sum test to compare the expression changes between AURKB^high and AURKB^low groups. 49 miRNAs were screened based on FDR < 0.05 and fold change > 1.5. Then mirWallk database was used to predict the target genes of these differential miRNAs and to construct mRNA-miRNA regulatory network, as shown in Figure 4C. Further analysis of the relationship between miRNA in the network and EGFR pathway showed that 82% of miRNAs predicted were significantly correlated with EGFR pathway (Supplementary Figure 2B). Based on the scores calculated by the degree algorithm, the following top 10 hub genes were screened: miR-93-5p, RAB3B, miR-17-5p, miR-20a-5p, WNK3, miR-7974, miR-6783-3p, PRR11, GINS4 and miR-128-3p (Table 3). miR-93-5p, which was reported to promote the proliferation of NSCLC cells and is indicative of a poor prognosis, was found to enhance the transcription of several oncogenic genes, such as HMGA2, E2F2, KIF23, POLQ, PRR11, GINS4, WNK3, LYPD6, KIF14, RAB3B and downregulate the expressions of tumor suppressors such as ATP1A2 and CYBRD1. Taken together, these observations provided rational evidence for the prognostic role of AURKB in lung adenocarcinoma.

To further validate the association of identified miRNAs with AURKB in vitro, we tested the protein expression levels of AURKB among 4 lung adenocarcinoma cell lines, including A549, H1299, PC9 and H1975 (Supplementary Figure 3A). Lentiviral shRNA-mediated ARUKB inhibition in AURKB high-expressed PC9 cells and lentiviral AURKB-mediated overexpression in AURKB low-expressed H1975 cells were confirmed by western blotting and quantitative PCR (Supplementary Figure 3B, 3C). Knockdown of AURKB in PC9 cells significantly downregulated the expressions of miR-93-5p, miR-17-5p, miR-20a-5p, miR-130b-3p and upregulated the expressions of miR-101-3p (Figure 4D). While overexpression of AURKB in H1975 cells significantly upregulated the expressions of miR-93-5p, miR-17-5p, miR-20a-5p, miR-130b-3p, and miR-128-3p (Figure 4E).

Whole-genome methylation profiles associated with AURKB

DNA methylation is an important epigenetic mechanism regulating gene expressions through three DNA methyltransferases (DNMT1, DNMT3A, and DNMT3B) and affects cancer cell behaviors. We analyzed differences in the transcription of the three methyltransferases, and as shown in Figure 5A–5C, expressions of the three were significantly higher in AURKB^high than AURKB^low groups. 476 samples with methylation profiles were included to further examine the methylation differences between the AURKB^high and AURKB^low groups. Methylation sites were excluded if more than 30% of the samples showed no values or happened to be the cross-reactive CpG sites in the Illumina Infinium Human Methylation 450 microarray. After filling up missing values using the KNN method in the R package ‘impute’, a total of 195,622 methylation sites were obtained and then compared with the R package ‘limma’, according to the thresholds of fold change > 1.3 and FDR < 0.05. The volcano map in Figure 5D showed that the hypermethylation sites were significantly more than hypomethylation sites (1900 and 411, respectively). Heatmap of these differential methylation sites is shown in Figure 5E.

By examining the position distribution of the 2311 differential methylation regions (DMR) around the CpG islands, as shown in Figure 5F, all the positions had more hypermethylated DMRs, although the differences were not significant in N_Shelf or S-Shelf (Island: P<0.001, N_Shore: P=0.0019, S_Shore and Others: P<0.001). As expected, most hypermethylated DMRs centered on the CpG islands. Finally, we determined the distance of methylation sites to the transcription starting site (TSS), as presented in Figure 5G, most hypomethylated DMRs and hypermethylated DMRs were on the upstream of TSS (-600 – +200), and hypermethylated DMRs fell to -200 – +200 regions. Generally, hypermethylated sites tended to distribute around TSS, while hypomethylated sites tended to distribute on the upstream of TSS.

We extracted the corresponding genes of these DMRs based on the TSS and ran a KEGG enrichment analysis using R package ClusterProfiler. The data demonstrated that genes corresponding to hypermethylation were enriched into 36 pathways, such as cAMP signaling pathway, axon guidance, viral carcinogenesis, and chemokine signaling pathway (p<0.05, Figure 5H). However, the relative hypomethylated genes were only enriched into two pathways, which were olfactory transduction and viral protein in interaction with cytokine and cytokine receptors (Figure 5I).

Tumor suppressors DLC1 and HLF negatively correlated with ARUKB and its associated miRNAs and methylation modification

As we have obtained the genome-wide mRNAs, miRNA-regulated genes, and methylation-regulated genes profiles, we then went on to investigate whether these three profiles shared genes in common that may have close crosstalk with AURKB to promote LUAD aggressiveness. As shown in Figure 6A, two genes, DLC1 and HLF, were present in all the three profiles. The differential mRNA expression levels of both genes were further confirmed using paired lung adenocarcinoma tumor and adjacent tissues in TCGA (p<0.001 for both genes, Figure 6B). Gene correlation analysis using GEPIA web-tool (http://gepia.cancer-pku.cn/) also confirmed that AURKB was significantly negatively correlated with DLC1 and HLF (R= -0.29, -0.34, respectively, Figure 6C). Survival analysis of included LUAD patients from the TCGA database also showed an apparent difference between HLF^high and HLF^low patients (p=0.0029), and DLC1 was determined to be a protective factor to the survival of patients with LUAD (p=0.078) (divided by the median value, Figure 6D). LCE online survival analysis demonstrated that low expressions of DLC1 and HLF were associated with a shorter survival time of LUAD patients (HR[95%CI], 0.82[0.77; 0.87], 0.77[0.72;0.82], respectively, Figure 6E).

Data acquisition and preprocessing

The expression profiles of mRNAs and miRNAs were sequenced by RNA-seq and 450K genome-wide methylation profile. LUAD clinical data, including 513 tumor samples and 59 para-cancer tissues in the Cancer Genome Atlas (TCGA) were downloaded from Genomic Data Commons Data Portal (portal.gdc.cancer.gov, 05/23/19). According to the clinical data, samples acquired from patients with follow-up time more than 30 days were selected to match the expression profiles of mRNAs, miRNAs, and methylation information for final analyses. The detailed statistical description is listed in Table 4. The cancer samples were subsequently divided into AURKB^high expression and AURKB^low expression subgroups based on the median value and were used for further analyses. Lung Cancer Explorer (LCE, http://lce.biohpc.swmed.edu/lungcancer/) online analysis was used to compare the expression levels between tumor and normal tissues and survival changes based on AURKB expressions. All the results were calculated by the meta-analysis of multiple LUAD datasets from TCGA and Gene Expression Omnibus (GEO).

Cell lines

LUAD cell lines (A549, H1299, H1975, and PC9) were purchased from the American Type Culture Collection (ATCC) Cell Bank and cultured in RPMI 1640 (Gibco, catalog: 31800-022) at 37° C in a 5% CO₂ atmosphere. All the media were supplemented with 10% fetal bovine serum (serana, catalog: S-FBS-EU-015), 100 U/ml penicillin and 100 μg/ml streptomycin.

Quantitative RT-PCR

Total RNA from cell pellets was extracted using the Total RNA Miniprep Kit (AxyPrep, USA) according to the manufacturer’s protocol. The quantity and quality of extracted RNA were assessed with a Nanodrop 2000c (Thermo Scientific, USA). Then, complementary DNA was synthesized from 500 ng of total RNA using PrimeScript™ RT Master Mix kit (Takara, Japan). For the stem-loop reverse transcription of microRNAs, PrimeScript™ RT Reagent Kit (Takara, Japan) was used to syntheze cDNA with miRNA-unique stem-loop primers. PCR analysis was performed in a Applied Biosystems 7500 Fast Real-Time PCR System (Foster City, USA) with SYBR Premix Ex Taq II (Takara, Japan). The expressions of miRNAs and mRNAs were calculated using the 2^-ΔΔCT method, with U6 and ACTB served as normalization controls, respectively. All primers used are listed in Table 5.

Western blot

Total protein were extracted from cell pellets using RIPA buffer. Protein concentration was determined by BCA assay (Pierce, Rockford, IL, USA) using bovine serum albumin as a standard. 30 ug protein samples were separated in 10% SDS-PAGE gel and transferred to a PVDF membrane, which was blotted with antibodies (primary antibodies rabbit anti-AURKB (Cat# ET1610-25, 1:1000), rabbit anti-GAPDH (Cat# ER1706-83, 1:3000), rabbit anti-beta-actin (Cat#R1207-1, 1:3000)), then rinsed and incubated with secondary IgG-HRP antibodies goat anti-rabbit (Cat# HA1019, 1:3000) and goat anti-mouse (Cat#HA1020, 1:3000) according to the manufacturer’s instructions. All antibodies came from HuaBio in China.

AURKB overexpression or inhibition in lung adenocarcinoma cells

Briefly, pLenti-CMV-GFP-Puro-AURKB and pPLK/GFP+Puro-AURKB shRNA were constructed on the backbone of pLenti-CMV-GFP-Puro and pPLK/GFP+Puro plasmids, respectively, and each was then co-tranfected with lentiviral packaging plasmids pMD2.G and psPAX2 (Addgene, USA) into HEK 293T cells to generate lentivirus particles, followed by the transfection of the particles into target cells. The cells were incubated with lentivirus for 24 hours (h) and then replaced with fresh medium. When cells stably grew to proliferative status, cell pellets were collected for protein and total RNA extraction.

Statistical analysis

The chi-squared test was used to examine the correlation between AURKB expression and clinical or pathological variables. The endpoint event overall survival (OS) was defined as the time starting from the occurrence of LUAD to death due to any cause. Progression-free survival (PFS) was defined as the time starting from the onset of LUAD to progression after the first-line treatment or death. Kaplan-Meier analysis with the log-rank test was applied to compare the survival differences between two groups of patients, while multivariate cox regression analysis was used to examine potential independent factors related to survival. Student’s t-test and false discovery rate (FDR) were used to identify differences in mRNAs, miRNAs, and methylation data between AURKB^high and AURKB ^low expression groups. All the analyses were performed using the R 3.6 software packages.