Genome-wide identification and characterization of long noncoding RNAs (lncRNA) in individual immune cell lineages helps us better understand the driving mechanisms behind melanoma and advance personalized patient treatment. To elucidate the transcriptional landscape in diverse immune cell types of peripheral blood cells (PBC) in stage IV melanoma, we used whole transcriptome RNA sequencing to profile lncRNAs in CD4+, CD8+, and CD14+ PBC from 132 patient samples. Our integrative computational approach identified 27,625 expressed lncRNAs, 2,744 of which were novel. Both T cells (i.e., CD4+ and CD8+ PBC) and monocytes (i.e., CD14+ PBC) exhibited differential transcriptional expression profiles between patients with melanoma and healthy subjects. Cis- and trans-level coexpression analysis suggested that lncRNAs are potentially involved in many important immune-related pathways and the programmed cell death receptor 1 checkpoint pathways. We also identified nine gene coexpression modules significantly associated with melanoma status, all of which were significantly enriched for three mRNA translation processes. Age and melanoma traits closely correlated with each other, implying that melanoma contains age-associated immune changes. Our computational prediction analysis suggests that many cis- and trans-regulatory lncRNAs could interact with multiple transcriptional and posttranscriptional regulatory elements in CD4+, CD8+, and CD14+ PBC, respectively. These results provide novel insights into the regulatory mechanisms involving lncRNAs in individual immune cell types in melanoma and can help expedite cell type-specific immunotherapy treatments for such diseases.
Significance: These findings elucidate melanoma-associated changes to the noncoding transcriptional landscape of distinct immune cell classes, thus providing cell type-specific guidance to targeted immunotherapy regimens. Cancer Res; 78(15); 4411–23. ©2018 AACR.
Melanoma is responsible for the majority of skin cancer–related death with limited effective treatment available for stage IV patients (1). Despite recent treatment successes of immunotherapy and targeted agents, the majority of patients develop resistance to treatment for reasons that are still largely unknown. To better understand the underlying the disease mechanism, instead of traditional approaches focusing on antibody and cytokine production or the expression of selected cell-associated molecules, more comprehensive approaches to investigate the full spectrum of the transcriptome of individual cell types in peripheral blood cells (PBC) are desired (2).
Long noncoding RNAs (lncRNA) are highly heterogeneous molecules with functional versatility due to their diverse structures and interactions with other molecules (3). Recently, deregulation of lncRNAs has been reported to be greatly involved in various human diseases such as cancer (4–6) and heart disease (7, 8). LncRNAs are also widely expressed in immune cells serving as important regulators of gene expression throughout the immune system (9). However, the contribution of lncRNA-related biological activities in various immune compartments in response to disease treatment still remains little known. Integrative analysis of lncRNA expression profiling in diverse immune cell types (e.g., CD4+, CD8+, and CD14+ cells) will greatly expedite the discovery of lncRNA regulatory mechanisms and shed light on new treatment options for patients with melanoma.
Previous melanoma studies mostly focused on melanoma located in the skin tissue and the peripheral blood (10), whereas few have investigated genome-wide transcriptome profiling of individual immune cell types in PBCs of patients with melanoma (11). Furthermore, many studies only focused on protein-coding genes (12, 13), overlooking the vast landscape of noncoding genes such as lncRNAs. Hence, an in-depth investigation of lncRNA roles in individual immune cell types of PBCs will provide a much richer understanding of molecular regulatory mechanisms in patients with melanoma.
In this study, we developed an integrative computational functional genomic approach to interrogate the lncRNA expression profiles and their potential regulatory roles in CD4+, CD8+, and CD14+ PBCs of patients with stage IV melanoma. Our study presented a most comprehensive landscape of lncRNA profiling in CD4+, CD8+, and CD14+ PBCs of patients with melanoma, which holds great potential to provide clues for future cell type–specific immunotherapy treatment.
Materials and Methods
Sample collection and isolation of various immune cell types
Participant recruitment and blood sample collection were conducted using protocols approved by the Mayo Clinic Institutional Review Board (IRB#12-002580 and IRB#13-002293). Twenty-three normal healthy subjects were selected from a Southeastern Minnesota community cohort of adults self-identifying as having no autoimmunity, allergy, immunodeficiency, immunosuppressive treatments, or cancer. Eleven patients with a diagnosis of stage IV melanoma to receive immune therapy were selected. Of these donors, two normal healthy subjects and 10 patients with melanoma commonly contain CD4+, CD8+, and CD14+ PBC samples (Supplementary Table S1).
Isolation of CD4+, CD8+, and CD14+ PBCs
Blood was drawn into sodium heparin tubes (APP Pharmaceuticals, NDC 63323-540-11). CD4+, CD8+, and CD14+ cells were isolated from whole blood using anti-CD4/CD8/CD14 microbeads (130-090-877/130-090-878/130-090-879; Miltenyi Biotec) and automatics purification as per the manufacturer's instructions. A small portion of the cell isolate was set aside for flow cytometry using anti-CD3 APC-H7, anti-CD4 PE-Cy7, anti-CD11b-APC, anti-CD8 PE, anti CD14-FITC (Miltenyi Biotec) to assess the purity. All antibodies were purchased from BD Pharmingen unless otherwise indicated. Isolated cells were resuspended at 1 million/mL in RPMI containing 10% FBS and either immediately lysed in 0.7 mL QIAzol lysis reagent (Qiagen) or incubated for 4 hours at 37°C with 25 μL/mL anti-CD3/anti-CD28 Human T-Activator Dynabeads (111.32D; Invitrogen) for CD4 and CD8 cells or 1 μg/mL pI:C (528906; Calbiochem), LPS (L439; Sigma) CPG ODN2006; Invivogen), and PGN (Sigma; 77140) for CD14 cells. Lysates were stored at −80°C until processing for RNAseq.
RNA and cDNA library preparation
RNA was prepared using the Qiagen miRNeasy Kit according to the manufacturer's directions. RNA samples with integrity values of ≥6.0 by Agilent Bioanalyzer were processed into TruSeq libraries according to the manufacturer's instructions (RNA Prep Kit v2; Illumina) by the Mayo Clinic Medical Genome Facility Gene Expression Core. Paired-end DNA adaptors (Illumina) with a single “T” base overhang at the 3′end were immediately ligated to the “A-tailed” cDNA population. Unique indexes, included in the standard TruSeq Kits (12-Set A and 12-Set B), were incorporated at the adaptor ligation step for multiplex sample loading on the flow cells. Libraries (8–10 pmol/L) were loaded onto paired-end flow cells to generate cluster densities of 700,000/mm2 following Illumina's standard protocol for cBot and cBot Paired-End Cluster Kit version 3. The flow cells were sequenced as 51 × 2 paired-end reads on an Illumina HiSeq 2000 using TruSeq SBS Sequencing Kit version 3 and HCS v2.0.12 data collection software. Base-calling was performed using RTA version 18.104.22.168.
RNA-seq data analysis
Raw RNA-seq reads generated from each RNA-seq library were assessed for duplication rate and gene coverage using FastQC. The reads were aligned to the human reference genome (GRCh38) with the GENCODE annotation (gencode.v23.annotation.gtf) by TopHat (Version 2.0.12). The transcriptomes were then reconstructed by the Cufflinks tool (Version 2.2.1).
Identification of known and novel lncRNAs
We compared our reconstructed transcriptome with the GENCODE annotation V23 using the Cuffcompare script. Based on the comparison with the GENCODE annotation, we kept the transcripts in the “u” category (i.e., unknown intergenic transcripts), which were considered as candidates of novel long intergenic noncoding RNAs (lincRNA). Based on the comparison with the GENCODE lncRNA annotation, the transcripts in the “=” category were kept as known lncRNAs, and other transcripts in the non “=” category (i.e., the “i,” “j,” “o,” “p,” and “x” categories) were considered as candidate novel lncRNAs. These candidate novel lncRNAs were then filtered based on (i) the lncRNA features (i.e., sequence length larger than 200 nt and exon number greater than 2), and (ii) the coding potential using the Coding Potential Calculator algorithm (14), the Pfam database (15), and the Cording Potential Assessing Tool (16). The putative novel lncRNAs were then obtained (Fig. 1A).
Identification of expressed lncRNAs, uniquely expressed lncRNAs, and differentially expressed lncRNAs
For each PBC type, one lncRNA was considered expressed if its mean fragments per kilobase of transcript per million mapped reads (FPKM) value was greater than 0.1 in either healthy subject group or stage IV melanoma patient group; one lncRNA was considered uniquely expressed if it expressed only in patients with stage IV melanoma or normal healthy subjects; one lncRNA was considered differentially expressed if its expression level had greater than a 1.5-fold change and an false discovery rate (FDR) less than 0.1 between two groups using the Cuffdiff script.
The principle component analysis (PCA) was conducted using the prcomp function in the R software. Hierarchical clustering was performed using the pvclust R package (18) with the multiscale bootstrap resampling (i.e., 1,000 time), and P values were computed for each of the clusters (e.g., AU, approximately unbiased; BP, bootstrap probability).
The permutation test was used to evaluate the statistical significance of the epigenetic marks associated with lncRNAs. The P values of the permutation test were calculated by P value = (E + 1)/(R + 1), where R is the number of permutations (equal to 1,000 in our study), and E is the number of permutation test statistics (i.e., number of lncRNAs associated with selected epigenetic marks) that are greater than or equal to the observed test statistic (i.e., number of identified lncRNAs associated with epigenetic marks).
Pathway and network enrichment analysis
The Ingenuity Pathway Analysis (IPA) tool (QIAGEN) was used to annotate the genes of interest, map and generate putative biological processes/functions, networks, and pathways based on the manually curated knowledge database of molecular interactions extracted from the public literature. The enriched pathways and interaction networks were generated using both direct and indirect known relationships/connectivity. These pathways and networks were ranked by the enrichment score, which measures the probability that the genes of interest were included in a network by chance.
Weighted gene coexpression network construction and gene module detection
The R package “WGCNA” was used to construct the weighted gene coexpression network (19). The transcripts with the mean FPKM value great than 0.1 in CD4+, CD8+, and CD14+ PBCs of patients with stage IV melanoma or normal healthy subjects were used for the analysis. After transforming the FPKM value by log2(1+FPKM), weighted gene coexpression analysis was performed for CD4+, CD8+, and CD14+ PBCs subgroups, respectively. To minimize the noise and spurious associations, transcripts with similar expression pattern(s) were clustered using the module detection function “blockwiseModules” with the parameters networkType = “signed” and TOMType = “signed.” The expression profile of a given module was represented by its first principal component (i.e., module eigengene), which can explain the most variation of the module expression levels. Modules with highly correlated module eigengenes (correlation > 0.75) were merged together. The module membership (kME) of each gene was calculated by correlating the gene expression profile with the module eigengene, representing the extent of a gene close to a given module. The function “overlapTableUsingKME” was used to assess whether two modules were preserved based on a hypergeometric test of module kME.
Chromatin immunoprecipitation sequencing data analysis
We used the publicly available chromatin immunoprecipitation sequencing (ChIP-seq) datasets (CD4, CD8, and CD14 primary cells) from the Human Epigenome Atlas Release 9 (http://www.epigenomeatlas.org) generated by the Epigenomics Roadmap project (20). The histone mark peaks (i.e., H3K4me3 and H3K36me3) were detected using the MACS2 (version 2.1.1; ref. 21) with default parameters (FDR < 0.05). The liftover tool was used to convert the genome coordinates of the significant peaks from the GRCh37 to the GRCh38 annotation. We then identified the lncRNAs with H3K36me3 or H3K4me3 mark peaks in the region of ±3 kb transcription start sites of lncRNAs. We downloaded the super-enhancer datasets of the CD8 and CD14 PBCs from the dbSUPER database and identified the super-enhancer-associated lncRNAs as the ones whose transcription start sites were assigned to super-enhancers within a 50 kb window.
Transcription factor binding analysis
The human transcription factor (TF) binding motifs were retrieved from the ENCODE database (http://compbio.mit.edu/encode-motifs/). Using the FIMO tool (22), we searched TF binding motifs (P < 0.0001) within the regulatory regions of lncRNAs (i.e., 1 kb upstream and 1 kb downstream from the transcription start sites of lncRNAs).
Prediction of lncRNA interactions
We used the RNAplex software (23) to identify potential lncRNA–lncRNA interactions based on the RNA duplex energy prediction. The RNAplex parameter was set as -e-20.
To investigate potential interactions between miRNAs and lncRNAs, we downloaded all human miRNAs in the miRBase v21 and predicted potential target lncRNAs for each miRNA using the miRanda software (24) with default parameters.
To identify potential protein–RNAs interactions, we used the catRAPID algorithm (25) to predict potential protein–lncRNA interactions based on the information of protein and RNA domains involved in the macromolecular recognition. The protein–lncRNA interactions with the highly ranking distribution (i.e., the stringent parameters: the star rating > 2, RNA-binding motifs, and domains) were considered as significant hits.
The sequencing data set has been deposited in the NCBI Gene Expression Omnibus (GEO) database (GEO accession number GSE104744).
Systematic identification and characterization of lncRNAs in diverse immune cell types
To obtain a comprehensive lncRNA catalog of individual immune cell types, we developed an integrative transcriptome assembly pipeline designed to identify expressed lncRNAs in our dataset (Fig. 1A). In total, we reconstructed 345,417 transcripts from 65,247 loci across all 132 samples (Supplementary Table S1) using our transcriptome assembly pipeline. Of these transcripts, 24,881 were known lncRNAs (Supplementary Table S2), and 2,744 were putative novel lncRNAs (Supplementary Table S3; https://sites.google.com/view/yujizhanglab/data-scripts?authuser=0). Meanwhile, 19,476 protein-coding genes (98.7%) in the GENCODE annotation (V23) were also recovered in our assembled transcriptome. The high degree of overlap provides an independent measure of the comprehensive coverage in our dataset.
We observed that patients with stage IV melanoma and healthy subjects shared thousands of expressed lncRNAs in CD4+, CD8+, and CD14+ PBCs respectively (Fig. 1B), that is, the percentage of lncRNAs expressed in both groups is more than 86%. These lncRNAs were widely distributed across all the chromosomes in the human genome (Fig. 1C). We further investigated the differential expression profiling between patients with stage IV melanoma and healthy subjects. In CD4+ PBC, we identified 781 differentially expressed genes (including 98 lncRNA genes) and 445 differentially expressed transcripts (including 10 lncRNA). In CD8+ PBC, we identified 789 differentially expressed genes (including 88 lncRNA genes) and 292 differentially expressed transcripts (including 10 lncRNA). In CD14+ PBC, we identified 186 differentially expressed genes (including 22 lncRNA genes) and 32 differentially expressed transcripts (including two lncRNAs).
The cell type–specific expression patterns of lncRNAs
To investigate the cell-specific expression profilings of lncRNAs in CD4+, CD8+, and CD14+ PBCs, we compared their expression profiles between patients with melanoma and healthy subjects. Overall, lncRNAs were expressed at relatively lower levels in comparison with mRNAs in all three immune cell types. Approximately 75% lncRNAs showed low expression level (mean FPKM value < 1) in both patients with stage IV melanoma (Fig. 2A) and healthy subjects (Fig. 2B). We also observed that hundreds of lncRNAs were uniquely expressed in CD4+, CD8+, or CD14+ PBCs of patients with melanoma (Fig. 2C) and healthy subjects (Fig. 2D), respectively. In addition, there were many differentially expressed lncRNAs between CD4+, CD8+, and CD14+ PBCs (Fig. 2E and F), particularly between the T-cell subsets (CD4+ or CD8+ PBCs) and monocytes (CD14+ PBC).
To evaluate the cell type- and tissue-specific expression patterns of lncRNAs, we explored the expression of some known lncRNAs associated with melanoma in our dataset. First, we examined the expression profiles of lncRNAs differentially expressed in the skin tissue of melanoma. Among 426 differentially expressed lncRNAs in the skin tissue of melanoma cataloged in the MiTranscriptome database, 63 were also expressed in our dataset. However, they were not differentially expressed between patients with stage IV melanoma and healthy subjects in CD4+, CD8+, or CD14+ PBCs (Supplementary Fig. S1A). Second, we examined the expression of seven key immune-associated lncRNAs from one previous report (26) in our datasets (Supplementary Fig. S1B), all of which showed no significant expression differences between patients with stage IV melanoma and healthy subjects in CD4+, CD8+, or CD14+ PBCs. Third, we examined the expression of eight key melanoma-associated lncRNAs from the experimentally supported Lnc2Cancer database (27) and one previous study (28) in our dataset (Supplementary Fig. S1C). Out of eight lncRNAs, MALAT1 and SPRY4-IT1 were expressed in patients with stage IV melanoma (mean FPKM value > 0.1) and showed differential expression patterns between patients with stage IV melanoma and healthy subjects in each of three cell types. Taken together, lncRNAs exhibited tight cell-/tissue-specific expression patterns between patients with melanoma and healthy subjects in all three immune cell types.
lncRNAs cis-regulate their neighboring protein-coding transcripts enriched in important immune-associated pathways
To investigate potential cis-regulatory roles of lncRNAs on protein-coding transcripts, we computed the pairwise expression correlations between lncRNAs and their neighboring protein-coding transcripts (i.e., a neighboring protein-coding transcript needs to be within a 10 kb distance to a lncRNA). Overall, we identified 205, 236, and 228 highly correlated lncRNA–mRNA pairs (i.e., Spearman correlation coefficient |r| > 0.7, and P < 1E−07) in CD4+, CD8+, and CD14+ PBCs, respectively (Fig. 3A; Supplementary Table S4). Most of these lncRNAs are known lncRNAs, of which 7, 12, and 12 are novel lncRNAs in CD4+, CD8+, and CD14+ PBCs, respectively. Specifically, 192 lncRNA–mRNA pairs (94%) in CD4+ PBC, 223 lncRNA–mRNA pairs (95%) in CD8+ PBC, and 222 lncRNA–mRNA pairs (97%) in CD14+ PBC contained immune-associated genes based on the immunologic signatures in the MSigDB database (29).
The functional enrichment analysis by IPA revealed that the lncRNA–mRNA pairs in CD4+, CD8+, and CD14+ PBCs were significantly enriched in immune-related pathways (Fig. 3B; Supplementary Table S4). In CD4+ PBC, three top upstream regulators including TCF7 (P =1.40E−04), IL27 (P =1.93E−04), and EBI3 (P = 2.60E−04) were associated with the immune response pathways. In CD8+ PBC, the top upstream regulator TGFB1 (P =1.57E−04) belongs to the wnt/β-catenin signaling pathway [−log(P) = 2.07], which plays a critical role in T-cell immunity. In CD14+ PBC, the top upstream regulator MAPK8 (P = 1.43E−05) was associated with the immune response pathway.
We also investigated the relationship of the correlation and the transcription directions between lncRNAs and their neighboring mRNAs. Over 95% of lncRNA–mRNA pairs were positively correlated (Fig. 3C), consistent with the previous report on cis-acting noncoding RNAs (30). Over 70% of lncRNA-mRNA pairs had opposite transcriptional directions in CD4+, CD8+, or CD14+ PBCs (Fig. 3D). These findings are consistent with previous studies suggesting that lncRNAs have potential cis-regulatory roles on their neighboring protein-coding genes (31, 32).
lncRNAs trans-regulate protein-coding transcripts enriched in immune and mRNA translation pathways
To investigate potential trans-regulatory roles of lncRNAs on protein-coding transcripts, we performed a weighted gene coexpression network analysis (WGCNA; ref. 19) to identify groups of coexpressed transcripts (i.e., modules). We adopted the coding-to-noncoding strategy to identify the lncRNA–mRNA coexpression associations in a module (33): for a module of interest, we identified statistically significant correlations between each lncRNA and mRNAs that are important in specific known functional processes. The functions of each lncRNA were then inferred by the “guilty-by-association” approach (34).
Overall, we identified 13, 17, and 17 coexpression modules in CD4+, CD8+, and CD14+ PBCs, respectively (Supplementary Fig. S2A). Notably, four modules (|correlation coefficient| > 0.7, P < 1E−6) in CD4+ PBC, four modules (|correlation coefficient| > 0.68, P <1E−6) in CD8+ PBC, and one module (correlation = 0.74, P =10−8) in CD14+ PBC were statistically correlated with the melanoma status (Fig. 4A–C). Interestingly, these modules were also significantly correlated with age (|correlation coefficient| > 0.57, P < 2E−5; Supplementary Fig. S2A), whereas none of them were associated with the patient status (M1b, metastasis to the lung; M1c, other metastasis), disease response [stable disease (SD), partial response (PR); and progressive disease (PD); assessed at 3-months post-initiation of ipilimumab therapy), or sex.
We also explored lncRNAs associated with important immune processes in these nine coexpression modules (trans-regulatory lncRNAs). Specifically, we identified 361, 151, and 3 trans-regulatory lncRNAs in CD4+, CD8+, and CD14+ PBCs significantly associated with immune responses pathways respectively [−log(P) > 2; Supplementary Tables S5–S7]. In addition, these nine coexpression modules were also all enriched in three mRNA translation pathways, including EIF2 signaling, mTOR signaling, and regulation of eIF4 and p70S6K signaling pathway [−log(P) > 2; Fig. 4A–C]. In total, we identified 228, 158, and 8 trans-regulatory lncRNAs associated with these three pathways in CD4+, CD8+, and CD14+ PBCs, respectively (Supplementary Tables S5–S7). Of these trans-regulatory lncRNAs, 184 (81%) and 134 (85%) trans-regulatory lncRNAs were in turquoise modules of CD4+ and CD8+ PBCs (Supplementary Tables S5 and S6), suggesting that these modules contain major molecular regulators in these pathways.
Among these nine coexpression modules correlated with the melanoma status, the transcripts of six modules were significantly overlapped between CD4+ and CD8+ PBCs (P < 1E−10), the transcripts of two modules were significantly overlapped between CD4+ and CD14+ PBCs (P < 1E−10), and no transcript was significantly overlapped between CD8+ and CD14+ PBCs (Fig. 4D; Supplementary Fig. S2B). These results further confirmed that the CD4+ and CD8+ T cells shared more similar expression profiles than the CD14+ monocyte.
The lncRNAs associated with programmed cell death receptor 1 checkpoint pathways
Using the immune checkpoint blockers to block the interactions between programmed cell death receptor 1 (PD-1) and its ligands is a widely acknowledged approach for cancer immunotherapy. To investigate the impact of PD-1 checkpoint in diverse immune cell types, we explored PD-1 expression and related pathways, including PI3K/PTEN/Akt/mTOR signaling pathway and the RAS/RAF/MEK/ERK signaling pathway (35). We observed that the PD-1/PDCD1 gene was upregulated in T-cell subsets (CD4+ or CD8+ PBCs) but not monocytes (CD14+ PBC) of patients with melanoma compared with healthy subjects, dominantly expressed by the PDCD1-201 transcript (Supplementary Fig. S3A and S3B). We also identified 11 lncRNAs coexpressed with the PDCD1-201 transcript in CD8 PBC. In eight coexpression modules of T cells correlated with the melanoma status, 435 and 167 trans-regulatory lncRNAs associated with PD-1 checkpoint pathways in CD4 and CD8 PBCs, respectively (|correlation coefficient| > 0.7, P < 1E−6; Supplementary Table S8). Intriguingly, all these modules contain the mTOR signaling pathway (Supplementary Fig. S3C and S3D). The fact that 335 (77%) and 133 (80%) trans-regulatory lncRNAs were in turquoise modules of CD4+ and CD8+ PBCs suggested that the lncRNAs in these modules were potentially associated with PD-1 checkpoint pathways.
The potential regulation of lncRNA expression at the transcriptional level
To investigate potential regulatory mechanisms mediating lncRNA expression at the transcriptional level, we explored the association between the lncRNA regions and various biological annotations including histone modification sites, super enhancers, TF binding sites, and single-nucleotide polymorphisms (SNP).
The histone H3 lysine 4 trimethylation (H3K4me3) is a hallmark of actively transcribed promoters, and the histone H3 lysine 36 trimethylation (H3K36me3) is enriched not only within actively transcribed gene bodies but also active enhancers (36). To explore the potential epigenetic modification on lncRNA expression, we classified lncRNAs into the H3K4me3 and H3K36me3 marked groups based on the chromatin status of putative regulatory regions of lncRNAs in the ChIP-seq datasets generated by the Epigenomics Roadmap project (20). By focusing on the biologically significant lncRNAs in our dataset (i.e., immune-associated, translation-associated, and differentially expressed lncRNAs) in CD4+, CD8+, and CD14+ PBCs, we found more lncRNAs marked with H3K4me3 than those marked with H3K36me3 (Fig. 5A and B). CD4+ PBC had the largest number of lncRNAs commonly marked with H3K4me3 and H3K36me3 (Fig. 5C). Our permutation test suggested that our lncRNAs were indeed enriched in epigenetic marks (P < 0.001 based on 10,000 permutations; Supplementary Fig. S4 and Table S9). The fact that many lncRNAs were marked with the epigenetic modifications (e.g., H3K4me3 and H3K36me3), indicated that these lncRNAs were transcribed in CD4+, CD8+, and CD14+ PBCs.
Super-enhancers are large clusters of transcriptional enhancers that regulate expression of genes associated with specific diseases and cell identity (37). Within the biologically significant lncRNAs, we identified 105 lncRNAs annotated by the super-enhancers dbSUPER database (38) in CD8+ and CD14+ PBCs (Fig. 5D), suggesting that these lncRNAs may play important transcriptional regulatory roles in immune cells of stage IV patients with melanoma.
To investigate whether lncRNAs can be potentially regulated by TFs, we searched the regulatory regions of the biologically significant lncRNAs (i.e., upstream 1 kb and downstream 1 kb to the start sites of lncRNAs) in CD4+, CD8+, and CD14+ PBCs for potential TF binding motifs based on the annotation of the ENCODE database (Fig. 5E). We found that all except three of these lncRNAs can be potentially bound by more than 200 TFs (Fig. 5F). The fact that these lncRNAs have many TFs in CD4+, CD8+, and CD14+ PBCs suggests that they can be similarly regulated by TFs as mRNAs do.
The genome-wide association studies have identified thousands of SNPs associated with cancer susceptibility, over 90% of which are located in noncoding regions of the human reference genome (39). Using the noncoding somatic mutation annotation (CosmicNCV, release v80) from the COSMIC database (40), we explored whether the biologically significant lncRNAs contain any cancer risk-associated SNPs in CD4+, CD8+, or CD14+ PBCs (Fig. 5G). In total, 743 lncRNAs (93%) contained at least one cancer risk-associated SNP except the translation-associated lncRNAs in CD14+ PBC, 96% of which contain more than two cancer risk-associated SNPs and 14 contain more than 1,000 cancer risk-associated SNPs (Fig. 5H). These results suggested that most of these lncRNAs are potentially involved in cancer risk-associated diseases.
The potential regulation of lncRNA expression at the posttranscriptional level
To investigate potential posttranscriptional regulation through molecular interaction mechanisms, we explored the possibilities of RNA–RNA interactions, miRNA–lncRNA interactions, and protein–lncRNA interactions involving our lncRNAs.
First, we investigated the lncRNA–lncRNA interactions for the biologically significant lncRNAs in CD4+, CD8+, and CD14+ PBCs using a RNA duplex energy prediction approach (23). We found that these lncRNAs can potentially interact with each other, except two differentially expressed lncRNAs in CD4+ PBCs, three in CD8+ PBCs, and one in CD14+ PBCs (Fig. 6A), each of which has more than 20 potentially interacting lncRNAs and over 97% have more than 100 potentially interacting lncRNAs (Fig. 6B). Their potential interactions with each other suggest that some lncRNAs can potentially regulate the expressions of other lncRNAs and mRNAs.
The miRNA–lncRNA interactions can regulate lncRNAs either as inhibitory decoys or as regulatory targets of miRNAs (41). We investigated potential miRNA–lncRNA interactions between 2,588 annotated human microRNAs (miRBase v21) and our biologically significant lncRNAs in CD4+, CD8+, and CD14+ PBCs using the miRanda approach (24). We found that all the immune-associated and translation-associated lncRNAs have potential miRNA binding sites, except two differentially expressed lncRNAs in CD4+ PBCs, three in CD8+ PBCs, and one in CD14+ PBCs (Fig. 6C). Of these lncRNAs, each lncRNA had potential interaction with more than 80 miRNAs, over 98% of which have potential interactions with more than 100 miRNAs (Fig. 6D). These findings suggested that lncRNA expression can be potentially regulated by many miRNAs at the posttranscriptional level.
The protein–lncRNA interactions also play important roles in posttranscriptional regulation of the immune system (42). We used the catRAPID algorithm (25) to investigate whether there are any potential such interactions involving 16 most highly expressed lncRNAs (i.e., mean FPKM value >100 in patients with melanoma). The results showed that all 16 lncRNAs were potentially bound by multiple RNA-binding proteins (Supplementary Table S10). Of these lncRNAs, TCONS_00050166, TCONS_00212088, and TCONS_00050136 could potentially interact with 2204, 1802, and 414 RNA-binding proteins, respectively. These findings suggested that the protein–lncRNA interactions could serve as one important regulatory role for these lncRNAs.
Many patients with stage IV melanoma develop resistance upon their routine treatments such as immunotherapy and targeted therapy. The underlying mechanisms leading to resistance remain largely unknown. During the last decade, lncRNAs have been recognized as one of the largest regulatory RNA classes encoded in eukaryotic genomes. Our premise is that genome-wide analysis of lncRNA expression profiling in individual immune cell types in PBCs will expedite the understanding of the immune cell type–specific transcriptional regulatory mechanism in patients with melanoma, thus providing resources to investigate potential biomarkers or therapeutic targets for patients with melanoma in a noninvasive clinical setting in the future. As melanoma is often characterized by immune infiltration disease, we explored lncRNA expression patterns in CD4+, CD8+, and CD14+ PBCs of patients with stage IV melanoma. Overall, lncRNAs exhibited widespread expression patterns in the human genome and tight cell-type-specificity in CD4+, CD8+, and CD14+ PBCs. Our analysis suggested that lncRNAs can not only mediate the expression of protein-coding transcripts in a cis- or trans-regulatory manner in the immune system (26), but also be regulated by similar regulatory mechanisms as mRNAs do. For instance, almost biologically significant lncRNAs can be regulated by both transcriptional regulatory elements and posttranscriptional regulatory elements (e.g., miRNAs and lncRNAs; Fig. 7A). By comparing CD4+ and CD8+ PBCs, we identified 29 overlapping lncRNAs in T cells, including eight immune-associated lncRNAs, 20 translation-associated lncRNAs, and one differentially expressed lncRNA (Fig. 7B). TCONS_00060410, TCONS_00227034, and TCONS_00325476 were identified as both the immune-associated and translation-associated lncRNAs. Sixteen most highly expressed lncRNAs in CD4+, CD8+, and CD14+ PBCs were potentially regulated by multiple transcriptional elements and posttranscriptional elements (Fig. 7C). Through the functional enrichment analysis, we found that our lncRNAs were associated with several immune-related pathways in CD4+, CD8+, and CD14+ PBCs. Interestingly, most of the trans-regulatory lncRNAs involved in immune response pathways, PD-1 checkpoint pathways, and three mRNA translation processes were grouped in turquoise modules of CD4+ and CD8+ PBCs. Further functional investigation of these lncRNAs will be greatly appreciated to understand their biological roles in these pathways.
T cells (i.e., CD4+ and CD8+ PBCs) and monocyte (i.e., CD14+ PBC) exhibited differentially transcriptional expression profiles between patients with melanoma and healthy subjects (Supplementary Fig. S5 and S6). Specifically, the expression profiles of protein-coding transcripts in CD4+ and CD8+ T cells were able to distinguish patients with stage IV melanoma and healthy subjects. Furthermore, in nine coexpression modules significantly associated with the melanoma trait, CD4+ and CD8+ T cells share more overlapping modules than CD14+ monocytes.
Age is closely related to melanoma status in diverse immune cells. Using the coexpression network analysis, we found that age and melanoma traits were both significantly correlated in the nine coexpression modules in CD4+, CD8+, and CD14+ PBCs. Patient age has been reported as an independent prognostic factor for melanoma (43). With the increasing age, the immune system exhibits age-related changes (44), mainly due to the dysregulation of T-cell function resulting in T-cell immunosenescence. In T cells (CD4+ and CD8+ PBCs), the mTOR signaling pathway was enriched all eight coexpression modules. PD-1 is an immune-checkpoint receptor that mainly expresses in T cells negatively regulating human immune response (45). Inhibition of mTOR signaling pathway can extend the life span in many different species (46) and delay the onset of age-related diseases in mice (47). Because of potential age-related changes in immune cells (48), patients with melanoma in different age groups may need more personalized immunotherapeutic treatments.
There are some limitations in our study. First, CD4+, CD8+, and CD14+ PBCs primarily contain CD4 T cell, CD8 T cell, and CD14 monocyte, respectively. However, they are still heterogeneous cell types that may have distinct transcriptional expression profiles. Second, our RNA-seq dataset only contains lncRNAs with ploy-A tails, which account for a partial lncRNA transcriptome. Third, the prediction software in our analysis could introduce some computational errors. Although the multi-omics integrative analysis provides a wide understanding of the original cause of disease (genetic, environmental, or developmental; ref. 49), the limited sample size and expensive experimental cost prevent such multi-omics experiments to be conducted in individual diseases. As the cost of single-cell sequencing technology becomes more affordable (50), we expect that an integrative analysis of diverse immune cell types at single cell level will promote a much deeper understanding of the immune response mechanism and more personalized treatment to patients with melanoma.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Conception and design: L. Wang, S.J. Felts, V.P. Van Keulen, A.D. Scheid, L.R. Pease, Y. Zhang
Development of methodology: L. Wang, V.P. Van Keulen, Y. Zhang
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): V.P. Van Keulen, M.S. Block, S.N. Markovic, L.R. Pease, Y. Zhang
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis):L. Wang, V.P. Van Keulen, A.D. Scheid, L.R. Pease, Y. Zhang
Writing, review, and/or revision of the manuscript: L. Wang, S.J. Felts, V.P. Van Keulen, L.R. Pease, Y. Zhang
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): V.P. Van Keulen, S.N. Markovic, Y. Zhang
Study supervision: Y. Zhang
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.