This article is mentioned in:

Introduction

At present, 14 human peroxisomal biogenesis factor (PEX) genes that encode peroxin proteins that serve a role in various stages of peroxisome biogenesis have been identified, including peroxisome matrix protein input, membrane formation and peroxisome proliferation (1,2). Peroxisome membrane proteins are also thought to function as signaling platforms in reactive oxygen species (ROS)-induced autophagy and antiviral immunity (3,4). Peroxisome is a metabolic organelle involved in cellular redox balance and lipid metabolism (5,6). The function of peroxisome is important for ether phospholipid synthesis, fatty acid oxidation, bile acid synthesis and ROS homeostasis (7). Although peroxisomes are involved in cell metabolism, their functional effects in cancer remain unclear compared with those of other metabolic organelles such as mitochondria (8,9).

In the last decade, peroxisome has also been recognized as a central regulator of immunity (9). Lipid metabolites of peroxisomes, such as polyunsaturated fatty acids (PUFAs), are precursors of important immune mediators, including leukotrienes (LTs), and suppressors (10–12). Peroxisome redox metabolism regulates cellular immune signaling, such as activation of B-cell activated NF-κB (13). In addition, the development and activation of innate and adaptive immune cells are regulated by peroxisome β-oxidation and ether lipid synthesis (12,14). These findings open up avenues for targeting peroxisome interventions for immune disorders, inflammation and cancer. In addition, Lee et al (2) found that PEX13 is required for selective autophagy (virus autophagy) of Sindbis virus and damaged mitochondria (mitochondrial autophagy), and PEX13 mutant I326T and W313G associated with the disease are defective in mitochondria autophagy. The mitophagy function of PEX13 is shared with another peroxin family member, PEX3, but not with the other two peroxins required for general autophagy, namely PEX14 and PEX19 (15). Since PEX13 gene has a significant role in various stages of peroxisome biogenesis and autophagy, and peroxisome or autophagy has a crucial role in tumor occurrence and development, the present study aimed to investigate the biological functions that PEX13 may play in multiple cancers, especially in pancreatic adenocarcinoma (PAAD) (16,17).

In the present study, systematic bioinformatics analysis was conducted on patient data retrieved from The Cancer Genome Atlas (TCGA) database to verify the biological function and prognostic significance of PEX13 in various tumors. The present study comprehensively explored the functional significance of alterations in PEX13 expression level in various tumors in terms of PEX13 expression level, prognosis, genetic alterations, PEX13-related gene enrichment analysis, weighted correlation network analysis (WGCNA), protein interaction, long non-coding (lnc)RNA/circular (circ)RNA-micro (mi)RNA network and tumor immunity. As a general surgery team, the present authors are particularly interested in PAAD and focused. Finally, the relationships between the expression of PEX13 and the biological functions of PAAD were verified through in vitro experiments.

Materials and methods

PEX13 mRNA and protein expression and prognostic analyses

The PEX13 mRNA expression data for 33 tumors retrieved from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) database were analyzed using the SangerBox portal (http://SangerBox.com/Tool) (18). These tumors included adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large b-cell lymphoma (DLBC), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck cancer (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian cancer (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), sarcoma (SARC), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma (UCS) and uveal melanoma (UVM). PEX13 mRNA expression distributions were visualized in violin plots. Subsequently, the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) portal (genepattern.org) was applied to further explore the PEX13 mRNA expression in 33 tumors (19). The immunohistochemical (IHC) staining data for the analysis of PEX13 expression in various tumor tissues and corresponding normal tissues were obtained from The Human Protein Atlas (THPA) portal (https://www.proteinatlas.org/) (20). The analyzed tumors included glioma and breast, colorectal, liver, lung, pancreatic, skin, stomach, thyroid and renal cancer. Survival prognosis analysis was performed using GEPIA2 and SangerBox. Through GEPIA2, the overall survival (OS) map data for PEX13 in multiple tumors of TCGA database was obtained. The low- and high-expression cohorts of PEX13 were obtained using the expression thresholds of cutoff-high (>50%) and cutoff-low (<50%) values. Cox regression analysis was used to analyze the OS, disease-free interval (DFI), progression-free interval (PFI) and disease-specific survival (DSS) data in various cancers using SangerBox portal.

Genetic alteration analysis

The characteristics of PEX13 genetic alterations were explored in multiple tumors using the cBioPortal database (https://www.cbioportal.org/), including mutation, amplification and deep deletion (21). The structural variant data, mutation data and copy number alteration (CNA) data were obtained from the cBioPortal database. Data of somatic CNAs and somatic mutations of PAAD were obtained from TCGA datasets (gdc.cancer.gov/). The CNAs were related to the PEX13 expression and the threshold CNA peaks were analyzed through GISTIC 2.0 (https://cloud.genepattern.org/) (22). According to the PEX13 expression level, the patients with PAAD were divided into the first 25% PEX13low (n=46) and the last 25% PEX13high (n=46) groups. To download and visualize the somatic mutations of patients with PAAD, the ‘maftools’ package (bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html#7_Visualization) was applied in R software (4.2.1) (https://www.r-project.org/) (23).

PEX13-related genes enrichment analysis

Through the GEPIA2 portal, the top 200 PEX13-related genes were explored in TCGA database. Subsequently, the TIMER2 (http://timer.cistrome.org/) portal was used to generate a heatmap of the top 10 PEX13-related genes in various cancers, including KCMF1, KIAA1841, MPHOSPH10, MRPL19, PPP1CB, PUS10, RAB1A, SLC30A6, SMEK2 and YIPF4 (24). The correlation coefficient (R) and P-value were calculated. The Gene Ontology (GO)-biological process (BP) subontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 200 PEX13-correlated genes in pan-cancers were explored via DAVID (https://david.ncifcrf.gov/home.jsp) and SangerBox portal. Moreover, the top 50 PEX13 negatively related genes and top 50 PEX13 positively related genes were explored and heatmaps were generated in PAAD using LinkedOmics (www.linkedomics.org/login.php) (25). The GO-BP and KEGG enrichment analysis of the top 50 PEX13 negatively related genes and top 50 PEX13 positively related genes in PAAD were explored via DAVID and SangerBox portal. Linkedomics (linkedomics.org) portal was used to analyze the enrichment pathways of PEX13 positively and negatively related genes in PAAD (25).

WGCNA

The hub genes affected by PEX13 expression level were obtained through WGCNA using SangerBox portal. First, the PAAD gene expression profile in TCGA database was used to calculate the median absolute deviation (MAD) of each gene. The outlying genes and samples were removed using the GoodSamplesGenes method of the ‘WGCNA’ R software package (genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA) (26). WGCNA was then applied to build a scale-free co-expression network. As a soft-thresholding parameter, β could emphasize strong gene-gene associations and penalize weak correlations. To ensure a scale-free network, the power of β=7 (0.85) was selected as the soft-thresholding parameter and the power of β=7 (16.23) was selected as the soft-thresholding parameter in mean connectivity. Moreover, the PAAD samples were clustered using Pearson's correlation and average linkage method. The average linkage hierarchical clustering identified 23 modules and the differences among these were analyzed according to the expression level of PEX13, among which darkgrey module had the strongest association with PEX13 expression. Consequently, the darkgrey module was chosen as the important module for extracting hub genes.

Protein-protein interaction and ncRNA-miRNA-PEX13 regulatory networks analyses

Through the STRING database (https://string-db.org/), 50 PEX13 interacting proteins were identified in pan-cancers. Using Cytoscape (3.10.1; cytoscape.org/download.html), a PEX13 interacting protein network map was generated. The key protein-protein interaction network of PEX13 was extracted in this protein-protein interaction network using Cytoscape. Molecular docking technique was applied to predict the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3. By using the HDOCK server (http://hdock.phys.hust.edu.cn/), the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 were analyzed. Using the PyMOL tool (2.5.5; pymol.org/2/), the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 were visualized (27).

In addition, lncRNA/circRNA-miRNA-PEX13 networks were investigated. Through the niRabel (http://bioinfo.univ-rouen.fr/mirabel/index.php), miRDB (http://mirdb.org/), TargetScan (http://www.targetscan.org/vert_72/) and ENCORI (https://starbase.sysu.edu.cn/) databases, 16 miRNAs that might target PEX13 mRNA were identified. Subsequently, four miRNAs were further screened and selected through miRabel, miRDB and miRWalk scores, including hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p. Using the ENCORI online tool, 105 lncRNAs and 96 circRNAs that targeted hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p were predicted, and the lncRNA-miRNA-PEX13 and circRNA-miRNA-PEX13 regulatory networks were visualized via Cytoscape tool.

Immune-related analysis

The SangerBox portal was used to investigate the relationships between PEX13 mRNA expression and tumor infiltration immune cells (TIICs) in various cancers, including T_cells_CD8, T_cells_CD4_memory_resting, B_cell_memory, Tregs, NK_cells_activated, macrophages and dendritic_cells_activated. Subsequently, the correlations between PEX13 mRNA expression and immune checkpoint (ICP) genes, microsatellite instability (MSI), tumor mutation burden (TMB) and tumor purity in multiple TCGA tumors we investigated using SangerBox. Furthermore, TIMER2 portal was applied to analyze the associations between PEX13 and the infiltration levels of B cell and macrophage in pan-cancers, especially in PAAD via different algorithms, such as XCELL, CIBERSORT and EPIC. Moreover, the relationships between PEX13 and stromal score, immune score and ESTIMATE score in multiple cancers were explored using SangerBox, especially in lung adenocarcinoma (LUAD), stomach and esophageal carcinoma (STES), sarcoma (SARC) and lung squamous cell carcinoma (LUSC).

Cell culture

The human PAAD cell lines of PANC1 and BxPC3 were obtained from the Shanghai Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Comprehensive cell line authentication was performed and cells were periodically checked for mycoplasma contamination. The DMEM (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (FBS; Satorius AG) was used to culture the PAAD cell lines. The PAAD cells were incubated in a cell incubator at 37°C with 5% CO2.

Small-interfering (si)RNAs delivery, RNA isolation and reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

siRNAs used to knock down the expression of PEX13 were obtained from Shanghai GenePharma Co., Ltd. A total of 2×105 PAAD cells were seeded in 6-well plates, 4.5 µl Lipofectamine™ RNAiMAX reagent (Thermo Fisher Scientific, Inc.) and 40 pmol siRNA were mixed for 10 min in each well at room temperature. Subsequently, the mixture was added to the cells cultured with FBS-free medium and the cells were replaced with full medium at 24 h after transfection. The PAAD cells were collected after 48 h with TRIzol® reagent (Takara Bio, Inc.). The GoScript reverse transcription system (Promega Corporation) was applied to generate PEX13 cDNA and the GoTaq® qPCR Master Mix (Promega Corporation) was used to detect the knockdown efficiency of PEX13 siRNAs through qPCR using the ABI QuantStudio 3 system. The siRNA sequences were as follows: si-PEX13-1 sense, 5′-ACUUGAUUCCACCGUUUUCCUTT-3′ and antisense, 5′-GAAAACGGUGGAAUCAAGUAATT-3′; si-PEX13-2 sense, 5′-UUAUUCUAGUCCUGAAUUCCUTT-3′ and antisense, 5′-GAAUUCAGGACUAGAAUAAAGTT-3′; and siRNA negative control (si-NC) sense, 5′-UUCUCCGAACGUGUCACGUTT-3′ and antisense, 5′-ACGUGACACGUUCGGAGAATT-3′. The thermocycling conditions used for qPCR were: Initial denaturation at 95°C for 10 min; 40 cycles of amplification for 15 sec at 95°C and 45 sec at 60°C, and an extension step at 72°C for 2 min. The relative expression levels of the PEX13 mRNA were calculated using the 2−ΔΔCq method and normalized to the internal control GAPDH (28). The primer sequences were as follows: GAPDH forward, 5′-GGTGGTCTCCTCTGACTTCAACA-3′ and reverse, 5′-GTTGCTGTAGCCAAATTCGTTGT-3′ reverse; PEX13 forward, 5′-GCCCCACTTTCCAATCTGCT-3′ and reverse, 5′-AGAATAGGTGGGGGCACTCT-3′.

Cell Counting Kit-8 (CCK-8) assay

A total of 2×103 PANC1 or BxPC3 cells were plated into 96-well plates. PEX13 expression was knocked down using the PEX13 siRNAs. Subsequently, the optical density (OD) at 450 wavelengths was measured for 5 consecutive days (at days 0, 1, 2, 3, 4 and 5) after incubation for 1 h at 37°C with the CCK-8 reagent (ApexBio Technology) using a multiplate reader.

Wound healing and Transwell assays

A total of 2.5×105 PANC1 or BxPC3 cells were plated in 6-well plates. A 200 µl pipette tip was used to scratch the cell monolayer. After 0, 12 and 24 h incubation in an FBS-free medium, the images of wound closure were captured and further explored via an inverted microscope (Olympus Corporation) and Image J software (1.8.0; National Institutes of Health). For the Transwell assay, 8-µm pore filters coated with Matrigel matrix glue (Corning, Inc.) were used. Matrigel was added into the upper chamber for 40 min at 37°C. A total of 5×104 PANC1 and BxPC3 cells in 300 µl serum-free DMEM were seeded into the upper chamber and 600 µl DMEM with 10% FBS medium was added to the lower chamber. After 24 h, cells were fixed with 4% paraformaldehyde at room temperature and stained with 1% crystal violet at room temperature. The images were captured using a microscope (BX53; Olympus Corporation).

Colony formation assay

In 6-well plates, 1×103 PANC1 and BxPC3 cells were seeded in each well with 2.5 ml DMEM containing 10% FBS. The cells were incubated for 2 weeks at 37°C in a humidified chamber with 5% CO2. The colonies were fixed with 4% paraformaldehyde for 30 min at room temperature before being stained with crystal violet for 10 min at room temperature. A cell mass of >50 cells was considered a colony, and Image J (1.8.0) was used for counting colonies.

Statistical analysis

All of the experimental data were statistically analyzed using GraphPad Prism 8.0 (Dotmatics). Comparisons between two groups were made using an unpaired t-test and the statistical method used for survival analysis was Kaplan-Meier. Cox regression analysis was performed using logrank test to obtain prognostic significance. The correlation analyses in the manuscript were applied using Pearson's correlation analysis. One-way ANOVA with Bonferroni's post-hoc test was used to compare the differences between multiple groups. Two-way ANOVA with Bonferroni's post-hoc test was used to compare the differences between CCK-8 groups (si-NC and si-PEX13-1/si-PEX13-2). The data are presented as mean ± standard deviation from three individual experiments and each experiment was repeated at least three times. P<0.05 was considered to indicate a statistically significant difference.

Results

PEX13 expression patterns in various cancers

Through TCGA and GTEx databases analysis, PEX13 was revealed to be highly expressed in multiple tumors, including GBM, LGG, GBMLGG, BRCA, UCEC, LUAD, CESC, ESCA, COAD, STES, READ, STAD, PRAD, HNSC, LUSC, LIHC, high-risk wilms tumor (WT), BLCA, THCA, PAAD, UCS, LAML, acute lymphoblastic leukemia (ALL) and ACC. In addition, the expression level of PEX13 was downregulated in KIRP, pan-kidney cohort (KIPAN), KIRC, SKCM, TGCT and KICH compared with the corresponding normal tissues (Fig. 1A). In GEPIA2 portal, it was discovered that PEX13 was remarkable upregulated in DLBC, GBM, LGG, LUSC, PAAD, STAD and THYM (Fig. 1B). In addition, the results of IHC showed that the expression of PEX13 protein in breast cancer, colorectal cancer, glioma, liver cancer, lung cancer, pancreatic cancer, skin cancer, stomach cancer, thyroid cancer and renal cancer was higher than that in normal tissues (Fig. S1). The present findings indicated that PEX13 mRNA and protein expression levels are increased in a variety of tumors, suggesting that PEX13 may function as an oncogene in various cancers, including PAAD.

Figure 1. PEX13 expression levels in normal tissue and tumors. (A) Comparison of PEX13 expression level among 33 tumor tissues and normal tissues via The Cancer Genome Atlas and Genotype-Tissue Expression databases. (B) The differences of PEX13 expression levels in tumor and normal samples by GEPIA2. **P<0.01, ***P<0.001 and ****P<0.0001. PEX, peroxisomal biogenesis factor.

PEX13 expression and prognosis across multiple types of cancer

The prognostic value of PEX13 in various tumors was analyzed by comparing the OS of patients with PEX13 expression. Analysis of the pan-cancer cohort using the GEPIA2 database demonstrated that high PEX13 mRNA expression levels was associated with a significantly lower OS in patients with ACC, KICH, LGG, LIHC and PAAD, and high PEX13 mRNA expression levels were associated with a significantly higher OS in patients with KIRC (Fig. 2A). Cox regression analysis in SangerBox database suggested that PEX13 mRNA expression level was related to OS, DFI, DSS and PFI in patients with various tumors. Regarding the OS, the findings revealed that high PEX13 mRNA expression was linked to a shorter OS in GBMLGG, LGG, PAAD, ACC, KICH, BRCA, LIHC and LAML (Fig. 2B). Regarding the DFI, high PEX13 mRNA expression was linked to a shorter DFI in PAAD, KIPAN, KIRP and ACC (Fig. S2A). In PFI, the high PEX13 expression was related to shorter PFI in GNMLGG, ACC, LGG, PAAD, LIHC, KICH and UVM (Fig. S2B). Furthermore, a high PEX13 mRNA expression was correlated with a shorter DSS in GBMLGG, LGG, PAAD, ACC, KICH and LIHC (Fig. S2C). It was also observed that high PEX13 mRNA expression was significantly correlated with shorter OS, DFI, PFI, and DSS in PAAD and ACC. The present results revealed that the level of PEX13 mRNA expression was related to the prognosis of various tumors. Furthermore, high PEX13 mRNA expression was linked to a worse prognosis in patients with PAAD.

Figure 2. Prognosis analysis of PEX13 in pan-cancers. (A) Survival maps of overall survival according to PEX13 expression in various cancers via Gene Expression Profiling Interactive Analysis 2. (B) Cox regression analysis according to PEX13 expression in multiple tumors. PEX, peroxisomal biogenesis factor.

PEX13 genetic alterations in various cancers

Mutations, deletions or amplification of oncogenes or tumor suppressor genes have been linked to the growth and progression of several tumors (29). Hence, the present study analyzed diverse types of genetic alterations in the PEX13 gene using the cBioPortal portal, including amplification, mutation and deep deletion. The most common genetic alteration of PEX13 gene in DLBC (18.75%), LUSC (5.75%), BLCA (2.68%), UCS (1.75%), ovarian serous cystadenocarcinoma (1.54%) and ESCA (1.1%) was amplification. In UCEC (2.46%), STAD (1.36%), COADREAD (1.35%) and KIRP (0.71%), the most common genetic alteration of PEX13 gene was mutation (Fig. 3A). Subsequently, the TCGA-PAAD dataset was interrogated to investigate the correlation between PEX13 and specific genomic features such as copy number variations (CNVs) and somatic mutations. Fig. 3B shows the comparison of the CNV profiles in the PEX13high (n=46) and PEX13low (n=43) groups. In the PEX13high group, deletion peaks were found in the 5p12.3, 6p25.3, 9p21.3, 10p15.3, 10q23.2, 18q12.1, 18q21.2 and 19p13.2 chromosomal locations, while amplification peaks were found in the 1p12, 6p21.31, 7q21.13, 8q24.21, 9p13.3, 12p13.33, 18q11.2 and 19q13.2 chromosomal regions (Fig. 3C-D). In the PEX13low group, deletion peaks were found in the 9p23 and 9p21.3 chromosomal locations, while frequent amplification peaks were found in the 8p11.23 and 19q13.2 chromosomal regions (Fig. 3C and D). These findings suggested that high PEX13 expression led to more deletions and amplifications of multiple sites in the genome of patients with PAAD. The PEX13low group indicated that high frequency of somatic mutations in the FRG1B (72%), KRAS (49%), TP53 (47%), UBBP4 (44%) and TTN (40%) genes and the PEX13high group suggested that high frequency of mutations in the KRAS (89%), FRG1B (76%), TP53 (72%), MAMLD1 (35%), MED12L (35%) and UBBP4 (30%) genes (Fig. 3E and F). These findings revealed that high PEX13 expression may lead to the increased mutation probability of KRAS, FRG1B, TP53 and other genes, and to the decreased mutation probability of UBBP4 gene in PAAD, which may be a potential mechanism leading to the occurrence and development of tumors.

Figure 3. Distinct genomic profiles associated with PEX13 expression in pan-cancers, including PAAD. (A) Genetic alterations of PEX13 in different cancers in TCGA database were discovered using the cBioPortal tool, including mutation, structural variant, amplification and deep deletion. (B) The copy number alteration profile analysis between 25% PEX13low and 25% PEX13high groups in a PAAD cohort was preformed using GISTIC2.0. Frequencies of (C) deletion and (D) amplification of PAAD associated with PEX13 expression (blue, deletion; red, amplification). Detection of somatic mutations in PAAD, including (E) 25% PEX13low and (F) 25% PEX13high group. PEX, peroxisomal biogenesis factor. TCGA, The Cancer Genome Atlas; PAAD, pancreatic adenocarcinoma.

These findings suggested the presence of PEX13 gene mutations, amplification and deletion in multiple tumors. The PAAD tissues revealed different somatic mutations and CNVs depending on the level of PEX13 expression. This suggested that alterations in the PEX13 gene could influence the occurrence and progression of multiple cancers, particularly PAAD.

Pathway enrichment analysis of PEX13-associated genes

The PEX13-related genes were analyzed in 33 cancers through GEPIA2 database and 200 PEX13-related genes were obtained. The top 10 genes with the highest correlation among 200 PEX13-related genes were generated correlation heat maps in 33 cancers via TIMER2, including KCMF1, KIAA1841, MPHOSPH10, MRPL19, PPP1CB, PUS10, RAB1A, SLC30A6, SMEK2 and YIPF4 (Fig. 4A). Fig. 4B shows that PEX13 expression was positively correlated with PUS10 (R=0.61), KIAA1841 (R=0.61), RAB1A (R=0.59), SMEK2 (R=0.58), YIPF4 (R=0.57), PPP1CB (R=0.57), KCMF1 (R=0.57), SLC30A6 (R=0.57), MPHOSPH10 (R=0.56) and MRPL19 (R=0.56) in GEPIA2. The GO-BP analysis of 200 PEX13-related genes in 33 tumors showed that PEX13-related genes were mainly enriched in mRNA metabolic process, cellular component biogenesis, cell division, intracellular transport, cytokinesis, mRNA processing, cell cycle and RNA splicing (Fig. 4C). PEX13-related genes were mainly enriched in Endocytosis, mRNA surveillance pathway, nucleocytoplasmic transport, T cell receptor, signaling pathway and hippo signaling pathway based on the results of the KEGG enrichment analysis (Fig. 4D).

Figure 4. The enrichment analysis of PEX13-related genes in pan-cancers. (A) Correlation heat map displayed the top 10 PEX13-related genes in various tumors via the TIMER2.0 portal. (B) Correlations between PEX13 expression and top 10 PEX13-related genes in pan-cancers via GEPIA2. The (C) Gene Ontology-Biological Process and (D) KEGG enrichment analyses of 200 PEX13-related genes in various cancers. PEX, peroxisomal biogenesis factor; TPM, transcripts per million; KEGG, Kyoto Encyclopedia of Genes and Genomes; UVM, uveal melanoma; UCS, uterine carcinosarcoma; UCEC, uterine corpus endometrial carcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; THYM, thymoma; STAD, stomach adenocarcinoma; SKCM, skin cutaneous melanoma; SARC, sarcoma; READ, rectum adenocarcinoma; PRAD, prostate adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PAAD, pancreatic adenocarcinoma; OV, ovarian serous cystadenocarcinoma; MESO, mesothelioma; LUSC, lung squamous cell carcinoma; LUAD, lung adenocarcinoma; LIHC, liver hepatocellular carcinoma; LGG, brain lower grade glioma; KIRC, kidney renal clear cell carcinoma; KICH, kidney chromophobe; HNSC, head and neck cancer; HPV, human papillomavirus; GBM, glioblastoma multiforme; ESCA, esophageal carcinoma; DLBC, lymphoid neoplasm diffuse large b-cell lymphoma; COAD, colon adenocarcinoma; CHOL, cholangiocarcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; BRCA, breast invasive carcinoma; BLCA, bladder urothelial carcinoma; ACC, adrenocortical carcinoma.

To further explore the biological functions of PEX13 in PAAD, the enrichment of PEX13 positively and negatively related genes and their enrichment pathways in PAAD were analyzed using the Linkedomics portal. Fig. S3A shows the heat maps of 50 positively and 50 negatively related genes of PEX13 in PAAD. Fig. S3B-E showed the GO-BP and KEGG enrichment analysis of 100 positively and 100 negatively related genes of PEX13, respectively. The results showed that GO-BP was mainly enriched in positive regulation of drug response, small GTPase mediated signal transduction, regulation of epidermal growth factor-activated receptor activity and ERBB signaling pathway among the positively related genes of PEX13 and in the regulation of glycogen (starch) synthase activity, intracellular protein transport, intracellular protein transport and autophagosome maturation among the negatively related genes of PEX13 (Fig. S3B and C). For the KEGG enrichment analysis, PEX13-positive related genes were mainly enriched in Rap1 signaling pathway, ErbB signaling pathway, EGFR tyrosine kinase inhibitor resistance, Neurotrophin signaling pathway and AMPK signaling pathway and PEX13-negatively related genes were mainly enriched in Vitamin B6 metabolism, Apelin and Insulin signaling pathway (Fig. S3D-E). These findings show that PEX13 could regulate mRNA processing, RNA splicing, cell cycle and the response to drugs in cancers. In PAAD, PEX13 expression was related to protein transport, tumor cell response to drugs, multiple cancer pathways and immune pathways.

Analysis of WGCNA in PAAD

The TCGA-PAAD cohort was analyzed via WGCNA to identify hub genes strongly associated with PEX13 expression in PAAD. Using the average linkage approach and Pearson's correlation method, the samples of TCGA-PAAD data were first clustered (Fig. 5A-C). With the help of the average linkage hierarchical clustering, a total of 23 modules were obtained (Fig. 5D). According to the expression level of PEX13, we analyzed the differences among 23 modules, among which darkgrey module had the strongest association with PEX13 expression (Fig. 5E-F). The darkgrey module was selected as the significant module for further analysis. Fig. 5G shows the hub genes of darkgrey module and GO-BP and KEGG enrichment analyses were used to analyze these hub genes. In GO-BP analysis, hub genes mainly enriched in cell migration, regulation of ERK1 and ERK2 cascade, wound healing, ephrin receptor signaling pathway, ERBB2-EGFR signaling pathway, positive regulation of MAP kinase activity and signal transduction (Fig. 5H). In KEGG enrichment analysis, hub genes mainly enriched in PI3K-Akt signaling pathway, extracellular matrix (ECM)-receptor interaction, focal adhesion, pathways in cancer and MAPK signaling pathway (Fig. 5I). These results suggested that the hub genes of darkgrey module with different PEX13 expression were mainly involved in cell migration, wound healing, PI3K-Akt and MAPK signaling pathway.

Figure 5. Weighted gene co-expression network analysis of PAAD associated with PEX13 expression. (A) Analysis of the scale-free fit index for various soft-thresholding powers (β). (B) Analysis of the mean connectivity for various soft-thresholding powers. (C) Clustering dendrogram of PAAD samples. (D) Module feature vector clustering. (E) Dendrogram of all differentially expressed genes clustered based on a dissimilarity measure. (F) Heatmap of the correlation between module eigengenes and PEX13 expression level of PAAD. (G) The hub genes diagram of darkgrey module. The GO-BP (H) and KEGG (I) enrichment analyses of hub genes of darkgrey module in PAAD. GO-BP, Gene Ontology-Biological process; KEGG, Kyoto Encyclopedia of Genes and Genomes; PEX, peroxisomal biogenesis factor; PAAD, pancreatic adenocarcinoma.

Construction of protein interaction and ncRNA-miRNA-PEX13 regulatory network

The 50 PEX13 interacting proteins were identified using STRING database and the PEX13 interacting protein network map was generated in Cytoscape (Fig. 6A). In this protein interaction network, key protein interaction network of PEX13 was extracted using Cytoscape, as shown in Fig. 6B. The PEX13 interacting proteins were scored in Cytoscape and a larger volume indicated higher interacting score of these proteins in the network diagram. Subsequently, four proteins had the highest interaction scores with PEX13, including PEX2, PEX12, PEX14 and ABCD3 (Fig. 6B). Molecular docking technique was applied to forecast the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 (Fig. 6C-F). The results showed that PEX13 may interact with PEX2, PEX12, PEX14 and ABCD3 to influence the biological functions and acted as a significant role in various tumors.

Figure 6. Protein-protein interaction analysis of PEX13. (A) Protein-protein interaction network of PEX13 by STRING portal. (B) Key protein-protein interaction network of PEX13 using STRING and Cytoscape tools. The binding sites between PEX13 and PEX2 (C), PEX12 (D), PEX14 (E) and ABCD3 (F) via molecular docking technique. PEX, peroxisomal biogenesis factor; ABCD3, TP binding cassette subfamily d member 3.

In addition, numerous studies found that lncRNA and circRNA regulated the expression of downstream mRNA by adsorbing target miRNAs, which acted as a crucial factor in tumorigenesis and the progression of tumors (30,31). Therefore, the present study explored the lncRNA/circRNA-miRNA network that may regulate mRNA expression in multiple cancers. LncRNA-miRNA regulatory networks of PEX13 were constructed to investigate the underlying molecular mechanisms of PEX13 in tumors. Through the miRabel, miRDB, TargetScan and ENCORI databases, 16 miRNAs that might target PEX13 mRNA were identified (Fig. 7A). Subsequently, four miRNAs, including hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p, were screened and selected using miRabel, miRDB and miRWalk. Using the ENCORI online tool, the lncRNAs and circRNAs that targeted their corresponding miRNAs were predicted and Cytoscape was used to visualize the lncRNA-miRNA and circRNA-miRNA regulatory networks (Fig. 7B and C). These findings suggested the upstream ncRNA-miRNA regulatory networks that may control the aberrant expression of PEX13 in cancers.

Figure 7. Construction of lncRNA/circRNA/miRNA-PEX13 network. (A) Venn diagram showed the 16 miRNAs that might target PEX13 mRNA by the miRabel, miRDB, TargetScan and ENCORI databases and 4 miRNAs were eventually identified via miRabel, miRDB and miRWalk databases. (B) Construction of lncRNA-miRNA-PEX13 network through Cytoscape tool. (C) Construction of circRNA-miRNA-PEX13 network through Cytoscape tool. LncRNA, long non-coding RNA; miRNA, microRNA; circRNA, circular RNA; PEX, peroxisomal biogenesis factor.

PEX13 regulated tumor immunity in multiple tumors

The relationships between PEX13 and the tumor immune microenvironment (TIM) were further analyzed in a variety of tumors, and the correlation between PEX13 and the composition of tumor infiltrating immune cells (TIICs) was investigated in multiple cancers. Studies showed that TIICs were a significant component of the tumor microenvironment and regulated the occurrence, development and metastasis of tumors (32,33). The Sangerbox database was applied to investigate the correlation between PEX13 and TIIC levels in diverse tumors. The expression of PEX13 in LAML, PRAD, BRCA, PAAD, STES, HNSC, LIHC, KIRC, THCA, STAD, SKCM, THYM and CESC was remarkably related to various immune cells, including T_cells_CD4_memory_resting, B_cell_memory, T_cells_CD8, macrophages_M1, macrophages_M2, Tregs, NK_cells_activated, dendritic_cells_activated and neutrophils (Fig. 8A). Fig. S4A demonstrated the relationship between PEX13 expression and B cell infiltration levels in different tumors. In PAAD, the expression level of PEX13 under different algorithms was mainly negatively related to the infiltration level of B cells (Fig. S4B). Fig. S4C presents a heat map of the relationship between PEX13 expression levels and macrophage infiltration levels in various tumors. PEX13 expression level was also negatively related to the infiltration level of macrophages under different algorithms in PAAD (Fig. S4D).

Figure 8. Correlations between PEX13 expression and immunity in pan-cancers. (A) Association between PEX13 expression level and immune cell infiltration level in pan-cancers via Sangerbox. (B) The correlations between PEX13 expression and immune checkpoint genes in multiple tumors. The correlations between PEX13 expression and microsatellite instability (C), tumor mutation burden (D) and tumor purity (E) in various tumors. *P<0.05. PEX, peroxisomal biogenesis factor; MSI, microsatellite instability; TMB, tumor mutation burden; ICP, immune checkpoint.

Furthermore, the correlation between PEX13 expression level and ICP genes, MSI, TMB and tumor purity were analyzed via SangerBox. The occurrence of MSI is caused by the functional defect of DNA mismatch repair in tumor tissues (34). The presence of MSI and ineffective DNA mismatch repair is a crucial clinical tumor marker (35). TMB has a strong relationship to the effectiveness of PD-1/PD-L1 inhibitors and can forecast their benefits in terms of efficacy (36). The present study discovered that PEX13 expression was positively linked with ICP genes in several tumor types (Fig. 8B). PEX13 expression was positively correlated with MSI in TGCT and KIRC and negatively correlated with MSI in GBMLGG, HNSC, PRAD, KIPAN, COAD, DLBC and THCA (Fig. 8C). In LUAD and BRCA and pheochromocytoma and paraganglioma, PEX13 expression was linked favorably with TMB (Fig. 8D). In tumor purity, PEX13 expression was negatively related to purity in KICH, KIPAN and PRAD and positively related to purity in GBMLGG, SARC, LUAD, LUSC, STES, LGG, STAD, BRCA, ESCA, etc. (Fig. 8E). Moreover, it was evaluated the relationship between the ESTIMATE score and PEX13 expression level in various cancers. ImmuneScore reflects the proportion of infiltrating immune cells in cancer tissues. The proportion of stromal cells in tumor tissue was reflected by the Stromal score (StromalScore). The ESTIMATE score (ESTIMATEScore) is the sum of immune and stromal scores and reflects the state of the TIM and tumor purity (18,37). The present study found significant relationships between PEX13 levels and StromalScore, ImmuneScore and ESTIMATEScore across multiple tumors, especially in LUAD, STES, SARC and LUSC (Fig. S5A-C). These findings indicated that PEX13 expression regulated the sensitivity and resistance of certain cancers to immunotherapy. Thus, PEX13 emerged from this analysis as a potential immunotherapy biomarker and predictor of tumor immunotherapy response.

PEX13 promoted cell proliferation, migration, invasion and colony formation in PAAD

PANC1 and BxPC3 cell lines were used to verify the biological functions of PEX13. The knockdown of PEX13 expression in PANC1 and BxPC3 cells transfected with PEX13 siRNAs was confirmed using RT-qPCR (Fig. 9A). The CCK-8 assay suggested that PEX13 knockdown inhibited the cell proliferation of PANC1 and BxPC3 cells in vitro (Fig. 9B). In addition, wound healing and Transwell experiments demonstrated that PEX13 knockdown suppressed the cell migration and invasion of PANC1 and BxPC3 cells in vitro (Fig. 9C-D). The colony formation assay suggested that PEX13 knockdown could inhibit the cell colony formation ability of PANC1 and BxPC3 cells in vitro (Fig. 9E). These findings showed that PEX13 could increase cell proliferation, migration, invasion and colony formation abilities of PANC1 and BxPC3 cell lines in vitro.

Figure 9. Biological functions of PEX13 in PAAD. (A) Reverse transcription-quantitative PCR assay verified the knockdown efficiency of PEX13 small interfering RNAs. (B) Cell Counting Kit-8 assay measured the cell proliferation after PEX13 knockdown in PANC1 and BxPC3 cells. (C) Wound healing, (D) Transwell and (E) colony formation assays measured cell migration, invasion and colony formation abilities after PEX13 knockdown in PANC1 and BxPC3 cells. *P<0.05, **P<0.01, ***P<0.001 and ****P<0.0001. si-, small interfering RNA; NC, negative control; PEX, peroxisomal biogenesis factor; si-PEX13, siRNA targeting PEX13.

Discussion

Metabolic reprogramming is one of the hallmarks of cancer and has a selective advantage in the occurrence and progression of cancer cells (38). Metabolites ingested through cell membranes are mostly initially processed in the cytoplasm and then typically transported to intracellular organelles for further processing, such as peroxisome and mitochondria (39). Although mitochondria are key metabolic hubs of cancer cells, the metabolic functions of other organelles in cancer cells are rarely studied and need further exploration.

In mammals, the peroxisome, a single membrane organelle, participates in >50 separate enzymatic functions (40). Defects in genes encoding peroxisome proteins are associated with peroxisome disease. Lipid metabolism is also critical to tumorigenicity (8). Lipid synthesis pathways (such as the production of lipogenic triglycerides and the synthesis of phospholipids and cholesterol) are critical for energy storage, cell membrane structure and cell signaling regulation. For lipid synthesis and activation of the corresponding signaling network to be coordinated in rapidly growing cancer cells, a significant volume of cell membrane formation is required (41). High levels of ether phospholipid in certain cancers suggest that elevated lipid synthesis of peroxisome is related to tumor progression (42). PEX13 is an indispensable key molecule in the occurrence of peroxisome, but, to the best of our knowledge, its functions and effects in tumors have not been reported in the literature. Therefore, the present study aimed to investigate the biological functions and clinical significance of PEX13 in various tumors, especially in PAAD.

The present study first revealed that PEX13 was upregulated in several types of tumors compared with the corresponding normal tissues, including GBM, LGG, GBMLGG, UCEC, CESC, BRCA, ESCA, LUAD, STES, COADREAD, COAD, PRAD, STAD, LUSC, HNSC, WT, LIHC, THCA, BLCA, PAAD, UCS, ALL, LAML and ACC. In ACC, KICH, LGG, LIHC and PAAD, patients with high PEX13 mRNA expression had significantly lower OS, and patients with high PEX13 mRNA expression had significantly higher OS in KIRC. Cox regression analysis indicated that PEX13 mRNA was related to OS, DFI, PFI and DSS in various tumors. In addition, high expression of PEX13 was associated with poorer prognosis in patients with PAAD. Subsequently, genetic alteration analysis of PEX13 showed the presence of PEX13 gene mutations, amplification and deletion in multiple tumors. The PAAD tissues displayed different somatic mutations and CNVs depending on the level of PEX13 expression, indicating that the PEX13 gene alteration may mediate the occurrence and progression of multiple tumors, especially in PAAD. The GO-BP and KEGG enrichment analyses of PEX13-related genes found that PEX13 expression was significantly correlated with endocytosis, mRNA surveillance pathway, nucleocytoplasmic transport, T cell receptor signaling pathway and hippo signaling pathway in cancers. In PAAD, GO-BP and KEGG enrichment analyses of PEX13-related genes revealed that PEX13 expression was significantly associated with Rap1, ErbB, Neurotrophin, AMPK, Apelin and Insulin signaling pathways. WGCNA analysis suggested that PEX13-related hub genes were mainly enriched in PI3K-Akt signaling pathway, ECM-receptor interaction, focal adhesion and MAPK signaling pathway via KEGG enrichment analysis. Moreover, the present study discovered that PEX13 may interact with PEX2, PEX12, PEX14 and ABCD3 to influence the biological functions and serve as a significant factor in pan-cancers and the upstream lncRNA/circRNA-miRNA regulatory networks may regulate the aberrant expression of PEX13 in various cancers.

Furthermore, PEX13-related genes were mainly enriched in immune and cancer-related pathways, and the correlation between PEX13 and tumor immunity was further analyzed. As an important component of the immune microenvironment, TIICs act as a key factor in the occurrence, progression and treatment of tumors. PEX13 mRNA levels were associated with TIICs in multiple tumors, including PAAD. The present study evaluated the effect of PEX13 on immunotherapy sensitivity and resistance in patients with cancer. The ICP, MSI, TMB and tumor purity analyses suggested that PEX13 may be a promising target for the treatment of patients with certain tumors, especially in immunotherapy. In vitro assays also demonstrated that PEX13 knockdown could significantly inhibit the proliferation, migration, invasion and colony formation of PAAD cells.

In conclusion, the present results demonstrated that PEX13 may be an invaluable prognostic indicator and a possible predictor of immunotherapy sensitivity and resistance in patients with malignant tumors, especially in PAAD. However, the present study has limitations because it only uses TCGA database, while further databases should be used for corroboration, and only verifies the relationship between PEX13 and PAAD. The possible mechanisms by which PEX13 functions in PAAD were not further explored.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the Natural Science Foundation of Shanxi Province (grant no. 202103021223011).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

JS designed the study. PD and XD performed the bioinformatics analyses and experiments. TY and DL analyzed the results. YW and YD performed data analysis and wrote the manuscript. JS and PD confirm the authenticity of all the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References