We retrieved mRNA expression profiles and clinical data of AML from TCGA. An overall count of 170 AML patient were included, out of which 139 patients had both clinical data and sequencing data available. We also collected normal bone marrow samples (n=70) from the Genotype-Tissue Expression databases. For analysis, mRNA expression was in HTSeq-FPKM format (level 3) and normalized to TPM reads.

The online NetworkAnalyst (https://www.networkanalyst.ca/) (Zhou et al., 2019) was employed to assess various mRNA levels of HERVs in AML, comparing data from the TCGA database with normal samples sourced from the GTEx database. The HERVs were identified following the criteria of P value< 0.05 and absolute log2-fold change > 1.

Survival analyses were performed by Kaplan-Meier estimates and followed by Cox regression model. Patients were divided into groups using the median expression level of Suppressyn and Syncytin-2 as the threshold.

Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were carried out to perform functional enrichment analyses by utilizing R clusterProfiler (version 4.4.4). The significance criteria were set at P value< 0.05 and absolute log2-fold change > 1. The AML cohort obtained from TCGA was divided into two groups based on the median expression scores of Suppressyn and Syncytin-2. The data was visualized using R ggplot2 (version 3.3.6). Moreover, gene set enrichment analysis (GSEA) was conducted using R gplot2 (version 3.3.6) and clusterProfiler (version 4.2.1). Function or pathway terms were considered significantly enriched in a statistical context, if they had an adjusted p-value< 0.05 and a false discovery rate (FDR)< 0.25.

The Spearman’s correlation analysis was employed to find out whether the diverse expression of Syncytin-2 and Suppressyn correlates with the ten most significant differentially expressed genes (DEGs). The DEGs were used to construct a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING) online database. The confidence score was considered high when greater than 0.7. Default values were used for all other parameters. Then, the resulting was edited via Cytoscape (version 3.9.1), a software for visualizing networks(Lotia et al., 2013). Ten most significant hub genes were identified using Cytoscape plugin CytHubba(Chin et al., 2014).

Immune Infiltration in AML was computed and analyzed through the ESTIMATE package in R, using ImmuneScore, StromalScore, and ESTIMATEScore(Yoshihara et al., 2013). We assessed the level of immune infiltration by examining 24 immune cells. Single-sample Gene Set Enrichment Analysis (ssGSEA) was able to compute the proportional enrichment of these cells in AML, implemented under the R package GSVA(Bindea et al., 2013). A Spearman’s correlation analysis was used to determine the relationship between Suppressyn, Syncytin-2, and these immune cells. The disparity in immune infiltrates between diverse expression patterns of Suppressyn and Syncytin-2 was evaluated using Wilcoxon rank-sum tests. The connection between immune checkpoints and expression of Suppressyn and Syncytin-2 were further examined by R package ggplot2 (version 3.3.6), aiming to explore their correlation with tumor immunity.

A multivariable Cox regression was used to determine the coefficients for DE-HERVs that showed statistical significance in univariable Cox regression. A risk-score formula was constructed as follows:

where N = 5, the Expi denotes the expression value of every five HERV-related genes, and the Coei represents the corresponding coefficient obtained from the multivariable Cox regression. The R package ggplot2 was employed by visualization. Based on the risk scores obtained, AML patients were divided into high- and low-risk groups using the median risk score as a threshold. The overall survival (OS) analysis was conducted between these two groups. We evaluated performance of the model in terms of accuracy of prediction by the receiver operating characteristic (ROC) curves.

The significance of diverse expression of Suppressyn and Syncytin-2 was assessed by the Wilcoxon rank-sum test. The association of clinical parameters and Syncytin-2 expression was explored by the Wilcoxon rank-sum test and logistic regression. These analyses were two-sided, conducted by R software version 3.6.3., and statistical significance was defined as P values< 0.05.

In the TCGA dataset we collected 173 AML patients, while in the GTEx dataset, we obtained 70 normal samples. Using web-based tool NetworkAnalyst, we identified a total of four DE-HERVs (ERVW-1, Syncytin-2, ERV3-1, ERVMER34-1) based on the cut-off criteria of absolute log2-fold change > 1 and an adjusted p-value< 0.05 ( Figure 1A ). In AML patients, heatmaps and violin plots demonstrated that 4 HERV-related genes were significantly upregulated compared to normal tissue ( Figures 1B, C ). Among them, ERVW-1, Suppressyn and Syncytin-2 displayed excellent predictive performance in distinguishing AML from normal samples, with AUC values of 0.977 (95% confidence interval [CI]=0.958-0.996), 0.990 (95% CI =0.977-1.000) and 0.965 (95% CI =0.940-0.990), respectively, by receiver operating characteristic (ROC) curve analysis ( Figure 1D ).

Next, we applied the Kaplan-Meier approach to observe whether diverse expression of DE-HERVs could affect patient survival in the progression of AML. Patients were classified into high and low expression clusters, stratified according to median values. Our results displayed that increased Suppressyn indicated a significantly improved prognosis for OS (HR = 0.65, 95% CI = 0.42-0.99, P = 0.045). Similarly, higher levels of Syncytin-2 expression were linked to an enhanced prognosis for OS (HR = 0.58, 95% CI = 0.38-0.90, P = 0.014) ( Figure 1E ). However, ERVW-1 and ERV3-1 present no significant correlation in terms of patient survival, as depicted in the corresponding Kaplan-Meier curves.

According to these results, the expression levels of Suppressyn and Syncytin-2 could function as essential indicators for AML diagnosis and prognosis.

TCGA cohors were used to investigate how different clinical features of AML samples ( Table 1 ) influenced the transcription of the two DE-HERVs. For Suppressyn, we observed statistically significant differences in the following categories: age groups (<= 60 vs > 60, P = 0.01), race groups (Asian and Black or African American vs White, P = 0.03), and NPM1 mutation (negative vs positive, P = 0.02) ( Figure 2A ). When dividing Syncytin-2 between high and low expression groups, significant statistical differences were found between subgroups in terms of cytogenetics (inv(16) & t(8;21) & t(15;17) vs normal, P = 0.03), cytogenetics risk (favorable vs intermediate/normal, P = 0.02), and NPM1 mutation (negative vs positive, P = 0.04) ( Figure 2B ). These findings supported that the transcription levels of Suppressyn and Syncytin-2 are associated with specific clinical features of AML, offering valuable insights into the potential involvement of these HERVs.

Furthermore, we evaluated their prognostic value in different subgroups of AML based on OS. Our findings showed that high Suppressyn expression consistently linked to favorable outcomes in subgroups, including race (white, HR = 0.55, 95% CI = 0.36-0.86, P = 0.009), BM blasts (<= 20%, HR = 0.43, 95% CI = 0.22-0.86, P = 0.018), cytogenetic risks (intermediate/normal, HR = 0.54, 95% CI = 0.31-0.94, P = 0.028), RAS mutation (negative, HR = 0.60, 95% CI = 0.39-0.93, P = 0.023), NPM1 mutation (negative, HR = 0.61, 95% CI = 0.37-1.00, P = 0.049) and IDH1 mutation (R172 negative, HR = 0.59, 95% CI = 0.38-0.90, P = 0.015) ( Figure 3 ).

Likewise, high Syncytin-2 expression linked to better survival outcomes in AML subgroups, including race (white, HR = 0.59, 95% CI = 0.38-0.92, P = 0.021), BM blasts (<= 20%, HR = 0.38, 95% CI = 0.19-0.77, P = 0.007), WBC counts (<= 20, HR = 0.43, 95% CI = 0.23-0.78, P = 0.006), cytogenetic risks (favorable & intermediate/normal, HR = 0.53, 95% CI = 0.32-0.89, P = 0.016), RAS mutation (negative, HR = 0.55, 95% CI = 0.35-0.86, P = 0.008), IDH1 mutation (R132 negative, HR = 0.64, 95% CI = 0.41-0.99, P = 0.046), IDH1 mutation (R172 negative, HR = 0.60, 95% CI = 0.39-0.92, P = 0.019), and IDH1 mutation (R140 negative, HR = 0.63, 95% CI = 0.40-0.99, P = 0.044) ( Figure 4 ).

Gene expression profiles of AML cohorts exhibiting diverse expressions of two DE-HERVs enabled us to pinpoint a distinct set of DEGs. The criteria were set at |Log2-FC| > 1 and p-value< 0.05. Volcano plots displayed DEGs between high and low Suppressyn groups ( Figure 5A ). The ten most significant DEGs (AC074389.2, CT45A1, AC109492.1, CT45A10, PPDPFL, LINC02059, TPRG1-AS1, RN7SKP169, MY018B and AC007091.1) are presented in Figure 5B and Supplementary Table 1 . Using the online STRING, we generated a network of protein crosstalk to explore predicted correlations between DEGs ( Supplementary Figure 1A ). Utilizing the same tool, our analysis also revealed from the PPI network an intricate group of hub genes, with the ten most significant identified as follows: IL10, CD4, ITGAM, CD86, CD163, MRC1, IL6, CD68, CCR5 and CCR1 ( Figure 5C ).

In a similar manner, we utilized volcano plots to compare DEGs in individuals with varying expression of Syncytin-2 ( Figure 5D ). We also delved deeper into the relationship between Syncytin-2 and the ten most significant DEGs, which consisted of PDPN, SLC10A2, RASL12, SLITRK6, AC116609.1, AL445209.1, NKX3-2, POU4F1, HMX3, and AL139351.2 ( Figure 5E and Supplementary Table 2 ). Moreover, we utilized STRING to develope a PPI network, with the objective of exploring possible interactions among DEGs ( Supplementary Figure 1B ). The ten most significant hub genes were identified as ADAMTS2, ADAMTS5, SEMA5B, ADAMTS15, ADAMTSL3, THSD7A, THSD4, ADAMTS20, LAMB3, and LAMB2 ( Figure 5F ).

Our analysis revealed that DEGs associated with Suppressyn showed significant enrichment in GO terms. These included the regulation of immune effector processes, collagen-containing extracellular matrix, and immunoglobulin binding. Moreover, the KEGG analysis highlighted significantly DEGs-enriched pathways, such as neutrophil extracellular trap formation ( Figure 6A and Supplementary Table 3 ). Furthermore, we conducted Gene Set Enrichment Analysis (GSEA) using GSEA/MSigDB, which demonstrated enriched immune-associated biological processes associated with Suppressyn. These processes included immunoregulatory interactions between a lymphoid and non-lymphoid cell, antigen activation of B cell receptors, and CD22 mediated BCR regulation ( Figure 6B and Supplementary Table 4 ).

When target switches to Syncytin-2 related DEGs, GO terms displayed a significant enrichment in response to cellular adhesion, such as molecules of plasma to membrane adhesion, collagen-containing extracellular matrix, and cell adhesion mediator activity ( Figure 6C and Supplementary Table 5 ). Similar aspects could be found in the KEGG pathway analysis. Interestingly, GSEA analysis also discovered the potential roles of Syncytin-2 in patterns of immunoregulation ( Figure 6D and Supplementary Table 6 ). Taken together, this functional exploration revealed both involvement of Suppressyn and Syncytin-2 in enhancing the immune phenotype in AML.

The ESTIMATE algorithm was applied to compute stromal, immune, ESTIMATE scores to determine the link between Suppressyn and Syncytin-2 and immune infiltrates in AML. Both Suppressyn and Syncytin-2 displayed a notable inverse correlation with immune scores and ESTIMATE scores, indicating a substantial influence on immune cell infiltration ( Figures 7A, B ).

In addition, we employed single-sample GSEA to demonstrate a positive link between Suppressyn with T-helper cells (r = 0.384) and Tcm (r = 0.257), while a negative link was oberved with macrophages (r = -0.333) and iDC (r = -0.331) ( Figure 7C ). With respect to Syncytin-2, we found it positively linked with multiple immune infiltrates, including T-helper cells (r = 0.385), Tcm (r = 0.335), and mast cells (r = 0.311), but a negative linke with macrophages (r = -0.304) ( Figure 7D ). All these findings were within strong significance, with P< 0.001.

Following this, we explored the relationship between Suppressyn and Syncytin-2 expression and immune checkpoints, grouped by immunoinhibitors and immunostimulators. Notably, both Suppressyn and Syncytin-2 showed similar expression patterns. They were positively correlated with immunoinhibitors, such as CD160, ADORA2A, CD274, CD96, TGFBR1 and CD244 ( Figures 8A, B ), while displaying a negative correlation with immunostimulators, including TNFSF13B, TNFSF13, CXCR4 and CD86 ( Figures 8C, D ), with more details showed in Supplementary Table 7 .

The ability to accurately classify patients as either high- or low-risk for AML progression is of utmost importance for patient management. After conducting an initial preselection step using univariate Cox regression analysis, Suppressyn and Syncytin-2 are significantly linked to patient survival (P< 0.05) ( Figure 9A ). To take into account the expression patterns and regression coefficients, we computed a risk score using the formula: Risk score = Suppressyn expression × (-0.19727) + Syncytin-2 expression × (-0.20209). According to their risk scores, AML patients sourced from the TCGA cohort were stratified as low- and high-risk groups ( Figure 9B ). Individuals in the low-risk category experienced lower mortality rates and longer survival periods compared to those in the high-risk category. According to the Kaplan-Meier curves, the high-risk category demonstrated a lower probability of survival in comparison to the low-risk category (HR = 1.83, 95% CI = 1.19 – 2.81, P = 0.006) ( Figure 9C ). The prognostic accuracy for OS stood at 0.611 after 1 year, 0.587 after 3 years, and 0.768 after 5 years based on time-dependent ROC analysis ( Figure 9D ). These results highlight the potential of the Suppressyn and Syncytin-2 -based risk-scoring model as a valuable tool for predicting AML patient outcomes.