Age-related transcriptome changes in melanoma patients with tumor-positive sentinel lymph nodes

Transcriptome changes in SLN genes in older patients (≥60 years old) versus younger patients (<60 years old) by microarray analysis

In our previous studies, we have found that expression of two SLN genes (PIGR and TFAP2A), when combined with clinicopathological features, correlated with prognosis in SLN-positive melanoma patients [12]. Since age is also an important factor in determining patient outcomes, we were interested in comparing gene expression profiles in older versus younger patients and in assessing whether there was a correlation with melanoma recurrence. Therefore, we analyzed the first microarray dataset from 97 melanoma patients with positive SLNs from the Sunbelt Melanoma Trial (SMT) and evaluated the transcriptome changes of the SLN by two defined age groups: the older and the younger groups. To ensure that we had a large enough sample size for a robust analysis, patients were defined as being older if they were ≥60 years old (yr⁶⁰⁺). Patients were defined as being younger if they were <60 years old (yr^60-).

Table 1 lists the clinical data of the 97 melanoma patients grouped by age. In this dataset, there were no significant differences between the two age groups in primary tumor site, Breslow thickness, Clark level, or ulceration presence. However, in younger patients, the recurrence rate was significantly higher when Breslow thickness was higher. In older patients, there were no significant differences in Breslow thickness, Clark level, and ulceration presence between groups of patients with recurrence (recur^yes) and those without recurrence (recur^no). Using microarray filter T3 and T4, we detected a total of 577 and 156 differentially expressed probe sets, in older versus the younger patients. Among them, there were 41 and 11 differentially expressed probe sets by filters T3 and T4 in the older versus younger groups (p <0.05). Probe sets without defined gene names by annotation from Partek Genomics Suite software were removed from the lists. There were 7 differentially expressed genes (DEGs) in the yr⁶⁰⁺ group versus the yr^60- group with a p value < 0.05 by T4 filter (Table 2). Among them, FBJ murine osteosarcoma viral oncogene homolog (FOS) and nuclear receptor subfamily 4, group A, member 2 (NR4A2), were the two genes that had significant higher expression in the yr⁶⁰⁺ group than in the yr^60- group. The DEGs between the yr^60- and the yr⁶⁰⁺ group had various biological functions, including toll-like receptor signaling pathway transduction, adaptive and innate immune response, autophagy, and transcription regulation (Table 2). The network connection of the 156 DEGs by T4 filter is shown in Supplementary Figure 1. The top canonical pathway that showed a difference in the yr^60- and the yr⁶⁰⁺ group was the peroxisome proliferator-activated receptor (PPAR) signaling pathway, which had a close interaction with toll-like receptor signaling pathway (Supplementary Table 1) [13, 14]. The list of all DEGs by T3 filter is listed in Supplementary Table 2.

Transcriptome changes of immune genes and immune pathway genes in the SLNs of older versus younger patients as assessed by NanoString analysis

Immune cells are a major component of lymph node structure. We then focused on immune genes and immune pathways associated with both age groups and assessed by NanoString analysis. This analysis in the second dataset found that 12 immune-related genes were differentially expressed in SLNs in older versus younger patients (Table 3). There were 17 immune pathway-related genes in SLNs that were differentially expressed in yr⁶⁰⁺ versus yr^60- patients (Table 4). Of note is that the NR4A2 gene was found to be differentially expressed in yr⁶⁰⁺ versus yr^60- patients from the first microarray dataset. The NR4A3 gene, which belongs to the same family members of NR4A2, was also found to have a higher fold change (FC) in yr⁶⁰⁺ patients, even though the p value is borderline (p=0.0517) (last row in Table 4). The immune gene, integrin subunit beta 1 (ITGB1), was found to be differentially expressed in yr⁶⁰⁺ versus the yr^60- patients (Table 3). Integrin subunit beta like 1 (ITGBL1) was also found to be differentially expressed in yr⁶⁰⁺ versus yr^60- patients by microarray analysis (Supplementary Table 2). The immune gene with the highest and lowest fold change in the yr⁶⁰⁺ versus the yr^60- patient group was melanoma antigen family A, 3 (MAGEA3) and leukemia inhibitory factor (LIF) (fold change =2.87 and -1.16) (Table 3). Among the three highest fold changes of the immune pathway genes, two of them were secreted frizzled-related protein 2 and 4 (SFRP2 and SFRP4) (fold change =1.93 and 1.78) (Table 4). Both genes belong to the Wnt pathway. A volcano plot of the immune genes and immune pathway genes that were differentially expressed in yr⁶⁰⁺ versus the yr^60- patients is presented in Supplementary Figures 2, and 3.

NanoString results suggested that NR4A and ITGB1 genes are highly expressed immune genes in older melanoma patients rather than their younger counterparts with lymph node metastasis. These genes, therefore, might be responsible for the age-related differences in response of SLN to the presence of nodal metastasis. The Wnt pathway might also be an important immune pathway associated with age-related immune response to melanoma metastasis to the SLN.

Transcriptome changes in SLN associated with recurrence in yr⁶⁰⁺ or yr^60- melanoma patients by microarray analysis

After we compared the transcriptome changes in SLN genes between the yr⁶⁰⁺ and the yr^60- melanoma patients, we studied whether there were any differences between patients who experienced recurrence versus those who remained disease free. We also evaluated these results by age categories. A multivariable linear regression model was fitted for each gene of each sample about age (yr⁶⁰⁺ or yr^60-), outcome (recur^yes or recur^no), and the interaction of age and outcome in the first microarray dataset. There were 100 differentially expressed probe sets with a statistically significant difference (p <0.05) after adjusting either by age or outcome or the interaction of age and outcome using filter T4. (Data not listed due to the long list of genes, but available upon request.) Among them, there were 11 differentially expressed probe sets with a significant difference adjusting by the interaction of age and outcome (p <0.05). Probe sets of the same gene were merged. There were 6 genes with statistically significant differences between groups (Table 5). We further analyzed the mean and 95% confidence interval (CI) of these 6 DEGs (Table 5). Means (95% CI) without overlapped values between each group were italicized. The non-overlapped values implied that there were statistically significant differences between the two groups. For example, NR4A2 was differentially expressed in yr⁶⁰⁺ versus yr^60- melanoma patients without recurrence. NR4A2 also showed significant differences in yr⁶⁰⁺ patients with (recur^yes) versus those without recurrence (recur^no) (Table 5).

Transcriptome changes of immune genes and immune pathway genes in SLNs associated with recurrence in yr^60- and yr⁶⁰⁺ melanoma patients by NanoString analysis

In the NanoString dataset, we first analyzed the differentially expressed immune genes between recur^yes and recur^no groups in younger melanoma patients (yr^60-). The results showed that there were 20 differentially expressed immune genes (p<0.05) in this comparison (Supplementary Table 3). Selected differentially expressed immune genes between the recur^yes and recur^no patients with p<0.05 and absolute fold change >0.5 in the yr^60- group were listed in Table 6. In yr^60- patients with positive SLNs, highly expressed C6, interleukin 23 receptor (IL23R), B melanoma antigen (BAGE), chemokine [C-C motif] ligand 16 (CCL16), and lower expression of S100 calcium binding protein B (S100B) were associated with recur^yes patients. A volcano plot of the immune genes that were differentially expressed in younger patients (yr^60-) in the recur^yes versus the recur^no groups is shown in Supplementary Figure 4.

In older patients, there were 20 differentially expressed genes between the recur^yes and recur^no group (p<0.05) (Supplementary Table 4). Table 7 lists the selected differentially expressed immune genes by recurrence status in the yr⁶⁰⁺ melanoma patients with p<0.05 and absolute fold change > 0.5. In yr⁶⁰⁺ patients with positive SLNs, highly expressed FOS and CCL18 were associated with recur^yes. A volcano plot of the immune genes that was differentially expressed in older patients in the recur^yes versus the recur^no group is shown in Supplementary Figure 5.

When comparing the difference in the DEGs by recurrence status in both age groups, we found that MAPK11 was highly expressed in the younger melanoma patients in the recur^yes versus the recur^no group (FC=2.84) (Table 6). A similar family member, MAP2K4, had marginal expression in older patients in the recur^yes versus the recur^no group (FC=0.25) (Supplementary Table 4). CCL16 had a higher expression in the younger patient cohort in the recur^yes versus the recur^no group (FC=3.46) (Table 6). Another family member, CCL18, also had a higher expression in older patients with recurrence (FC=1.8) (Table 7). C6 was highly expressed in younger melanoma patients with recurrence (FC=4.28) (Table 6), while C3 had marginal expression in older patients with recurrence (FC=0.83) (Table 7).

In terms of immune pathway genes, there were 18 differentially expressed genes with p<0.05 and absolute fold change > 0.5 in the younger patients when comparing recur^yes versus recur^no (Table 8). A complete list of the DEGs with p<0.05 is presented in Supplementary Table 5. Supplementary Figure 6 shows a volcano plot of the immune pathway genes in younger patients that were differentially expressed by recurrence status. In the group of older patients, there were 13 differentially expressed immune pathway genes with p<0.05 and absolute fold change > 0.5 by recurrence status (Table 9). Supplementary Figure 7 shows a volcano plot of the immune pathway genes that were differentially expressed in older patients based on recurrence status. All the DEGs with p<0.05 in the older group are listed in Supplementary Table 6. IRAK3 (interleukin-1 receptor-associated kinase 3) was the major immune pathway gene found in younger patients with recurrence (Table 8), while Wnt10b was the major pathway found in older patients with recurrence (Table 9). There were no overlapped immune pathway genes in either age group by recurrence status. These results suggested that, even though some immune genes have similar changes in older and younger patients, different pathways may be involved in recurrence in different age groups.

Verification of the DEGs in the third independent dataset by quantitative reverse transcriptase polymerase chain reaction (qRT-PCR)

After using microarray and NanoString technologies to identify the DEGs in the yr^60- and yr⁶⁰⁺ patients as well as the different outcome-associated age groups (recur^yes versus recur^no), we used qRT-PCR in another independent dataset to confirm the findings above. We selected genes that were differentially expressed in both microarray and NanoString analysis or had higher fold changes in either of the analysis. The results showed that FOS, NR4A2, PTGS2, LINC00518, IL1B, and Wnt10b were all highly expressed in older patients with recurrence (Table 10). These genes converged at the Wnt10b pathway (Figure 1).

Patient selection

This study used two different technologies in three independent datasets of RNA samples obtained from melanoma patients with positive SLNs to identify age-related transcriptome changes in SLN and their association with outcome.

Microarray analysis was performed in the first independent dataset to assess 97 samples obtained from the Sunbelt Melanoma Trial (SMT). The samples were randomly chosen from among 317 melanoma patients with positive SLNs. This patient cohort has been described previously [12]. Thirty-nine patients experienced recurrence melanoma in this cohort, and fifty-eight patients did not experience recurrence. Median follow-up was 93 months [12]. This study was approved by the institutional review boards (IRB) of each participating institution. Clinicopathologic factors, recurrence, and survival data were collected prospectively. Additional details of the SMT are described elsewhere [12, 32].

NanoString analysis was applied to the second patient cohort, which included 12 patients with tumor-positive SLNs from the James Graham Brown Cancer Center Biorepository at University of Louisville. This study followed an approved IRB protocol. There were 6 patients who experienced recurrence (3 of each at age <60 and ≥ 60 years old) and 6 patients who did not experience recurrence (3 of each at age <60 and ≥ 60 years old). Median follow-up was 34 months.

The third independent dataset of 36 samples from the James Graham Brown Cancer Center Biorepository was used to validate the differentially expressed genes (DEGs). The SLN tissue was acquired from patients at the time of surgical treatment of cutaneous melanoma, including staging with SLN biopsy between 2003 and 2017. Median follow-up of this cohort was 33.2 months. Patient characteristics such as age and outcome from all three datasets are summarized in Supplementary Table 7.

Definition of age groups

To ensure that we had a large enough sample size for a robust analysis, we grouped patients into two age groups. Patients were defined as being older if they were ≥60 years old (yr⁶⁰⁺). Patients were defined as being younger if they were <60 years old (yr^60-).

Microarray experiments

GeneChip Human HG-U133 plus 2.0 array (Affymetrix, Santa Clara, CA) was used in the first microarray dataset according to the manufacturer’s guidelines. Details of RNA isolation, microarray experiment, and quality control were described in detail previously [12]. This set of microarray data is accessible through NCBI’s Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) by accession number GSE 43081.

NanoString analysis of mRNA expression of immune panel genes and immune pathway panel genes

The second dataset of 12 RNA samples were isolated from fresh-frozen human SLN tissues from melanoma patients using RNeasy Plus Mini Kit (Qiagen). RNA quality control/quantity assessment (QC/QA) was checked by Agilent bioanalyzer. The RNA concentration was measured by Qubit. Total RNA (100 ng per sample) were analyzed on the nCounter MAX system. Two gene expression assays were used: PanCancer immune profiling and PanCancer immune pathway profiling (NanoString Technologies, Seattle, WA, USA). PanCancer immune profiling assay comprised 730 immune-related genes and 40 internal reference genes. Immune pathway profiling assay comprised 730 genes from 13 canonical pathways and 40 selected reference genes. Raw counts for each assay were collected using the NanoString data analysis software (nSolver).

Quantitative reverse transcriptase polymerase chain reaction (qRT-PCR)

The third dataset of 36 RNA samples were isolated from fresh-frozen human SLN tissues from melanoma patients using RNeasy Plus Mini Kit (Qiagen). Total SLN RNA (1000 ng) from each sample was reverse-transcribed with the SuperScript III First-Strand Synthesis System. mRNA primers were purchased from Life Technologies (Carlsbad, CA). Quantitative RT-PCR reactions were completed on a 7500 Fast Real Time PCR system (Life Technologies). The relative quantity of the target mRNA was normalized to endogenous gene (B₂M). The fold changes (FC) of each mRNA in the qRT-PCR experiments were calculated with the 2^-ΔΔCt method.

Statistical analysis

For microarray analysis, a fold change outlier (FCO) filter was applied independently to reduce the dimension of the data before determining the DEGs between the two age groups (yr⁶⁰⁺ and yr^60-) as well as between patients with recurrence (recur^yes) and those without recurrence (recur^no) [37, 38]. For each of 54,675 probes on the array, the fold change (FC) was calculated and four filters (T1, T2, T3 and T4) were used. T1={μ(FC) ± 1.5σ(FC)}, T2={μ(FC) ± 2σ(FC)}, T3={μ(FC) ± 3σ(FC)}, and T4={μ(FC) ± 4σ(FC)}, where μ(FC) is the mean of fold changes (FC) and σ is the standard deviation of FC from all 54,675 probes in the array. The genes that fell inside T1, T2, and T3 were filtered from the differential data. After filtering the data, a t-test for normal gene expression data and a Wilcoxon test for non-normal expression data were applied [39]. The Benjamini-Hochberg method [40] was employed to adjust the p values. When comparing the changes of the SLN gene expressions in the yr⁶⁰⁺ versus yr^60- patients, a multivariable linear regression model was fitted for each gene about age (<60 or ≥ 60). The equation used is below:

Gene Expression = α + β1 age, where age=1 if ≥60 years old and 0 otherwise.

The estimates and p values are presented by filter T2, T3, and T4. When assessing the changes of the SLN gene expressions that are associated with recur^yes versus recur^no in the yr⁶⁰⁺ and yr^60- melanoma patients, a multivariable linear regression model was fitted for each gene of each sample about age (<60 years or ≥60), outcome (recur^yes or recur^no), and the interaction of age and outcome. The equation used is below:

Gene Expression = α + β1 age+ β2 outcome+ β3 age*outcome, where age=1 if ≥60 years old and 0 otherwise, outcome=1 if recur^yes and 0 otherwise.

The estimate and p values are also presented by filter T2, T3, and T4. Statistical Analysis System (SAS) was used to perform the regression analysis [41, 42]. p values of FC were calculated using ANOVA.

For the NanoString results analysis, positive control normalization was performed by using gene expression data normalized to the mean of the positive control probes for each assay. RNA content normalization was performed by using gene expression data normalized to the geometric mean of housekeeping genes in the CodeSet. Raw data are also analyzed using the nSolver Advanced Analysis module. More information on the Advanced Analysis package can be found at http://www.nanostring.com/products/nSolver.

Ingenuity Pathway Analysis (IPA) software (Ingenuity Systems, Redwood City, CA) was used for gene network and pathway analysis. The statistical score of a pathway is defined as –log (p value) from Fisher’s exact test analysis.