PC gene expression patient cohorts utilized for this study included (i) the Taylor et al. (7) cohort, used for model training/discovery; (ii) The Cancer Genome Atlas Prostate Cancer Adenocarcinoma (TCGA PRAD) cohort (6), used for test/validation; and (iii) the Deutsches Krebsforschungszentrum (DKFZ, The German Cancer Research Center) PRAD cohort (8), used for additional test/validation ( Table 1 ).

The Taylor et al. cohort (7) is an MSKCC dataset composed of primary tumor prostatectomy samples (n = 138), profiled on the Affymetrix Human Exon 1.0 ST Array, and obtained from the Gene Expression Omnibus (GEO) GSE21032. Gene expression data were normalized using the robust multi-array average (RMA) algorithm for background correction, quantile normalization, and log2-transformation. The clinical endpoint was defined as a time to biochemical recurrence (BCR) or time to the latest follow-up (censored data).

TCGA PRAD cohort (6) is composed of primary tumor prostatectomy samples (n = 384), profiled on Illumina HiSeq 2000, with its RNA-sequencing and clinical data obtained from dbGap phs000178 and cBioPortal. TCGA dataset was mapped to hg19 using the STAR aligner (9). Raw counts were normalized, and variance was stabilized using DEseq2 normalization (10). Two clinical endpoints were considered: (i) time to BCR (defined as a rise of prostate-specific antigen in patient’s blood)/latest follow-up (censored data); and (ii) time to progression-free survival (PFS, defined as the progression of cancer, local recurrence, distant metastases, or PC-specific death)/latest follow-up.

The DKFZ PRAD cohort (8) is composed of primary tumor prostatectomy samples (n = 100) profiled on Illumina HiSeq 2000, with RNA-sequencing and clinical data obtained from cBioPortal. Gene expression data were obtained from cBioPortal as RPKM values, z-scored relative to all samples. The clinical endpoint was defined as a time to BCR or time to the latest follow-up (censored data).

To identify PC-specific biological aging genes associated with prostate cancer progression, we employed Cox proportional hazards analysis using the coxph function from the R survival package (11) on the Taylor et al. training cohort. The 52 biological aging genes were used as input variables, and time to BCR was used as a clinical endpoint in the univariable Cox proportional hazards analysis on the Taylor et al. patient cohort, adjusted for Gleason score (combined major plus minor Gleason was used for adjustment purposes) and chronological age variables (adjusted hazard p-value of < 0.05 was considered significant, Supplementary Table S2 , Supplementary Materials ).

To assess for multicollinearity among the identified (four) PC genes, we employed the variance inflation factor (VIF) analysis using the vif function from the R car package (11). The threshold of 5 was utilized to evaluate the presence of multicollinearity (all evaluated genes had values below 1.5).

For multivariable additive Cox proportional hazards modeling, expression levels of four identified genes were utilized as input variables and were adjusted for Gleason score and chronological age (as above). An adjusted Wald test p-value of < 0.05 was considered significant.

To investigate the predictive ability of the biological aging genes in subsets of patients defined by Gleason scores and chronological age, we employed cohort stratification. For this, the training Taylor et al. cohort was divided into two Gleason score groups: a group consisting of patients with Gleason scores 3 + 3 and 3 + 4, considered least aggressive clinically (n = 94), and the group consisting the rest of the patients (n IL = 44). For chronological age-stratified analysis, the Taylor et al. cohort was divided into two age groups: a group consisting of patients aged less than equal to a median age (in the Taylor et al. training cohort) of 60 (n = 87) and a group consisting of patients aged over 60 (n = 51).

Multivariable additive Cox proportional hazards modeling was run utilizing four identified genes as input variables and time to BCR as a clinical endpoint for each stratified group separately.

For testing/validation purposes, a multivariable additive Cox proportional hazards model in a stratified setting was performed using TCGA PRAD and DKFZ PRAD patient cohorts. For Gleason stratification (similarly to the training cohort), two subgroups were defined: the group that included patients with Gleason 3 + 3 and 3 + 4 (TCGA n = 159, DKFZ n = 68) and the group that included the rest of the patients (TCGA n = 225, DKFZ n = 32). For chronological age stratification, two subgroups were defined as in training: the group that included patients aged ≤ 60 (TCGA n = 175, DKFZ n = 100) and the group that included patients aged > 60 (TCGA n = 209, DKFZ = NA as this dataset only included patients younger than 6 years of age). The Cox proportional hazards model and p-values were estimated as above (i.e., during training).

AML gene expression patient cohorts utilized in this study included (i) the BEAT AML cohort (5), used for model training/discovery; and (ii) TCGA AML (also referred to as TCGA LAML) (6), used for model test/validation ( Table 2 ).

The BEAT AML (5) dataset is composed of bone marrow samples from de novo AML (n = 127), profiled on Illumina HiSeq 2500, with its RNA-sequencing and clinical data obtained from dbGap phs001657.v2.p1 and downloaded through the Genomic Data Commons (GDC). BAM files were converted to FASTQ files, which were then mapped to hg38 using the STAR aligner (9). Raw counts were normalized, and variance was stabilized using DEseq2 normalization (10). Overall survival was utilized as a clinical endpoint.

TCGA AML (6) dataset is composed of bone marrow samples from de novo AML (n = 151), profiled on Illumina HiSeq2000, with its RNA-sequencing and clinical data obtained from TCGA phs000178.v11.p8 and downloaded through the GDC. BAM files were converted to FASTQ files, which were then mapped to hg38 using the STAR aligner (9). Raw counts were normalized, and variance was stabilized using DEseq2 normalization (10). Overall survival was utilized as a clinical endpoint.

To identify AML-specific biological aging genes associated with AML progression, we employed Cox proportional hazards analysis using the coxph function from the R survival package (11) as above on the BEAT AML training cohort. The 52 biological aging genes were used as input variables, and overall survival was used as a clinical endpoint in the univariable Cox proportional hazards analysis, adjusted for cytogenetic risk score (based on ELN2017 guidelines (12), which utilize cytogenetics and mutation patient profiles and includes favorable, intermediate, and poor risk scores) and chronological age variables (an adjusted hazard p-value of< 0.05 was considered significant, Supplementary Table S3 , Supplementary Materials ). To ensure an appropriate sample size for each group, we combined intermediate and poor cytogenetic risk score categories into the “unfavorable” category.

Multicollinearity was assessed as above, with all evaluated genes having values below 1.5. For multivariable additive Cox proportional hazards modeling, expression levels of three identified genes were utilized as input variables and were adjusted for cytogenetic score and chronological age (as above). An adjusted Wald test p-value of < 0.05 was considered significant.

To investigate the predictive ability of the biological aging genes in subsets of patients defined by cytogenetic risk and chronological age, we employed cohort stratification. For this, the training BEAT AML cohort was divided into two cytogenetic risk groups (based on ELN2017 guidelines): consisting of patients with favorable cytogenetic risk (n = 53) and unfavorable cytogenetic risk (n = 74, which included patients with intermediate and poor cytogenetic risk scores, for sample size purposes). For chronological age-stratified analysis, the BEAT AML cohort was divided into two age groups: patients with a median age of ≤ 56 in the BEAT AML training cohort (n = 64) and patients aged > 56 (n = 63).

Multivariable additive Cox proportional hazards modeling was run utilizing three identified genes as input variables and overall survival as a clinical endpoint for each stratified group separately.

In AML, for testing/validation purposes, a multivariable additive Cox proportional hazards model in a stratified setting was performed using TCGA AML patient cohort. For cytogenetic stratification (similarly to the training cohort), two groups were defined: the group that included patients with favorable cytogenetic risk (n = 30) and a group that included patients with unfavorable cytogenetic risk (n = 121). For chronological age stratification, two groups were similarly defined: a group that included patients aged ≤ 56 (n = 76) and a group that included patients aged > 56 (n = 75). The Cox proportional hazards model and p-values were estimated as above (i.e., during training).

Since the 52 biological aging genes were nontissue-specific, our first objective was to assign their tissue specificity for prostate cancer. For this, we employed the Taylor et al. patient cohort (n = 138, see Methods), which comprised primary prostatectomy samples from patients with prostate cancer, with subsequent disease progression monitoring (defined as time to BCR). Taylor et al. dataset was specifically selected for training purposes as it presented a wide chronological age range (37.3–83.0) and contained 25% of patients with BCR failure events, making it an ideal dataset for training purposes. We utilized 52 genes as input variables in the univariable Cox proportional hazards model analysis, with time to BCR or the latest follow-up (if censored) as a clinical endpoint ( Figure 2A , left). This analysis identified 17 genes with unadjusted Cox proportional hazards p-value of< 0.05 ( Figure 2A , middle) that are significantly associated with prostate cancer progression in a tissue-specific manner. To ensure that our results are not due to the effect of the initial disease aggressiveness (i.e., Gleason score at diagnosis) or the effect of chronological age (at diagnosis), we adjusted our Cox proportional hazards model for these variables and identified four (4) genes that remained significant (adjusted Cox proportional hazards p-value of< 0.05) after adjusting for both Gleason score and chronological age ( Figure 2A , right; Supplementary Table S2 , Supplementary Materials ).

Four genes ( Figure 2B ) included those with HR ≤ 1 (negative association with BCR; the higher the gene expression, the poorer the prognosis the patient will have), such as CD44 (Cox proportional hazards model HR = 0.49, p-value of< 0.001), GADD45B (Cox proportional hazards model HR = 0.55, p-value = 0.008), STAT 3 (Cox proportional hazards model HR = 0.66, p-value = 0.017), and HR > 1 (positive association with BCR, the higher the gene expression, the more favorable the prognosis the patient will have) with BCR, such as GFAP (Cox proportional hazards model HR = 1.30, p-value = 0.036), demonstrating their significant association with BCR after subtracting the effect of Gleason score and chronological age ( Supplementary Table S2 , Supplementary Materials ).

Furthermore, to evaluate if the four genes could be utilized as a group for clinical purposes, we first assessed if there exists a multicollinearity among them using VIF analysis (see Methods, no multicollinearity identified), followed by the multivariable additive Cox proportional hazards model analysis (adjusted for Gleason and chronological age, as above), which demonstrated a significant benefit of evaluating the four genes as one additive model (Wald test p-value = 2.014−e9).

Finally, to evaluate pathway membership for the four identified genes, we employed the DAVID tool for pathway enrichment analysis (13) and identified two significant molecular pathways: the JAK-STAT and the FoxO signaling pathways (FDR-corrected chi-square test p-value of< 0.05). The JAK-STAT signaling pathway is involved in a wide range of cellular processes, including cell proliferation, differentiation, and survival (17, 18). The FoxO signaling pathway has been shown to regulate cellular processes such as apoptosis, cell cycle arrest, and DNA damage repair, which are important in maintaining genomic stability and preventing cancer development (19, 20). Both of these pathways have previously been shown to be implicated in biological age-related deterioration and disorders (21–23), and their further investigation might provide valuable insights into the underlying molecular mechanisms involved in prostate cancer biological age-related changes.

To further investigate if the significant association of the four biological aging genes with BCR is attributed to specific patient subgroups (that is driving significance of the association), we have performed stratified Cox proportional hazards analysis on the patient subgroups stratified by the (1) Gleason score, which is a commonly used prognostic factor for prostate cancer, and (2) chronological age, which is a known diagnostic factor for prostate cancer ( Figure 3A ).

For stratified analysis based on the Gleason score, we divided the Taylor et al. patient cohort into two patient subgroups: (i) lower Gleason score (3 + 3 and 3 + 4, commonly utilized in clinic as an indicator of favorable prognosis, n = 94) and (ii) a higher Gleason score (higher than 3 + 4, n = 44). A multivariable additive Cox proportional hazards model utilizing four identified genes as predictors and time to BCR as response was then performed for each subset separately ( Figure 3B ). This analysis demonstrated a significant association with BCR in a subset of patients with higher Gleason scores (Cox proportional hazards model Wald test p-value = 0.007).

For stratified analysis based on chronological age, we divided the Taylor et al. patient cohort into two patient subgroups: (i) younger chronological age (≤ 60, n = 87) and older age (> 60, n = 51), where 60 years of age is a commonly utilized threshold to estimate prostate cancer risk and also corresponds to the mean and median for the chronological age for this cohort. Multivariable additive Cox proportional hazards model on each chronological age subgroup separately ( Figure 3C ) demonstrated a significant association between the four identified biological aging genes and BCR (Cox proportional hazards model Wald p-value = 0.019 for younger chronological age group and p-value = 0.0045 for older chronological age group).

To validate our findings in independent patient cohorts, we have utilized (i) TCGA PRAD patient cohort (n = 384, see Methods ) and (ii) the DKFZ PRAD patient cohort (n = 100, see Methods), both with comparable to the Taylor et al. sample collection protocol (i.e., primary prostatectomy samples) and clinical endpoint (i.e., time to biochemical recurrence).

For testing the predictive ability of the identified four biological aging genes in TCGA PRAD dataset, it was stratified on Gleason score, similarly to training, into (i) Gleason scores 3 + 3 and 3 + 4 (n = 159); and (ii) Gleason scores higher than 3 + 4 (n = 225) ( Figure 4A ). Similarly to the Taylor et al. dataset, chronological age stratification on TCGA PRAD cohort defined (i) a group characterized by age ≤ 60 (n = 175) and (ii) a group characterized by age > 60 (n = 209) ( Figure 4A ).

Stratified analysis on the Gleason score demonstrated a significant association of the four biological aging genes to BCR in the higher Gleason score group in TCGA PRAD cohort (Cox proportional hazards model Wald test p-value = 0.049, Figure 4B ). Stratified analysis on chronological age, demonstrated a significant association of four biological aging gene panel with BCR in the younger age (≤ 60) in TCGA PRAD dataset (Cox proportional hazards model Wald test p-value = 0.029, Figure 4C ).

To confirm these findings in another PC patient cohort, we utilized the DKFZ PRAD dataset that was similarly stratified on Gleason score into (i) groups with Gleason scores 3 + 3 and 3 + 4 (n = 68); and (ii) groups with Gleason scores higher than 3 + 4 (n = 32) ( Figure 4A ). The DKFZ PRAD cohort only included patients aged below or at 60 (n = 100). The stratified analysis demonstrated significant association of the four biological aging genes to BCR in the higher Gleason score group (Cox proportional hazards model Wald test p-value = 0.011, Figure 4B ) and in the younger age group (≤ 60) (Cox proportional hazards model Wald test p-value = 0.0019, Figure 4C ), confirming results from TCGA PRAD dataset.

To make final conclusions, we have further integrated our findings from all three cohorts (Taylor et al., TCGA PRAD, and DKFZ PRAD) using FDR-corrected harmonic mean method ( Figures 4D, E ; Methods ), which nominated a significant association of the identified genes with BCR in the PC patient group with a higher Gleason score (FDR-corrected integrated p-value = 0.023) ( Figure 4D ) and younger age (≤ 60) subgroup (FDR-corrected integrated p-value = 0.019) ( Figure 4E ), indicating that the four biological aging genes could be utilized as markers to predict BCR in patients with a higher Gleason score and with a younger age, potentially allowing for personalized disease management and therapeutic planning in PC patients.

Furthermore, we expanded our investigations in TCGA dataset to evaluate if our four genes could predict PFS (defined as the progression of cancer, local recurrence, distant metastases, or PC-specific death). Our analysis demonstrated the ability of the four identified genes to predict PFS in the whole TCGA PRAD cohort in an additive Cox proportional hazards model (Wald test p-value = 0.034). Even though we observed a significant association with PFS, we anticipate that different biological aging genes might be involved in processes that govern BCR, PFS, or OS, requiring further investigations.

First, to ensure that the biological aging genes outperform the predictive ability of the genes chosen at random, we built a random model using the Taylor et al. cohort for prostate cancer and the BEAT AML cohort for AML. In such a model, genes (either four or three, respectively) were selected at random from the list of all available genes in each dataset and subjected to Cox proportional hazards model analysis as a group. The Wald p-value was reported, and the process was repeated 10,000 times ( Figure 6A ). For prostate cancer, this analysis confirmed a highly nonrandom ability of the identified four biological aging genes to predict prostate cancer progression ( Figure 6A , left, nominal random model p-value = 0.001), and for AML, it confirmed a highly nonrandom ability of the three identified genes to predict AML overall survival ( Figure 6A , right, nominal random model p-value = 0.048).

In PC analysis, to compare the predictive ability of the identified four biological aging genes to the predictive ability of the known transcriptomic markers of prostate cancer progression (27–39) ( Figure 6B ), we utilized Cox proportional hazards model analysis on the Taylor et al. dataset. Known prostate cancer genes associated with cancer progression were utilized as input variables, adjusted for Gleason score and chronological aging, with time to BCR as a clinical endpoint. The only known PC aggressiveness gene that remained significant after adjusted analysis was FOXM1 (adjusted Cox proportional hazards model p-value = 0.048), which was less significant compared to the identified four biological aging genes ( Figure 5B , right-most genes: CD44 adjusted p-value = 0.00054, GADD45B adjusted p-value = 0.0084, STAT3 adjusted p-value = 0.0167, GFAP adjusted p-value = 0.036). Following this result, we investigated the possibility of utilizing FOXM1 [previously identified by our group as a transcriptional regulator of prostate cancer progression (40, 41)] and the four biological aging genes together for enhanced predictive ability using a multivariable Cox proportional hazards model. Interestingly, when combined into a group of five, in the crude multivariable analysis, FOXM1 remained significant (hazard p-value = 0.01) but did not significantly reduce the predictive ability of the four biological aging genes, indicating their independent ability to predict prostate cancer progression. Yet, in the adjusted multivariable analysis for the group of five genes, FOXM1 lost its significance, while the biological aging genes remained significant. Furthermore, in the stratified multivariable analysis (i.e., analysis for the group of patients with higher Gleason scores and for the group of patients of younger age), FOXM1 did not show significance for the higher Gleason score (while the four biological aging genes remained significant) but demonstrated significance for the younger patient group, alongside biological aging genes. Taken together, this PC-specific analysis indicates that biological aging genes in PC offer independent predictive evidence compared to known markers of PC progression, especially for patients with higher Gleason scores, and could be effectively utilized alongside FOXM1 for patients of younger age for enhanced predictions of PC progression.

In the AML analysis ( Figure 6C ), we compared the predictive ability of the identified three biological aging genes to the predictive ability of the known markers of AML aggressiveness (42). Several genes have demonstrated significance in the adjusted analysis: IL-1R2 (adjusted Cox proportional hazards model p-value = 0.0024), RBP7 (adjusted Cox proportional hazards model p-value = 0.0269), and ZNF750 (adjusted Cox proportional hazards model p-value = 0.0467). Following this result, we investigated the possibility of utilizing these three genes (IL-1R2, ZNF750, and RBP7) and the three biological aging genes together to enhance their predictive ability using a multivariable Cox proportional hazards model. Interestingly, when utilized together, IL-1R2, ZNF750, and RBP7 showed significant ability to predict overall survival in a younger age group, as opposed to the CDC42EP2, CDC42, and ALOX15B, which showed a stronger predictive ability in older age group, indicating that both gene groups provide independent evidence and predictive ability for different patient populations. Taken together, this AML-specific analysis indicates that biological aging genes in AML offer independent predictive evidence, compared to known markers of AML progression, especially for older patients.