Atherosclerosis-related biomarker PABPC1 predicts pan-cancer events

Serum samples

Nineteen serum samples were obtained from 19 patients with AS-confirmed TIA who were treated at the hospital of Japan. The remaining four AS-related TIA serum samples were obtained from hospital of China. All patients were diagnosed with TIA according to the latest TIA diagnostic criteria of the American Stroke Association.35 36 All of them understood the purpose and steps of the experiments, and agreed to participate in the study before sera collection. Written informed consent was obtained from all participants. Thereafter, peripheral blood (5 mL) was collected from each patient and centrifuged at 3000×g for 10 min after being kept at room temperature for 2 hours. The supernatant was regarded as a serum sample and stored in a −80°C refrigerator until further use. Repeated freezing and thawing was avoided.

Immunoscreening of cDNA expression clones

Sera from 23 patients with AS-related TIA were subjected to 3 rounds of SEREX immunological screening of a human aortic endothelial cell cDNA expression library (Stratagene, USA; concentration: 2.0×10⁹ pfu/mL). XL1-Blue strain (2 µL) was extracted and added to 10 mL of Lysogeny Broth (LB) medium containing 12.5 g/L tetracycline. It was placed in an incubator at 30°C and a rotation speed of 210 rpm for overnight incubation. A concentration of 2–5×10⁴ pfu containing the human aortic endothelial cell cDNA library was mixed with 600 µL Escherichia coli(E-coli) XL1-Blue bacterial solution (A600=0.5). The plate was incubated until plaques formed, and preheated upside down at 37°C for 1 hour. The nitrocellulose (NC) membranes were pretreated with 10 mmol/L isopropyl β-D-thiogalactoside (IPTG) solution for 30 min and then the surface of the plate was covered. The NC membranes were pretreated with 10 mmol/L IPTG solution for 30 min, covered on the surface of the plate and incubated upright at 37°C for 2.5 hours to let the plaques be transferred to NC membranes. The NC membranes were gently peeled off and washed three times with Tris Buffered Saline with Tween 2o (TBS-T) buffer. After the NC membrane was blocked with 10 mg/g bovine serum albumin (Nacalai Tesque, Japan), it was transferred to TIA patient serum (primary antibody), which had been diluted to 1:2000, and shaken slowly for 1 hour. The NC membrane was washed with TBS-T again and soaked in 1:5000 diluted goat anti-human immunoglobulin (Ig) G H&L (HRP)(secondary antibody) (Jackson Immunology Research Laboratory, Japan) for 1 hour. The NC membrane was immersed in nitro blue tetrazolium and 5-bromo-4-chloro-3-indolylphosphate solution in dark for colour reaction to occur. To obtain purified positive single clones, the second and third rounds of screening were performed by diluting titre of cDNA expression library according to the first round of screening.

TBS-T buffer. After the NC membrane was blocked with 10 mg/g bovine serum albumin (Nacalai Tesque, Japan), it was transferred to TIA patient serum (primary antibody), which had been diluted to 1:2000, and shaken slowly for 1 hour. The NC membrane was washed with TBS-T again and soaked in 1:5000 diluted HRP-labelled goat anti-human immunoglobulin (Ig) G (secondary antibody) (Jackson Immunology Research Laboratory, Japan) for 1 hour. The NC membrane was immersed in nitro blue tetrazolium and 5-bromo-4-chloro-3-indolylphosphate solution in dark for colour reaction to occur. To obtain purified positive single clones, the second and third rounds of screening were performed by diluting titre of cDNA expression library according to the first round of screening.

Intraclonal excision and restriction enzyme digestion identification

A mixture of 100 µL XL1-Blue (A600=1.0), 125 µL positive monoclonal colony dilution and 0.5 µL ExAssist helper phage (>1.0×10¹⁰ pfu/mL) were prepared. The mixture was added to 3 mL LB medium and incubated at 30°C with shaking for 2.5–3.0 hours. The mixture was heated and incubated in a dry thermostatic metal bath at 65°C–70°C for 20 min, and then centrifuged at 1000×g for 15 min. The pBluescript phagemid was obtained from the supernatant, transfected with E-coli SOLR strain and spread on a petri dish, which was incubated upside down overnight at 37°C. The transformed positive monoclonal pBluescript plasmid was obtained using ampicillin selection. The plasmid was identified using EcoRI and XhoI double digestion to determine whether it contained target gene fragment and size of the inserted fragment.

DNA sequencing

After adding the universal forward primer T3 and reverse primer T7 to the pBluescript SK (+) plasmid containing the cDNA insert, forward and reverse sequencings of target gene were performed. The primer sequences were as follows: T3-20mer forward primer (100 µM), ATTAACCCTCACTAAAGGGA; and T7-20merreverse primer (100 µM), TAATACGACTCACTATAGGG. Homology analysis of the determined sequences to known genes was performed using the National Center for Biotechnology Information Basic Local Alignment Search Tool (NCBI BLAST) search tool in the GenBank database (http://blast.ncbi.nlm.nih.gov/Blast.cgi).

Function enrichment analysis

Gene Ontology (GO)37 is a powerful bioinformatics tool for annotating genes, and Kyoto Encyclopedia of Genes and Genomes (KEGG)38 is a set of databases for systematic analysis of biological pathway information. To study the pathogenesis of AS, the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/)39 was used to facilitate functional annotation. Using Metascape (http://metascape.org/gp/index.html)40 to perform functional enrichment analysis of target genes, we gained information on BPs, cytogenetics and signalling pathways. Using these two tools, a correlation between the function of screened genes and AS was found.

Gene–gene and protein–protein interactions

The BisoGenet and CytoNCA in Cytoscape were used to construct and analyse gene interactions. A bisogenet network can extract subnetwork topological features, based on nodes and connections to construct a gene screening network.41–43

GEO database (http://www.ncbi.nlm.nih.gov/geo) contains extensive data on disease-related gene expression. GEO database was used to determine whether identified genes were AS related. We searched GEO database and identified five series related to AS (GSE16561, GSE66724, GSE58294 GSE22255 and GSE230214). Peripheral blood samples of GSE16561 were from 39 patients with MRI-diagnosed ischaemic stroke and 25 controls who were non-stroke neurologically healthy.44 GSE66724 included eight stroke patients and eight controls, and microarray analysis of peripheral blood cells identified gene expression patterns associated with stroke.45 Blood samples from 69 patients with stroke and 23 controls were collected in GSE58294 and whole genome U133 Affymetrix arrays were performed.46 GSE22255 reported data on gene expression levels observed in peripheral blood mononuclear cells from 20 patients and 20 controls.47 The GSE230214 dataset reflected the genetic differences between patients with AS and healthy controls. After comparing gene expression levels of AS-related stroke patients and controls, the differentially expressed genes observed were used and significance level was set at 0.05.

The STRING V.11.0 Database (https://string-db.org/)48 was used to predict protein–protein associations. The minimum required interaction score was set to 0.4. The nodes and lines in the network graph represented target proteins and their interactions. The stronger interaction between two proteins, the thicker their connection.

Protein interaction function of Metascape was also used to synchronise protein–protein associations. The resulting network contains a subset of proteins that physically interact with at least one other member of the list.

The ConsensusPathDB database (http://consensuspathdb.org)49 integrates protein interactions, genetic interactions, metabolomes and drug-target interaction information. Interaction and function analyses of the screened genes were performed using ConsensusPathDB.

Biological General Repository for Interaction Datasets (BioGRID) (http://thebiogrid.org), formed by the aggregation of a great number of individual proteins or genetic interactions and interactions of RNA, DNA and other molecule metabolites, is fundamental to understanding gene–phenotype relationships and cellular functions.50

Expression analysis of PABPC1 in pan-cancer

BioGPS (http://biogps.org/%23goto=welcome) is an extensible and customisable gene annotation portal tool.51 We used it to identify the expression of newly discovered PABPC1 gene in various normal human organs and tissues.52

Gene Expression Profile Interaction Analysis (GEPIA) (http://gepia2.cancer-pku.cn/) can analyse gene expression in tumour and normal samples from The Cancer Genome Atlas (TCGA) and Genotype Tissue Expression databases.53 Pathological stages and ages of patients were selected to explore their association with PABPC1 expression. Patients were divided into pathological stages I, II, III and IV while 60 years old was set as a cut-off value to divide patients into two groups.

Human Protein Atlas (HPA) (https://www.proteinatlas.org/) (Protein Atlas V.21.0) is divided into three sections, those are, cell, tissue and pathology, which show expression of proteins in cells, normal tissues and cancerous tissues, respectively.54–56

TIMER V.2.0 (http://timer.cistrome.org/) uses six state-of-the-art algorithms to provide more reliable estimates of immune infiltration levels for tumour profiles provided by TCGA.57 58 TIMER V.2.0 offers immune associations, cancer exploration and immune estimation to study association among immune infiltration, genetic and clinical characteristics.

Kaplan-Meier plotter (http://kmplot.com/analysis/) assesses the impact of 54 000 genes (mRNA, miRNA and protein) on survival in pan-cancer59 and its main purpose is to discover and validate meta-analysis-based survival biomarkers. Overall survival (OS), disease-specific survival (DSS) and progression-free interval (PFI) were used to evaluate the relationship between PABPC1 expression and patient prognosis. Kaplan-Meier plotter was used to analyse the relationship between PABPC1 gene expression and survival rate of pan-cancer.60 61

UALCAN (http://ualcan.path.uab.edu) is a comprehensive and interactive web resource for analysing cancer OMICS data, which can provide access to publicly available cancer OMICS data. UALCAN was used to evaluate the expression of PABPC1 in normal and tumour samples.62

Tumour RNA-seq data (level 3) and clinical information were obtained from TCGA database (https://portal.gdc.com). Statistical analysis was performed using R software V.4.0.3 and immune correlations were assessed using immunedeconv.57 We used Assistant for Clinical Bioinformatics (https://www.aclbi.com) to determine the association between PABPC1 and pan-cancer immune correlations or tumour mutation burden (TMB).63

Statistical analyses

The relationship of PABPC1 expression level between tumour and normal tissue was described by boxplots using TIMER V.2.0, UALCAN and GEPIA V.2, which were calculated based on the Wilcoxon testing. Immune infiltration score and PABPC1 gene expression in multiple tumour tissues were analysed by Spearman correlation analysis. Significance was set as *p<0.05, **p<0.01 and ***p<0.001. OS, DSS and PFI were characterised using COX regression analysis with 95% CI for the HR. Survival curves were constructed using the Kaplan-Meier method. Survival analysis compared HR and log-rank p value of PABPC1 expression in pan-cancer. P value<0.05 was considered statistically significant and all tests were two tailed.

Selection and identification of immunoreactive clones from sera of patients with TIA

Serum samples from 23 patients (online supplemental table 1) were subjected to 3 rounds of screening. The plaques were transferred to NC membranes, in which negative clones showed pale plaque shadows, whereas the positive clones showed a clear purple-blue colour (figure 1A). One hundred and ninety-one mono-clones were screened after the first (figure 1A–a), second (figure 1A–b) and third (figure 1A–c) screening rounds. After internal shearing, 191 pBluescript SK (+) plasmids carrying SEREX antigen gene fragments were obtained and the concentration of all 186 plasmid DNA was >50 µg/mL. EcoRI and XhoI double digestion was used to identify 181 plasmid vectors with target cDNA fragments (figure 1B).

Request permission

A total of 181 plasmids were subjected to cDNA sequence detection using NCBI BLAST. Finally, 83 independent cDNA clones were obtained that showed high homology with known genes in RefSeq database (online supplemental table 2). Some of these 83 genes, such as BRAT1, WD repeat domain 1 (WDR1),26 matrix metalloproteinase 1 (MMP1), chromobox homolog 1 (CBX1), chromobox homolog 5 (CBX5),64 aldolase A (ALDOA) and fumarate hydratase (FH)65 had been reported by our research team.

TIA-screened genes consistent with AS-related functions

We conducted GO and KEGG pathway functional enrichment analyses using DAVID and found that 83 TIA-screened genes were mainly associated with components of some AS and tumour related biological pathways, such as focal adhesion, cell–substrate junction and platelet alpha granule in cellular components (CC), cadherin binding and actin binding in molecular function and viral process and negative regulation of response to external stimulus in BP (figure 2A). KEGG suggested two important pathways, GLYCOLYSIS/GLUCONEOGENESIS KEGG and PYRUVATE METABOLISM KEGG pathway maps, which were key signals of two chronic metabolic diseases, cancers and AS66–69 (figure 2B and C).

Request permission

In addition, Metascape demonstrated that 83 genes showed abundant gene function expression that was consistent with AS and cancers, such as haemostasis, glycolysis/gluconeogenesis, toxin transport and cellular response to growth factor stimulus (figure 2D). Top-level GO BPs, including localisation, metabolic processes and response to stimulus, were also viewed (figure 2E). Each node represented an enriched term and was coloured first by gene function expression cluster ID (figure 2F) and then by p value (figure 2G).

Therefore, we confirmed that the genes screened from sera of patients with AS-related TIA were consistent with AS and cancer related functions, making solid preparations for further research on identifying AS biomarkers.

PABPC1 expression examined through gene–gene and protein–protein interactions

After above studies confirmed that the genes screened were closely related to AS, it was necessary to further understand whether there were marker genes that could be targeted for AS treatment among these 83 genes to explore therapeutic mechanism of AS. Therefore, we used multiple gene–gene and protein–protein interaction methods to simultaneously screen for biomarkers.

First, 3 advanced screenings of 83 genes were performed by Bisogenet. The resulting hub network consisted of 196 proteins in total, which included 5 screening proteins, that were EEF1A1, HDAC2, PABPC1, VIM and IFI16 (figure 3A).

Request permission

Next, to simultaneously verify the reliability of 83 genes, we combined these genes as samples with 4 AS-related gene expression comprehensive datasets from GEO (GSE16561, GSE66724, GSE22255 and GSE58294), and the intersections between them represented coexistence genes. Among these, 64 genes in the centre existed in common positions in these 5 parts (figure 3B). The 64 genes were further processed by STRING and Cytoscape to present remaining genes with more centralised characteristics, including PABPC1, CBX1, MMP1, HDAC2, EEF1A1, VIM, ALDOA, IFI16 and CBX5 (figure 3C). To better use protein–protein interaction relationship to clarify marker genes of AS, we used Metascape and STRING again to respond to the corresponding proteins of these 83 genes (figure 3D and E).

Through intersection of above methods, we found that PABPC1, EEF1A1 and HDAC2 coexisted in various screening results and were likely to be marker genes of AS. Using ConsensusPathDB to conduct interaction and functional expression of these three genes, we found that their common intersection existed in RNA transport, metabolism of proteins, translation factors and translation (figure 3F). Combined with known research reports, we found that EEF1A1 and HDAC2 were common marker genes in cardiovascular diseases and various cancers,70–73 which further verified the reliability of our methods for screening AS and pan-cancer common biomarkers. Finally, we focused on the expression of PABPC1 and used BioGRID to produce an independent biological interaction network centred on it. We found that among the most related and adjacent genes around PABPC1, many genes had adverse effects on the occurrence and development of cancers, such as SNRNP70, LUC7L2, RPS3 and STAU1, which were circled in red (figure 3G). Therefore, we hypothesised that PABPC1 might play a role in the pan-cancer disease mechanism.

To confirm the reliable expression of PABPC1 in AS-related diseases, we performed reverse validation using GSE58294. In group comparison, there was a particularly significant statistical difference in the expression of PABPC1 between control group and stroke samples, ***p<0.001 (figure 4A). In the receiver operating characteristic curve, area under the curve was 0.771 and 95% CI was 0.666 to 0.875, p<0.001 (figure 4B). In addition, in order to further verify the difference of PABPC1 expression between AS and healthy people with external cohorts, the GSE230214 dataset was added and PABPC1 expression was indeed significantly increased in the AS group (figure 4C).

Request permission

The above results showed that we screened the AS disease marker gene PABPC1 using a variety of gene–gene and protein–protein interaction methods.

PABPC1 exhibited significant performance in pan-cancer

After identifying PABPC1 as an AS-related target gene, we were considered its expression in different contexts. First, we aimed to understand the mRNA expression of PABPC1 in various organs, tissues, and cell types of human body by BioGPS. BioGPS analysis showed that PABPC1 was expressed in a variety of organs and tissue cells (figure 5A). Among these, PABPC1 exhibited significantly higher expression levels in lymphoma burkitts. PABPC1 also exhibited a relatively high expression of clear aggregation in CD34⁺ and CD105⁺ endothelial cells, B lymphoblasts, CD19⁺ B cells, dendritic cells, CD8⁺ T cells, CD4⁺ T cells and so on.

Request permission

Because the inflammatory response, immune response and endothelial activation were all closely related to cancer, and these were also the responses shared by AS and cancer, the results of BioGPS suggested that we should further explore whether AS-related PABPC1 also had a predictive role in pan-cancer. Therefore, we applied various methods to explore the performance of PABPC1 in pan-cancer analysis. Differences in PABPC1 expression between all TCGA tumours and normal tissues were analysed and identified using the TIMER V.2.0 database (figure 5B). Relative to adjacent normal healthy tissue, PABPC1 was highly expressed in breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), oesophageal carcinoma (ESCA), etc.

Moreover, we applied UALCAN to perform simultaneous pan-cancer analysis of normal and tumour samples. We observed that the expression of PABPC1 in tumour tissues was significantly upregulated including breast cancer, colon cancer and ovarian cancer compared with normal tissues (figure 5C).

Simultaneously, we used GEPIA V.2 to generate PABPC1 expression profile of all tumour samples and paired normal tissues. We found that PABPC1 was significantly expressed in COAD, lymphoid neoplasm diffuse large b-cell lymphoma, ESCA, etc (shown in red font) (figure 5D).

These results suggested that we used numerous tumour databases and networks to determine the relationship between PABPC1 and pan-cancer, and high expression of PABPC1 in most cancers indicated a significant risk of carcinogenesis.

PABPC1 played a role in pan-cancer through tumour immune infiltration

Because existing studies had suggested a possible impact of cardiovascular disease on tumour immune infiltration, we aimed to further explore whether AS-associated PABPC1 was involved in cancer infiltration. TIMER V.2.0 was used to demonstrate the association of PABPC1 with various immune infiltrations in human tumours (figure 6A). It was significantly positively correlated with immune infiltration levels of several immune infiltrates, such as B cells, CD4⁺ T cells, neutrophils and mast cells. The positive correlation of PABPC1 in monocytes and myeloid dendritic cells was relatively low, and there was no special expression in CD8⁺ T cells, endothelial cells, macrophages, etc. However, abundance was negatively correlated with natural killer and plasmacytoid dendritic cells. It is noteworthy that tumours, including CHOL, kidney renal papillary cell carcinoma (KIRP) and kidney chromophobe (KICH), showed a positive correlation trend in almost every immune infiltration situation.

Request permission

Subsequently, to clarify whether there was a correlation between expression level of PABPC1 and TMB, which had an essential connection with immune checkpoint, we studied the role of PABPC1 in pan-cancer. The results showed that in prostate adenocarcinoma (PRAD), stomach adenocarcinoma (STAD) and adrenocortical carcinoma (ACC), etc, the expression of PABPC1 was significantly positively correlated with TMB (figure 6B). Correspondingly, COAD, ovarian serous cystadenocarcinoma (OV), uterine carcinosarcoma (UCS), etc, had significant negative correlation with it. figure 6C illustrated the correlation between PABPC1 expression levels and immune checkpoint expression levels, such as TIGIT. A significant negative correlation was found between PABPC1 and skin cutaneous melanoma (SKCM), Sarcoma (SARC) and so on, whereas a positive correlation was found with uveal melanoma (UVM), KIRP, etc.

Figure 7 showed the relationship between PABPC1 and tumour microenvironment across pan-cancer. PABPC1 was negatively associated with immune scores in ACC, bladder cancer (BLCA), SARC and so on, while it was positively correlated with BRCA, CHOL, glioblastoma (GBM), KIRP, etc.

Request permission

These findings suggested that AS disease-related PABPC1 played a role in cancer immune infiltration and affected TMB.

Association between PABPC1 expression and clinical features of patients

In the HPA, we observed significant differences in tissue expression of PABPC1 in several cancers (figure 8A). PABPC1 was found to express at low levels in glomerular cells and moderately expressed in renal tubular cells in normal kidney tissue. However, for renal cancer, the number of cytoplasmic/membranous staining was generally greater than 75%. For liver, PABPC1 was low expressed in normal tissues and moderately expressed in tumour tissues. For colon and testis conditions, PABPC1 was expressed at low levels in normal tissues but expressed at high in tumour ones.

Request permission

In addition, pathological stages and age of patients were selected to explore their association with PABPC1 expression. Patients were divided into four groups: I, II, III and IV of pathological stages (figure 8B). PABPC1 was significantly correlated with pathological stages in four cancer types, including COAD, KICH, KIRP and lung adenocarcinoma. In KIRP, PABPC1 expression increased as the pathological stage increased. In COAD and KICH, significant differences in PABPC1 expression were observed in stage IV.

Moreover, 60 years of age was set as a cut-off value to divide patients into two groups (figure 8C). PABPC1 was found to be highly expressed in patients more than 60 years old in BRCA, COAD, KIRP, liver hepatocellular carcinoma (LIHC), OV and SKCM, as well as highly expressed in less than 60 in PRAD.

The effect of different PABPC1 expression on the survival of patients with pan-cancer

Finally, knowing that PABPC1 had widespread effects on pan-cancer through immune infiltration and TMB, we studied its impact on the survival of cancer patients. Kaplan-Meier plotter was used to analyse the survival curves of patients with various types of cancer. We investigated the role of PABPC1 in pan-cancer prognosis using the COX model and found that PABPC1 played various prognostic roles in OS (figure 9A). PABPC1 acted as a protector in brain lower-grade glioma and Thymoma, and it was a detrimental prognostic factor of ACC, LIHC, KIRP, etc. For DSS, PABPC1 played a detrimental prognostic role in ACC, KIRP and SARC. In PFI, PABPC1 predicted poor outcomes in ACC, KICH, KIRP, etc.

Request permission

The survival curve showed that many tumours shared one characteristic, that was, compared with low expression of PABPC1, high expression made downward trend of survival curve steeper (figure 9B). HR was a ratio of 2 mortality risk rates in patients with high and low PABPC1 expression per unit time. All the PABPC1-expressing cancer survival curves above had HR between 1.23 and 5.75.

The above studies supported that high expression of PABPC1 increased the risk of death in patients with several different cancers and that the PABPC1 gene had a clear role in pan-cancer and AS.