Abstract
We employed whole-genome RNA-sequencing to profile mRNAs and both annotated and novel long noncoding RNAs (lncRNAs) in human naive, central memory, and effector memory CD4+ T cells. Loci transcribing both lineage-specific annotated and novel lncRNA are adjacent to lineage-specific protein-coding genes in the genome. Lineage-specific novel lncRNA loci are transcribed from lineage-specific typical- and supertranscriptional enhancers and are not multiexonic, thus are more similar to enhancer RNAs. Novel enhancer-associated lncRNAs transcribed from the IFNG locus bind the transcription factor NF-κB and enhance binding of NF-κB to the IFNG genomic locus. Depletion of the annotated lncRNA, IFNG-AS1, or one IFNG enhancer-associated lncRNA abrogates IFNG expression by memory T cells, indicating these lncRNAs have biologic function.
Introduction
Immunologic memory is a critical feature of the ad+aptive immune response conferring long-lasting protective immunity to the host. Memory T (TM) cells (CD4+ or CD8+) can be broadly divided into central memory T (TCM) and effector memory T (TEM) cells (1, 2). TCM cells are CD45RO+, constitutively express CCR7 and CD62L, recirculate through lymphoid organs, express high levels of IL-2 after restimulation, proliferate vigorously compared with TEM cells, express low levels of IFN-γ, and have limited effector function (3–5). TEM cells are also CD45RO+, but express low levels of CCR7 and expression of CD62L is heterogeneous. TEM cells preferentially reside in nonlymphoid tissues and rapidly produce effector cytokines in response to Ag stimulation (6). Despite our clear understanding of the critical importance of TM cells to the adaptive immune system, our understanding of mechanisms driving differentiation into TEM and TCM lineages and maintaining these phenotypes is incomplete (3, 7).
The vast majority of the genomic sequence is transcribed into either protein-coding RNAs or noncoding RNAs (ncRNAs) (8). Long ncRNAs (>200 bp to distinguish them from short RNAs such as micro-RNAs or tRNAs) can be subdivided into long noncoding RNAs (lncRNAs), enhancer RNAs (eRNAs) and other species. Certain lncRNAs exhibit many features of mRNAs and are 5′ capped, spliced, and polyadenylated but exhibit little if any coding potential (9–12). About 20,000 lncRNA genes have been identified in vertebrate genomes, similar to the number of protein-coding genes. A growing array of regulatory functions has been ascribed to this new class of RNAs (9, 10, 13). The eRNAs, so-called because they are transcribed near transcriptional enhancer elements and lack promoters in the classic sense, are divided into two classes: 1d-eRNAs and 2d-eRNAs (14–17). The 2d-eRNAs are transcribed in both sense and antisense orientations, tend to be relatively short (500–2000 bp), are relatively unstable, and are not polyadenylated. The 1d-eRNAs are transcribed in one direction only, are polyadenylated, and are relatively stable. Thus, if strand-specific information is available, it is relatively straightforward to distinguish 2d-eRNAs from 1d-eRNAs. At this point, a clear set of rules has not emerged to distinguish lncRNAs from 1d-eRNAs, although types of epigenetic modifications and cellular locations may aid classification. In contrast to lncRNAs, biologic functions of 1d- and 2d-eRNAs are less well understood and it has been argued that enhancer transcription may simply represent transcriptional noise without true biologic function ascribed to the transcribed RNA (18, 19).
DNA enhancer elements are defined by their ability to direct cell type–specific gene expression programs determining cell identity. Stretch-enhancers, also referred to as superenhancers (SE), are a subset of enhancers exhibiting proximity to certain protein-coding genes that drive cell identity (20, 21). SEs differ from typical enhancers (TE) in size, transcription factor density, and ability to activate transcription. Both are marked by histone acetyltransferases, p300 and CBP, and H3K27-acetylation (H3K27Ac) (22, 23). Bromodomain-containing proteins (BRD), such as BRD2, BRD3, BRD4 and BRDT, are recruited to TE and SE H3K27Ac marks to facilitate recruitment of RNA pol2 and activate transcription (24). Genes with SE structures are also enriched near single nucleotide polymorphisms linked to autoimmune disorders, including rheumatoid arthritis, and expression of risk genes is modulated by the JAK inhibitor tofacitinib, which is used to manage rheumatoid arthritis (22).
In our present study, we sought to identify mRNAs, and annotated and novel lncRNAs preferentially expressed in naive T (TN), TCM, and TEM cells. Our results support the idea that lineage-specific mRNAs, known lncRNAs, and novel lncRNAs are coexpressed, and these loci are colocalized in the genome with SEs and TEs. Further, inhibition of SE or TE function abrogates expression of annotated and novel lncRNAs and mRNAs at key gene loci required for lineage-specific cell phenotypes. One mechanism by which these novel lncRNAs act is to bind the transcription factor, NF-κB, and reinforce binding of NF-κB to chromatin. We also show that depletion of either the annotated lncRNA, IFNG-AS1, or one of the IFNG-associated novel lncRNAs abrogates IFNG transcription by TM cells, a critical function of the adaptive immune response.
Materials and Methods
Cell isolation and culture
Human PBMC were isolated from healthy control subjects with no chronic or acute conditions using Ficoll-Hypaque centrifugation. All subjects included in the study, both in initial RNA-sequencing (RNA-seq) and follow-up analyses, were male, of Caucasian descent, and between ages of 25–32. TN, TCM, and TEM cells used for RNA-seq studies were isolated using the BD Pharmingen Human Naive/Memory T cell panel (561438; BD) according to the manufacturer’s supplied protocol, followed by FACS. The relative purity of cell populations was assessed prior to RNA-seq. For certain studies as described in the text, TN cells and total TM (TCM + TEM) cells were also isolated from total PBMC using the EasySep Naive CD4+ T cell isolation kit (19555; StemCell Technologies) or the EasySep Memory CD4+ T cell enrichment kit (19157; StemCell Technologies) using the provided protocol. Tissue culture experiments were conducted in RPMI 1640 medium supplemented with 10% FBS, penicillin streptomycin, and l-glutamine. I-BET was obtained from EMD Chemicals (cat no. 401010; Calbiochem) and (±)-JQ1 (SML0974) was obtained from Sigma-Aldrich. The study was approved by the institutional review board at Vanderbilt University Medical Center. Written informed consent was obtained at the time of blood sample collection.
Chromatin immunoprecipitation
Chromatin immunoprecipitation (ChIP) procedures were followed as described using anti–H3K27-Ac Ab (ab6002; Abcam), anti-RNA polymerase II (ab24759; Abcam), anti–NF-κB, p65 subunit (ab47423), anti–T-bet (ab91103; Abcam), or anti-mouse IgG (10 μg; sc-2025; Santa Cruz Biotechnology) (25, 26). Quantitative real-time PCR with SYBR Green was used to measure DNA after ChIP and DNA purification. Primers used to amplify DNA are listed in Supplemental Tables I and II. Quantitative ChIP assays were performed in triplicate and results from three or more independent experiments were averaged. We also performed a modified ChIP assay. Cells were solubilized in nonionic detergents to isolate chromatin. Chromatin was treated with RNase for 10 min, fixed with formaldehyde, sonicated, and processed for ChIP assay using the same procedures described above.
RNA isolation and quantitative RT-PCR
Total RNA was isolated with TRI Reagent (MRC) and purified with the RNeasy MinElute Cleanup kit (Qiagen) using an on-column DNAse treatment to ensure the absence of genomic DNA contamination. cDNA was synthesized from total RNA using the SuperScript III First-Strand Synthesis Kit (Life Technologies) using oligo-dT primers, and purified using the Qiagen QiaQuick PCR purification kit. Transcript levels were measured in duplicate by quantitative RT-PCR using TaqMan assays (ABI 7300 Real Time PCR System; Life Technologies) or using SYBR Green (PowerUp SYBY Green; ABI 7300; Life Technologies). Expression levels were calculated relative to GAPDH using the ΔΔ cycle threshold method. Primer pairs used in SYBR Green reactions are listed in Supplemental Table II.
Immunoprecipitations
Cell transfections
Flow cytometry experiments involving FACS and/or intracellular cytokine measurements small interfering RNA (siRNA) transfection were performed as described previously (26).
RNA-seq sample preparation and data analysis
RNA was isolated using Tri-Reagent. Library preparation was performed using the Illumina Tru-Seq Stranded Total RNA kit. Whole genome RNA sequencing was performed by the Vanderbilt Technologies for Advanced Genomics. An Illumina HiSeq2500 instrument was used to generate 100 bp paired-end reads. A quality control step was performed at all stages of sequencing analysis including raw data, alignment, and expression quantification using tools such as QC3 and MultiRankSEquation (28–30). The RNA data were aligned with TopHat 2 and gene expression levels were quantified using Cufflinks (31–33). Differentially expressed genes measured using Cufflinks/Cuffdiff were expressed as fragments per kilobase per million reads (FPKM) and a cutoff of 0.5 FPKM and 1.0 FPKM was employed for lncRNA and mRNA measurements, respectively. Gene definition lists and results using the GRCh37/hg19 and GRCh38/hg38 were compared and generated similar results. Results from the GRCh37/hg19 builds are presented to permit a more direct comparison with previously described TE and SE locations (21), which were also derived from the GRCh37/hg19 build. False discovery rate (FDR < 0.05) was used to correct for multiple testing. De novo transcriptome assembly was performed on whole genome RNA-seq data from Illumina Tru-Seq Stranded Total RNA libraries with TopHat 2 and Cufflinks using upper quartile normalization (34) and fragment bias correction (35). Novel transcripts were assembled from reads prealigned to human genome 19 (GRCh37/hg19) using TopHat. Identification of novel lincRNA transcripts was accomplished using established methodologies (8, 36). Summarized pipelines for discovery of novel lincRNAs and prediction of transcripts with protein-coding potential were reported previously (26).
For example, to identify TN-specific lncRNAs, we calculated lineage-specific lncRNAs using the equation (FPKM in TN)/[FPKM in (TCM + TEM)/2] and so on for TCM- and TEM-specific memory lncRNAs. We required that each lncRNA exhibit a >5-fold difference between groups. Relative lncRNA expression levels are shown for a portion of these lncRNAs using a heatmap, with each column representing an individual lncRNA and each row representing expression levels of each of three replicate TN, TCM, or TEM cell samples (Fig. 1B). Expression of these lncRNAs clustered according to cell type with TEM cells expressing the greatest number of cell type–specific lncRNAs followed by TN cells and TCM cells.
Cell type–specific expression of mRNAs and lncRNAs. (A) Left panel, total numbers of mRNAs (GRCh38 release, total protein-coding genes = 19,950), annotated lncRNAs (GRCh38 release, long noncoding genes = 15,767), and novel lncRNAs detected in this study by de novo assembly with average FPKM > 1 in at least one cell lineage; middle panel, average expression levels of the different RNA classes in TN, TCM, and TEM cells; right panel, total expression levels of the indicated RNA classes in TN, TCM, and TEM cells (n = 3). (B) Volcano plots identifying lineage-specific RNA classes after FDR correction, p < 0.05. The x-axis represents the indicated expression ratios, log2, and the y-axis is the p value, −log10. (C) Total numbers of mRNAs, annotated lncRNAs, and novel lncRNAs that show lineage-specific expression [from (B)]. (D) Hierarchical clustering showing expression levels of mRNAs expressed by human CD4+ T cells determined by whole genome RNA-seq aligned to hg38 (n = 3). A scale bar is adjacent to the heat map. MSigDB analysis of lineage-specific mRNAs (identified in the heatmap) as input, outputs from the MSigDB analysis are shown with corresponding p values. (E) Individual variation in lncRNA expression. Expression levels in FPKM of four annotated lncRNAs are shown in TN, TCM, and TEM subsets isolated from three individual subjects.
Cell type–specific expression of mRNAs and lncRNAs. (A) Left panel, total numbers of mRNAs (GRCh38 release, total protein-coding genes = 19,950), annotated lncRNAs (GRCh38 release, long noncoding genes = 15,767), and novel lncRNAs detected in this study by de novo assembly with average FPKM > 1 in at least one cell lineage; middle panel, average expression levels of the different RNA classes in TN, TCM, and TEM cells; right panel, total expression levels of the indicated RNA classes in TN, TCM, and TEM cells (n = 3). (B) Volcano plots identifying lineage-specific RNA classes after FDR correction, p < 0.05. The x-axis represents the indicated expression ratios, log2, and the y-axis is the p value, −log10. (C) Total numbers of mRNAs, annotated lncRNAs, and novel lncRNAs that show lineage-specific expression [from (B)]. (D) Hierarchical clustering showing expression levels of mRNAs expressed by human CD4+ T cells determined by whole genome RNA-seq aligned to hg38 (n = 3). A scale bar is adjacent to the heat map. MSigDB analysis of lineage-specific mRNAs (identified in the heatmap) as input, outputs from the MSigDB analysis are shown with corresponding p values. (E) Individual variation in lncRNA expression. Expression levels in FPKM of four annotated lncRNAs are shown in TN, TCM, and TEM subsets isolated from three individual subjects.
Genomic locations of total CD4+, naive CD4, and memory CD4 super enhancers and typical enhancers were previously reported (20).
Statistical analysis
Unless otherwise indicated, data are expressed as the mean ± SD of three or more independent experiments. Unless otherwise noted, significance was determined by Student t test using GraphPad Prism Software (La Jolla, CA). A p value <0.05 was considered significant. χ2 analysis was performed using GraphPad Prism Software. Fisher exact test was performed with Excel; see http://udel.edu/∼mcdonald/statfishers.html.
Data availability
The RNA-seq data generated by this study are accessible through National Center for Biotechnology Information’s Gene Expression Omnibus using the accession code GSE85294 (https://www.ncbi.nlm.nih.gov/geo).
Results
Lineage-specific annotated lncRNAs in naive and memory CD4 T cells
To initiate our studies, we isolated TN, TCM, and TEM cells from age- and gender-matched healthy donors. We performed whole-genome RNA-seq to identify lineage-specific RNAs. We found that ∼12,000 of ∼20,000 mRNAs, ∼2,800 of ∼20,000 annotated lncRNAs, and ∼ 57,000 novel lncRNAs were expressed with means >1 FPKM in at least one T cell lineage (Fig. 1A, left panel). Average expression levels of annotated and novel lncRNAs were comparable and substantially lower than mRNAs (Fig. 1A, middle panel). About 70, 1, and 30% of total RNA detected were represented by mRNA, annotated lncRNA, and novel lncRNA classes, respectively (Fig. 1A, right panel). We also analyzed differences in expression of individual RNA classes in TN, TCM, and TEM lineages after FDR corrections. We found that lineage-specific quantitative differences in expression of novel lncRNAs were greater than lineage-specific differences in expression of either mRNAs or annotated lncRNAs (Fig. 1B). After FDR correction, we found that total numbers of lineage-specific novel lncRNAs were much greater than total numbers of lineage-specific mRNAs or annotated lncRNAs (Fig. 1C). As an initial quality assessment, we further analyzed lineage-specific expression of mRNAs. The heatmap identifies mRNAs selectively expressed in TN, TCM, and TEM lineages (Fig. 1D). As expected (37–39), BACH2 expression was highest in TN cells compared with TCM and TEM cells. In contrast to TN cells, both TCM and TEM cells expressed high levels of IFNG and CCR4. EOMES and genes that encode granzymes (GZMA, GZMH, and GZMK) were highly expressed by the TEM lineage. We identified only a small number of TCM-specific mRNAs compared with TN- and TEM-specific mRNAs. Next, we input our lineage-specific TN, TCM, and TEM gene lists into the Molecular Signature Database (MSigDB, Broad) (40) and the outputs for these lists were TN, TCM, and TEM, respectively. Thus, our isolated cells showed cell type–specific expression of genes encoding lineage-specific cytokines, transcriptional regulators, and effector proteins.
Although the current study was not designed for this purpose, we performed a preliminary analysis to examine lncRNA profiles that may vary due to individual identity. We identified those lncRNAs whose expression varied most according to individual identity and not according to expression in a given TN, TCM, or TEM subset. We found that expression of certain lncRNAs seemed to vary more according to individual identity, independent of expression in a specific T cell subset (Fig. 1E). Thus, individual identity also contributes to expression levels of certain lncRNAs in T cell subsets.
We assessed genomic relationships between genes encoding lineage-specific lncRNAs and genes encoding lineage-specific mRNAs as lncRNAs are known to regulate expression of neighboring protein-coding genes and lncRNA - protein-coding gene pairs exhibit cell type–specific coexpression (10). Examples include IFNG-AS1 and IFNG and TH2LCRR and IL4, IL5, and IL13 (26, 41, 42). We aligned positions of TN-, TCM-, and TEM-specific lncRNA genes with TN-, TCM-, and TEM-specific protein-coding genes across the genome. We used Pearson’s linear regression to determine whether lncRNA:mRNA gene pairs were coexpressed in TN, TCM, and TEM subsets (see Supplemental Table III for sample calculations). We next analyzed a window of up to 250 kb from the transcription start site of each lncRNA gene to determine genomic distance between coexpressed lncRNA:mRNA gene pairs. For this window, we considered the absolute distance of each lncRNA gene from the nearby coexpressed mRNA gene. We found the highest proportion of lncRNA genes were intragenic with coexpressed mRNA genes (Fig. 2A). We also examined this relationship relative to the number of genes falling within this genomic window and obtained similar results (Fig. 2B). Many lncRNA genes were coexpressed with multiple lineage-specific protein-coding genes in the genome. We found that 51 and 49% of naive- and memory-specific lncRNA genes were transcribed in antisense or sense directions relative to adjacent coexpressed protein-coding genes, respectively (Fig. 2C). Thus, the majority of lineage-specific lncRNA genes were coexpressed with lineage-specific protein-coding genes.
Relationships between genomic positions of coexpressed lineage-specific protein-coding genes and annotated lncRNA genes. (A) Results are expressed as the % of lineage-specific mRNA genes (○, TN; ☐, TCM; △, TEM) coexpressed with lineage-specific annotated lncRNA genes relative to total mRNA genes within the indicated distance in kilobase from a lineage-specific annotated lncRNA gene. (B) As in (A) except results are expressed relative to the number of genes in the genome away from the lineage-specific annotated lncRNA gene. (C) Percentage of coexpressed lineage-specific mRNA genes transcribed in antisense or sense orientations relative to adjacent coexpressed lineage-specific annotated lncRNA genes. (D) Genes encoding lineage-specific mRNAs are enriched in the genome near lineage-specific annotated lncRNA genes. The y-axis is the percent of total lineage-specific TN, TCM, and TEM mRNAs (Fig. 1, heatmaps). The x-axis is the distance of the gene encoding a lineage-specific mRNA from the nearest gene encoding a lineage-specific annotated lncRNA. p values determined by χ2 analyses, TN, p = 10−22; TCM, p = 10−12; TEM, p = 10−14.
Relationships between genomic positions of coexpressed lineage-specific protein-coding genes and annotated lncRNA genes. (A) Results are expressed as the % of lineage-specific mRNA genes (○, TN; ☐, TCM; △, TEM) coexpressed with lineage-specific annotated lncRNA genes relative to total mRNA genes within the indicated distance in kilobase from a lineage-specific annotated lncRNA gene. (B) As in (A) except results are expressed relative to the number of genes in the genome away from the lineage-specific annotated lncRNA gene. (C) Percentage of coexpressed lineage-specific mRNA genes transcribed in antisense or sense orientations relative to adjacent coexpressed lineage-specific annotated lncRNA genes. (D) Genes encoding lineage-specific mRNAs are enriched in the genome near lineage-specific annotated lncRNA genes. The y-axis is the percent of total lineage-specific TN, TCM, and TEM mRNAs (Fig. 1, heatmaps). The x-axis is the distance of the gene encoding a lineage-specific mRNA from the nearest gene encoding a lineage-specific annotated lncRNA. p values determined by χ2 analyses, TN, p = 10−22; TCM, p = 10−12; TEM, p = 10−14.
We performed the opposite analysis by determining the frequency that lineage-specific protein-coding genes were near lineage-specific lncRNA genes in the genome to ask whether lineage-specific protein-coding genes were randomly distributed in the genome or were near lineage-specific lncRNA genes. We found that 58, 27, and 33% of genes encoding TN-, TCM-, or TEM-specific mRNAs were within 300 kb of a gene encoding a TN-, TCM-, and TEM-specific lncRNA, respectively (Fig. 2D). We conclude from these studies that lineage-specific annotated lncRNA and mRNA genes are not randomly distributed across the genome but are colocalized in the genome (TN; p = 10−22, TCM; p = 10−12, TEM; p = 10−14, χ2 analysis).
Lineage-specific novel lncRNAs, SEs, and TEs
We expanded our analysis to examine novel or nonannotated lncRNAs expressed in these cell types in a lineage-specific fashion. These transcripts fulfilled criteria of lncRNAs (>200 bp in length and lacked coding potential according to analysis using the coding potential assessment tool, CPAT, and PhyloCSF to estimate evolutionary coding potential). In contrast to the total number of annotated lncRNAs, de novo assembly of our RNA-seq datasets uncovered numerous TN, TCM, and TEM lineage-specific lncRNAs as illustrated in the heatmap (Fig. 3A). We also evaluated the extent to which novel lncRNA genes were coexpressed with adjacent protein-coding genes. For comparison, we determined the distances between mRNA genes selectively expressed in the different populations and found a broad distribution between lineage-specific protein-coding genes in the genome (Fig. 3B, left panel). In contrast, when we examined genomic positions of lineage-specific lncRNA loci and coexpressed lineage-specific protein-coding genes, we found that ∼50% of lineage-specific lncRNA loci were within 0–50 kb of a lineage-specific mRNA gene and another ∼20% were within 50–100 kb of each other (Fig. 3B, right panel). Thus, as with annotated lncRNA loci, lineage-specific novel lncRNA loci exhibited a nonrandom distribution in the genome positioned very near lineage-specific protein-coding genes. We also performed the opposite analysis by determining the frequency at which lineage-specific protein-coding genes were near lineage-specific novel lncRNA genes in the genome to ask whether lineage-specific protein-coding genes were randomly distributed in the genome or were near lineage-specific novel lncRNA genes. We found that 45, 37, and 57% of genes encoding TN-, TCM-, or TEM-specific mRNAs were within 300 kb of a gene encoding a TN-, TCM-, or TEM-specific novel lncRNA, respectively (Fig. 3C). We conclude from these studies that lineage-specific novel lncRNA and mRNA genes are not randomly distributed across the genome but are colocalized in the genome (p = 10−18, 10−14, 10−19 for TN, TCM, and TEM, respectively, χ2 analysis).
Discovery of lineage-specific novel lncRNAs. (A) Hierarchical clustering showing expression levels of novel lncRNAs expressed by human CD4+ TN, TCM, and TEM cells determined by whole genome RNA-seq aligned to hg19 (n = 3). (B) Colocalization of genes encoding lineage-specific mRNAs and novel lncRNAs in the genome. Left panel: % of lineage-specific protein-coding genes that are within the indicated distance in kilobase from another lineage-specific protein-coding gene; right panel: % of lineage-specific protein-coding genes that are within the indicated distance in kilobase from a coexpressed lineage-specific novel lncRNA-coding gene. (C) Genes encoding lineage-specific mRNAs are enriched in the genome near lineage-specific novel lncRNA loci. The y-axis is the percent of total lineage-specific TN, TCM, and TEM mRNAs (Fig. 1, heatmaps). The x-axis is the distance of the gene encoding a lineage-specific mRNA from the nearest loci encoding a lineage-specific novel lncRNA. p values determined by χ2 analyses, TN, p = 10−21; TCM, p = 10−14; TEM p = 10−19.
Discovery of lineage-specific novel lncRNAs. (A) Hierarchical clustering showing expression levels of novel lncRNAs expressed by human CD4+ TN, TCM, and TEM cells determined by whole genome RNA-seq aligned to hg19 (n = 3). (B) Colocalization of genes encoding lineage-specific mRNAs and novel lncRNAs in the genome. Left panel: % of lineage-specific protein-coding genes that are within the indicated distance in kilobase from another lineage-specific protein-coding gene; right panel: % of lineage-specific protein-coding genes that are within the indicated distance in kilobase from a coexpressed lineage-specific novel lncRNA-coding gene. (C) Genes encoding lineage-specific mRNAs are enriched in the genome near lineage-specific novel lncRNA loci. The y-axis is the percent of total lineage-specific TN, TCM, and TEM mRNAs (Fig. 1, heatmaps). The x-axis is the distance of the gene encoding a lineage-specific mRNA from the nearest loci encoding a lineage-specific novel lncRNA. p values determined by χ2 analyses, TN, p = 10−21; TCM, p = 10−14; TEM p = 10−19.
Relationship between novel lncRNAs and transcriptional enhancers
Transcriptional enhancers are known to express RNAs. We performed genome-wide analyses to further explore the relationship among protein-coding gene loci marked by high levels of lineage-specific novel lncRNAs and lineage-specific transcriptional SEs (43) and TEs identified by levels of H3K27-acetylation marks (20) shown in the Circos plot (Fig. 4A). We employed the following scheme: 1) chromosomes with banding patterns and length in megabases in the outer ring; 2) positions of CD4 TN and TM cell TEs and SEs; and 3) positions of TN, TCM, and TEM lineage-specific novel lncRNAs depicted in that order by the three inner rings as a function of expression levels. Although somewhat obscure in the Circos plot, we found a large overlap between genomic positions of SEs and TEs present in the different lineages, lineage-specific novel lncRNA loci, and lineage-specific mRNAs. We next ranked lineage-specific mRNAs according to total number of associated lineage-specific novel lncRNAs (Fig. 4B and Supplemental Table IV), SEs according to total number of associated lineage-specific lncRNAs (Fig. 4C) and TEs according to total number of associated lineage-specific lncRNAs (Fig. 4D). Lineage-specific mRNA genes and lineage-specific SEs and TEs had the highest abundance of associated lineage-specific novel lncRNA loci (linear regression analysis, R2 > 0.9, p < 0.0001). Genes with a high abundance of lineage-specific lncRNA loci and associated SE and TE density encoded key lineage-specific transcriptional regulators (BACH2, KLF12, RORA, GATA3, ATXN1, LEF1), cytokines (IFNG), and receptors (ICAM2, EMB).
Relationships among lineage-specific lncRNAs and transcriptional enhancers. (A) Circos plot depicting genome-wide positions of CD4-naive typical enhancers and super enhancers (purple) along with CD4 memory typical enhancers (orange) and super enhancers (blue) in the outer rings (from 20). Bar graphs colored in black in the inner concentric circles show lineage-specific novel lncRNAs in naive (outermost), central memory, and effector memory (innermost) cells (B) Lineage-specific mRNAs ranked by number of lineage-specific novel lncRNAs per gene loci. Ranks of certain gene loci are noted in the graphs. (C) Lineage-specific SE are ranked by number of novel lncRNAs per genomic locus. (D) Lineage-specific TEs are ranked by number of novel lncRNAs per genomic locus.
Relationships among lineage-specific lncRNAs and transcriptional enhancers. (A) Circos plot depicting genome-wide positions of CD4-naive typical enhancers and super enhancers (purple) along with CD4 memory typical enhancers (orange) and super enhancers (blue) in the outer rings (from 20). Bar graphs colored in black in the inner concentric circles show lineage-specific novel lncRNAs in naive (outermost), central memory, and effector memory (innermost) cells (B) Lineage-specific mRNAs ranked by number of lineage-specific novel lncRNAs per gene loci. Ranks of certain gene loci are noted in the graphs. (C) Lineage-specific SE are ranked by number of novel lncRNAs per genomic locus. (D) Lineage-specific TEs are ranked by number of novel lncRNAs per genomic locus.
We next determined the distances between expressed novel lncRNA loci in the genome. We found that ∼50% of novel lncRNA loci were within 1000 bp of another novel lncRNA loci and >90% were within 10,000 bp of another novel lncRNA genomic loci, suggesting that novel lncRNA loci are clustered in the genome (Fig. 5A). We refer to these small genomic loci that transcribe dense numbers of novel lncRNAs as bins (Fig. 5B). We found that the vast majority of lineage-specific novel lncRNAs were transcribed from these bins that had an average size of <30 kb. Distances between bins were much greater, averaging ∼600 kb, indicating the nonrandom distribution of lncRNA loci in the genome. In large part, expression of all the lncRNAs transcribed from an individual bin exhibited lineage-specific expression patterns. Examples include a bin selectively expressed in TN (Fig. 5C), TCM (Fig. 5D), and TEM (Fig. 5E), respectively.
Genomic distributions of novel lncRNAs. (A) Distances between genomic loci transcribing unique novel lncRNAs. The x-axis ranks novel lncRNA loci according to distance to the next novel lncRNA loci and the y-axis shows the actual distance in base pair between neighboring lncRNA loci. (B) Characteristics of genomic bins transcribing multiple novel lncRNAs. (C–E) Examples of genomic bins selectively transcribing multiple novel lncRNAs in (C) TN relative to TM, (D) TCM relative to TN and TEM, and (E) TEM relative to TN and TCM. The x-axes show genomic positions from which novel lncRNAs are transcribed; lines illustrate sizes of the transcribed novel lncRNA. The y-axes are mean expression ratios, log2, between the indicated T cell lineages. (F) Genomic bins transcribing novel lncRNAs colocalize with TN and TM typical and super transcriptional enhancers (from 20). We considered that the genomic bins may not exactly overlap with enhancers so we performed the analysis to include the proportion of bins that colocalize with an enhancer (0), and the proportion of bins that colocalized with enhancers if we extended bin size by 5, 10, 15, 20 or 30 kb in both 5′ and 3′ directions, x-axis. The y-axis is the proportion of total bins with a TN or TM enhancer. The p values, −log10, were determined by χ2 analysis. (G) Number of bins within TM SE or TE (CD4MSE, CD4MTE) structures or TN SE or TE structures. The x-axis is as in (F), y-axis is the number of bins with enhancers. (H and I) Reanalysis of RNA-seq data using bins as definition lists. (H) After FDR correction, bins were identified that were overexpressed (blue lines) or underexpressed (brown lines) in TCM and TEM cells compared with TN cells. (I) Expression levels of bins differentially expressed in TCM relative to TN and TEM after FDR correction. The y-axes are expression ratios, log2; x-axes identify ratios; Tn = TN/TN, Tcm = TCM/TN, Tem = TEM/TN.
Genomic distributions of novel lncRNAs. (A) Distances between genomic loci transcribing unique novel lncRNAs. The x-axis ranks novel lncRNA loci according to distance to the next novel lncRNA loci and the y-axis shows the actual distance in base pair between neighboring lncRNA loci. (B) Characteristics of genomic bins transcribing multiple novel lncRNAs. (C–E) Examples of genomic bins selectively transcribing multiple novel lncRNAs in (C) TN relative to TM, (D) TCM relative to TN and TEM, and (E) TEM relative to TN and TCM. The x-axes show genomic positions from which novel lncRNAs are transcribed; lines illustrate sizes of the transcribed novel lncRNA. The y-axes are mean expression ratios, log2, between the indicated T cell lineages. (F) Genomic bins transcribing novel lncRNAs colocalize with TN and TM typical and super transcriptional enhancers (from 20). We considered that the genomic bins may not exactly overlap with enhancers so we performed the analysis to include the proportion of bins that colocalize with an enhancer (0), and the proportion of bins that colocalized with enhancers if we extended bin size by 5, 10, 15, 20 or 30 kb in both 5′ and 3′ directions, x-axis. The y-axis is the proportion of total bins with a TN or TM enhancer. The p values, −log10, were determined by χ2 analysis. (G) Number of bins within TM SE or TE (CD4MSE, CD4MTE) structures or TN SE or TE structures. The x-axis is as in (F), y-axis is the number of bins with enhancers. (H and I) Reanalysis of RNA-seq data using bins as definition lists. (H) After FDR correction, bins were identified that were overexpressed (blue lines) or underexpressed (brown lines) in TCM and TEM cells compared with TN cells. (I) Expression levels of bins differentially expressed in TCM relative to TN and TEM after FDR correction. The y-axes are expression ratios, log2; x-axes identify ratios; Tn = TN/TN, Tcm = TCM/TN, Tem = TEM/TN.
Given these data, we further analyzed the relationship between genomic positions of TN and TM cell enhancers and TN and TM cell bins, and determined the proportion of lineage-specific bins that contained lineage-specific enhancers. We found that lineage-specific bins were highly enriched in the lineage-specific SE and TE portion of the genome (Fig. 5F, 5G). In fact, almost 50% of bins contained at least one lineage-specific SE or TE. Considering that bin loci may not exactly overlap with TE and SE loci, we asked whether expanding bin size in either 5′ or 3′ directions increased the fraction of bins with TEs or SEs and found that increasing bin size by 20–30 kb increased the proportion of bins with TEs and SEs to ∼80%, which was also highly significant. We conclude that positions of lineage-specific bin loci in the genome overlap lineage-specific TE and SE loci.
We reanalyzed the RNA-seq data using bin loci as definition lists, and determined total FPKM expressed from a given bin in TN, TCM, and TEM lineages to determine whether bin loci were expressed in selective lineages. After FDR correction, we found that bins could be divided into three categories (Fig. 5H, 5I). About half the bins were highly expressed in TN lineages compared with TCM and TEM lineages and about half were highly expressed in TCM and TEM lineages compared with the TN lineage. We identified a small number of bins that were differentially expressed in the TCM compared with TN and TEM lineages. We conclude that lncRNAs expressed from these bins also exhibit lineage-specific expression.
IFNG-associated known and novel lncRNAs
We selected two genes for more detailed investigation: IFNG because it is not only selectively expressed by effector Th1 cells and silenced by effector Th2 and Th17 cells but is also expressed by TEM cells after TCR stimulation, and BACH2 because it is selectively expressed in TN cells and silenced after T cell activation. We found that the IFNG gene locus expressed a high density of novel lncRNAs in TEM and TCM cells (Fig. 6A). For IFNG, we examined RNA-seq tracks (∼300 kb) not only in TN, TCM, and TEM cells but also in primary and effector Th1, Th2, and Th17 cells (primary, 3 d culture under respective polarizing conditions; effector, 5–7 d culture under the respective polarizing conditions followed by restimulation with anti-CD3 and additional 2 d culture). Similar to TN cells, we found a paucity of lncRNAs transcribed by this region in Th2 and Th17 cells (Fig. 6A). In contrast, lncRNAs were highly transcribed from this region by Th1 cells as was IFNG, with levels of lncRNAs somewhat higher in effector than primary Th1 cells. LncRNAs transcribed from the IFNG locus were similarly elevated in TCM and TEM cells, and somewhat higher in TEM than TCM cells even though IFNG was not actively transcribed by resting TEM and TCM cells. Our de novo assembly also identified lncRNAs with defined sizes and genomic locations and these lncRNAs and genomic positions are identified below the RNA-seq tracks. Positions of TEs and SEs based on levels of H3K27-Ac modification in TM cells are also shown. We confirmed by RT-PCR that these novel lncRNAs were selectively expressed in TCM and TEM cells relative to TN cells (Fig. 6B). Each lncRNA transcribed from the IFNG locus was named relative to the IFNG transcriptional start site in the human genome. For example, for IFNG-R-163, the R is for RNA and -163 means it is transcribed from a position 163 kb from IFNG in the direction of the P terminus. RT-PCR was performed using oligo-dT for cDNA synthesis, so our conclusion is that these IFNG-R-lncRNAs are poly-adenylated.
Novel lncRNAs at the IFNG locus. (A) RNA-seq tracks spanning ∼300 kb from IFNG-AS1 past IL-26 in primary and effector Th1, Th2 and Th17 cells, naive CD4+ T cells and TCM and TEM cells. Below are positions of genomic loci encoding novel lncRNAs determined by de novo assembly, genomic positions of typical and super enhancers (TE and SE, respectively) identified in total CD4+ TM cells are also shown (from Ref. 20). (B) PCR validation of RNA-seq results normalized to GAPDH. Each novel lncRNA is named relative to the IFNG transcriptional start site, − indicates toward the p-terminus, + toward the q-terminus, differences between either TCM and TN or TEM and TN were statistically significant, n = 4, p < 0.05. (C) JQ1 treatment of memory T cells abrogates RNA polymerase II binding but not H3K27Ac across the IFNG locus. Indicated concentrations of JQ1 were added to TM cultures. Cells were stimulated with plate-bound anti-CD3 for 24 h and levels of H3K27-Ac (left panel) and RNA polymerase II (right panel) were determined by ChIP. For this analysis, we designed a tiling array. A total of 100 PCR primer pairs were designed at ∼3 kb intervals to span the 300-kb IFNG locus. Data are expressed as mean fraction of input (n = 3) p < 0.05, ANOVA. (D) Total memory CD4+ T cells were cultured for 24 h with JQ1 or I-BET to displace BRD-containing proteins from genomic loci. RNA was purified and levels of the indicated IFNG locus novel lncRNAs, IFNG-R-#, determined by PCR. Results are expressed relative to levels of GAPDH. p < 0.05, treated versus untreated (n = 4). (E) JQ1 or I-BET inhibit IFNG-AS1 and IFNG expression. Total memory CD4+ T cells were treated with and without JQ1 or I-BET for 24 h and stimulated with anti-CD3 for an additional 24 h. Levels of three unique isoforms of IFNG-AS1 designated 1, 2, 3 (A), and IFNG were determined by PCR and normalized to levels of GAPDH. *p < 0.05 (n = 4). (F) Individual siRNAs targeting IFNG-AS1 or IFNG-R-49 were transfected into total memory CD4+ T cells prior to stimulation with anti-CD3 for 24 h. Transfected cells were isolated by FACS prior to mRNA analysis or to measure intracellular levels of IFN-γ. Error bars are mean ± SD. The p values were determined by unpaired t test. *p < 0.05.
Novel lncRNAs at the IFNG locus. (A) RNA-seq tracks spanning ∼300 kb from IFNG-AS1 past IL-26 in primary and effector Th1, Th2 and Th17 cells, naive CD4+ T cells and TCM and TEM cells. Below are positions of genomic loci encoding novel lncRNAs determined by de novo assembly, genomic positions of typical and super enhancers (TE and SE, respectively) identified in total CD4+ TM cells are also shown (from Ref. 20). (B) PCR validation of RNA-seq results normalized to GAPDH. Each novel lncRNA is named relative to the IFNG transcriptional start site, − indicates toward the p-terminus, + toward the q-terminus, differences between either TCM and TN or TEM and TN were statistically significant, n = 4, p < 0.05. (C) JQ1 treatment of memory T cells abrogates RNA polymerase II binding but not H3K27Ac across the IFNG locus. Indicated concentrations of JQ1 were added to TM cultures. Cells were stimulated with plate-bound anti-CD3 for 24 h and levels of H3K27-Ac (left panel) and RNA polymerase II (right panel) were determined by ChIP. For this analysis, we designed a tiling array. A total of 100 PCR primer pairs were designed at ∼3 kb intervals to span the 300-kb IFNG locus. Data are expressed as mean fraction of input (n = 3) p < 0.05, ANOVA. (D) Total memory CD4+ T cells were cultured for 24 h with JQ1 or I-BET to displace BRD-containing proteins from genomic loci. RNA was purified and levels of the indicated IFNG locus novel lncRNAs, IFNG-R-#, determined by PCR. Results are expressed relative to levels of GAPDH. p < 0.05, treated versus untreated (n = 4). (E) JQ1 or I-BET inhibit IFNG-AS1 and IFNG expression. Total memory CD4+ T cells were treated with and without JQ1 or I-BET for 24 h and stimulated with anti-CD3 for an additional 24 h. Levels of three unique isoforms of IFNG-AS1 designated 1, 2, 3 (A), and IFNG were determined by PCR and normalized to levels of GAPDH. *p < 0.05 (n = 4). (F) Individual siRNAs targeting IFNG-AS1 or IFNG-R-49 were transfected into total memory CD4+ T cells prior to stimulation with anti-CD3 for 24 h. Transfected cells were isolated by FACS prior to mRNA analysis or to measure intracellular levels of IFN-γ. Error bars are mean ± SD. The p values were determined by unpaired t test. *p < 0.05.
TEs and SEs can be defined by levels of H3K27Ac. JQ1 and I-BET are cell-permeable small molecules that displace BRD-containing proteins from acetylated lysine motifs and disrupt function of TEs and SEs but do not directly reverse H3K27Ac marks (35, 44, 45). We also found that the IFNG locus spanning ∼300 kb (IFNG-AS1-IFNG-IL26) was marked by H3K27Ac in TM cells, which was not reversed by treatment with JQ1 (Fig. 6C, left panel). In contrast, RNA pol2 was also recruited to the IFNG locus in TM cells and RNA pol2 was displaced by treatment with JQ1 (Fig. 6C, right panel). To ask whether lineage-specific novel lncRNA expression at these gene loci was linked to TE or SE function, we studied whether displacement of BRD proteins by JQ1 or I-BET affected levels of novel lncRNAs transcribed from the IFNG locus. We found that both JQ1 and I-BET effectively inhibited expression of these lncRNAs in TM cells (Fig. 6D). Thus, these novel lncRNA transcripts are consistent with the definition of enhancer-associated RNAs or eRNAs.
Three unique isoforms of the known lncRNA, IFNG-AS1, exist in current human genome annotations. Levels of each of these transcripts were also reduced by treatment with JQ1 or I-BET (Fig. 6E). Although IFNG is expressed at low levels in resting TM cells, it is highly expressed after stimulation via the TCR. We also found that pretreatment with JQ1 or I-BET abrogated induction of IFNG expression in TM cells (Fig. 6E). Transfection of an siRNA targeting either IFNG-AS1 or IFNG-R-49 led to marked reduction of both IFNG mRNA and protein expression (Fig. 6F). Taken together, these data demonstrate that both IFNG-AS1 and the enhancer-associated RNA, IFNG-R-49, are required for efficient induction of IFNG by TM cells.
NF-κB binds novel lncRNAs transcribed from the IFNG locus
Recent studies demonstrate that RNAs transcribed from enhancers bind the transcription factor, YY1, as well as the histone acetyltransferase, CBP, to enhance function of these proteins and stimulate transcription of target genes (46, 47). The prosurvival, proinflammatory transcription factor, NF-κB, as well as T-bet, have multiple binding sites across the IFNG locus and both are critical for IFNG expression (48). Thus, we asked whether one function of the IFNG-R lncRNAs was to bind either T-bet or NF-κB. To explore this idea, we performed immunoprecipitations with either anti–T-bet or anti–NF-κB, p65 subunit, and measured IFNG-associated lncRNAs in the immunoprecipitates. We failed to detect IFNG-associated lncRNAs in the anti–T-bet immunoprecipitates but found that many IFNG-associated lncRNAs were detected in the anti-NF-κB immunoprecipitates indicating a portion of these enhancer-associated lncRNAs specifically associate with NF-κB (Fig. 7A, 7B). As a control experiment, we depleted cells of enhancer-associated lncRNAs by JQ1 treatment and found this treatment also resulted in loss of IFNG-associated lncRNAs from the anti-NF-κB immunoprecipitates (Fig. 7B). More importantly, using ChIP assays, we found that treatment of cells with JQ1 resulted in loss of NF-κB binding to chromatin across the IFNG locus (Fig. 7C). Next, we reasoned that if IFNG enhancer-associated lncRNAs contributed to affinity of NF-κB for IFNG locus genomic sites, then treatment of chromatin with RNase may reduce binding of NF-κB to IFNG genomic sites. To explore this possibility, we isolated chromatin without formaldehyde crosslinking, treated chromatin with RNase followed by formaldehyde crosslinking, and processed chromatin for ChIP assays. We found that treatment of chromatin with RNase resulted in loss of NF-κB binding to IFNG genomic sites (Fig. 7D). We also treated isolated chromatin from TM cells with RNAse and found that this treatment resulted in release of substantial NF-κB protein from chromatin (Fig. 7E). Finally, we transfected cells with the siRNA specific for IFNG-R-49 or a scrambled control siRNA. Transfection of the siRNA specific for IFNG-R-49 reduced NF-κB binding to the IFNG-D-49 genomic DNA site but not several other sites across the IFNG locus (Fig. 7F). Our conclusion is that IFNG-enhancer associated lncRNAs contribute to maintaining binding of NF-κB to the IFNG locus and in general chromatin-associated RNAs contribute to localization of NF-κB to chromatin, suggesting this notion may have general applicability.
IFNG-associated lncRNAs bind NF-κB. (A and B) Cells were cultured with or without JQ1 (150 nM, 4 h) and lysed. After lysis, (A) anti–T-bet and (B) anti–NF-κB, p65 subunit, immunoprecipitations (IP) were performed. Levels of each IFNG-R-lncRNA in the anti–T-bet, anti–NF-κB, p65 subunit, and isotype control immunoprecipitates were determined by PCR. Results are expressed as fraction of input of the indicated IFNG-associated novel lncRNAs relative to totals of each IFNG-R-lncRNA. Levels of IFNG-R-associated lncRNAs in the isotype control immunoprecipitates were negligible and therefore not included in the calculations. Error bars are SD of eight independent experiments, *p < 0.05. (C) As in (B) except ChIP assays were performed to measure NF-κB, p65, binding to IFNG genomic locus sites. Results are expressed as a fraction of input relative to an isotype control. Error bars are SD of five independent experiments, *p < 0.05. (D) Chromatin was isolated from TM cells by lysis in nonionic detergents and treated with RNase, 10 min, room temperature. Chromatin was cross-linked with formaldehyde and processed for ChIP assays. NF-κB, p65 binding is expressed as fraction of input, n = 3, *p < 0.05 comparing treated and untreated samples. (E) After RNAse treatment, chromatin was pelleted by centrifugation and supernatant fluids harvested and analyzed for NF-κB protein by Western blotting. The arrow indicates the p65 subunit of NF-κB. (F) Cells were transfected with the siRNA targeting IFNG-R-49 (+) or a scrambled control siRNA (−) and cultured for 24 h. Cultures were fixed with paraformaldehyde, harvested, and processed for ChIP assays using an NF-κB, p65 subunit, Ab. The y-axis shows NF-κB, p65 subunit, binding as fraction of input at the indicated IFNG genomic loci, x-axis.
IFNG-associated lncRNAs bind NF-κB. (A and B) Cells were cultured with or without JQ1 (150 nM, 4 h) and lysed. After lysis, (A) anti–T-bet and (B) anti–NF-κB, p65 subunit, immunoprecipitations (IP) were performed. Levels of each IFNG-R-lncRNA in the anti–T-bet, anti–NF-κB, p65 subunit, and isotype control immunoprecipitates were determined by PCR. Results are expressed as fraction of input of the indicated IFNG-associated novel lncRNAs relative to totals of each IFNG-R-lncRNA. Levels of IFNG-R-associated lncRNAs in the isotype control immunoprecipitates were negligible and therefore not included in the calculations. Error bars are SD of eight independent experiments, *p < 0.05. (C) As in (B) except ChIP assays were performed to measure NF-κB, p65, binding to IFNG genomic locus sites. Results are expressed as a fraction of input relative to an isotype control. Error bars are SD of five independent experiments, *p < 0.05. (D) Chromatin was isolated from TM cells by lysis in nonionic detergents and treated with RNase, 10 min, room temperature. Chromatin was cross-linked with formaldehyde and processed for ChIP assays. NF-κB, p65 binding is expressed as fraction of input, n = 3, *p < 0.05 comparing treated and untreated samples. (E) After RNAse treatment, chromatin was pelleted by centrifugation and supernatant fluids harvested and analyzed for NF-κB protein by Western blotting. The arrow indicates the p65 subunit of NF-κB. (F) Cells were transfected with the siRNA targeting IFNG-R-49 (+) or a scrambled control siRNA (−) and cultured for 24 h. Cultures were fixed with paraformaldehyde, harvested, and processed for ChIP assays using an NF-κB, p65 subunit, Ab. The y-axis shows NF-κB, p65 subunit, binding as fraction of input at the indicated IFNG genomic loci, x-axis.
Lineage-specific novel lncRNAs at the BACH2 locus
We used a similar strategy to analyze novel lncRNAs transcribed from the BACH2 locus in TN and TM lineages. RNA-seq tracks identified high levels of transcripts across the BACH2 gene in TN cells that were very reduced in TCM and TEM cells (Fig. 8A). Our de novo assembly defined >13 unique lncRNAs transcripts preferentially expressed in TN compared with TCM and TEM cells (Fig. 8B). As above, we confirmed these expression differences by PCR and obtained similar results (Fig. 8C). We treated TN cells with JQ1 and I-BET and found these inhibitors substantially lowered levels of BACH2-associated lncRNAs (Fig. 8D) and BACH2 mRNA (Fig. 8E). Thus, inhibitors that displace BRD proteins from H3K27Ac motifs reduced levels of both BACH2 enhancer–associated novel lncRNAs and BACH2 mRNA. Overall, we find these results consistent with the notion that genomic loci surrounding lineage-specific protein-coding genes are enriched with lineage-specific TEs and SEs that transcribe cell type–specific novel lncRNAs.
Novel lncRNAs at the BACH2 locus. (A) RNA sequencing tracks of the BACH2 gene in TN, TCM, and TEM cells. Positions of known TEs and SEs in total CD4+ T cells are shown below the sequencing tracks (from 20). Exon-intron structure of BACH2 is the bottom track. (B) Expression levels of the 13 BACH2-associated novel lncRNAs discovered by de novo assembly in TN, and TEM, and TCM expressed as average FPKM, n = 3, *p < 0.05. (C) PCR verification of results in (B) expressed relative to GAPDH in TN and TM cells, n = 4, *p < 0.05. (D) TN cells were cultured for 24 h with JQ1 or I-BET to displace BRD-containing proteins from genomic loci. RNA was purified and levels of 13 BACH2 associated novel lncRNAs were determined by PCR. Results are expressed relative to levels of GAPDH. n = 4, *p < 0.05. (E) As in (D), except levels of BACH2 mRNA were determined by PCR and are expressed relative to GAPDH. The p values were determined using an unpaired t test. Error bars are mean ± SD, *p < 0.05.
Novel lncRNAs at the BACH2 locus. (A) RNA sequencing tracks of the BACH2 gene in TN, TCM, and TEM cells. Positions of known TEs and SEs in total CD4+ T cells are shown below the sequencing tracks (from 20). Exon-intron structure of BACH2 is the bottom track. (B) Expression levels of the 13 BACH2-associated novel lncRNAs discovered by de novo assembly in TN, and TEM, and TCM expressed as average FPKM, n = 3, *p < 0.05. (C) PCR verification of results in (B) expressed relative to GAPDH in TN and TM cells, n = 4, *p < 0.05. (D) TN cells were cultured for 24 h with JQ1 or I-BET to displace BRD-containing proteins from genomic loci. RNA was purified and levels of 13 BACH2 associated novel lncRNAs were determined by PCR. Results are expressed relative to levels of GAPDH. n = 4, *p < 0.05. (E) As in (D), except levels of BACH2 mRNA were determined by PCR and are expressed relative to GAPDH. The p values were determined using an unpaired t test. Error bars are mean ± SD, *p < 0.05.
Discussion
The purposes of the investigations presented in this study were to 1) define differential expression patterns of mRNAs, annotated, and novel lncRNAs in human TN, TCM, and TEM cells; 2) explore relationships among their expression patterns; and 3) explore relationships among SEs, TEs, and novel lncRNAs. We found that most lineage-specific annotated lncRNAs are coexpressed with adjacent lineage-specific protein-coding genes in the genome, suggesting these annotated lncRNAs contribute to regulation of lineage-specific protein-coding genes and help imprint these very specific cell phenotypes. Previous studies have also defined lineage-specific expression patterns of both annotated and novel lncRNAs in T cells (26, 49, 50). Our de novo assembly identifies many more lineage-specific novel lncRNAs (∼15,000) compared with lineage-specific annotated lncRNAs (∼200). Expression levels of annotated and novel lncRNAs are comparable and these novel lncRNA loci are also colocalized and coexpressed with lineage-specific protein-coding genes, suggesting they also contribute to regulation of lineage-specific protein-coding genes and underlying cellular phenotypes. Many known lncRNA genes are composed of exons and introns and RNAs are spliced and polyadenylated to yield a mature lncRNA (10). Our de novo assembly reveals that most novel lncRNAs we discovered are polyadenylated but do not undergo splicing to yield mature lncRNAs. This probably explains why we uncovered many more novel lncRNAs compared with annotated lncRNAs described in previous reports (29, 36, 37). These novel lncRNAs are transcribed from enhancer regions found in hematopoietic cells, TEs, and SEs, and inhibition of enhancer function abrogates novel lncRNA expression, suggesting they more closely resemble the class of 1d-eRNAs (17). Lineage-specific expression patterns indicate they are biologically regulated. Transcription of these novel lncRNAs is prevented by inhibition of enhancer function; as is transcription of adjacent lineage-specific protein-coding genes, IFNG and BACH2, further supporting the notion that their transcription is derived from active enhancers.
Additional roles and mechanisms of regulation of these novel lncRNAs are suggested by our studies. The IFNG genomic locus possesses enhancers specific for IFNG and for IFNG-AS1, which also regulates IFNG expression (48). Both groups of enhancers transcribe novel lncRNAs and are marked by high levels of H3K27Ac and RNA pol2. Both groups of novel lncRNAs are highly regulated during Th cell differentiation in which IFNG is actively expressed or actively silenced. However, in resting TEM cells, where IFNG is not highly expressed, novel lncRNAs are still present at high levels as is the known lncRNA, IFNG-AS1. Our results also demonstrate that continued presence of IFNG-enhancer–associated novel lncRNAs and IFNG-AS1 is necessary for rapid expression of IFNG by TEM cells in response to TCR stimulation. IFNG-AS1 is known to be required for efficient expression of IFNG by murine effector CD4+ and CD8+ lymphocytes so it may not seem unexpected that IFNG-AS1 is also required for efficient expression of IFNG by human TEM cells (41, 42). However, our results extend the function of IFNG-AS1 to TEM cells. In contrast, it might seem somewhat unexpected that one novel lncRNA, IFNG-R-49, transcribed from the IFNG locus, is required for efficient expression of IFNG by TEM cells when numerous novel lncRNAs are transcribed from this locus. However, it is worth noting that multiple transcriptional enhancers have been identified across this locus and several examples exist where deletion of a single 1 kb enhancer is sufficient to abrogate IFNG expression by Th1 effector cells and/or TM cells (51–53). Novel IFNG-associated lncRNAs appear to act, in part, to increase binding of NF-κB to the IFNG genomic locus. In separate studies, enhancer-associated lncRNAs have been shown to bind the transcription factor YY1 and reinforce binding of YY1 to chromatin and bind the histone acetyltransferase, CBP, to stimulate histone acetylation and transcription of target genes (46, 47). Thus, it would appear that one common mechanism by which these novel lncRNAs function at other genomic loci is to reinforce transcription factor and histone acetyltransferase binding to chromatin and function.
One question not completely addressed in this study is the functional relationship between IFNG-AS1 and enhancer-associated RNAs transcribed from the IFNG locus. Depletion of either IFNG-AS1 or IFNG-R-49 is sufficient to reduce IFNG transcription by TM cells. Depletion of IFNG-R-49 also reduces NF-κB binding at IFNG-D-49 but not other enhancer loci. Previous studies suggest that there are distinct NF-κB response enhancers required for IFNG-AS1 expression and distinct NF-κB response enhancers required for IFNG expression and both groups of enhancers contribute to expression of IFNG by targeting either IFNG-AS1 or IFNG genes (48, 54). Two alternate hypotheses that cannot be excluded are 1) that the major role of NF-κB is to drive IFNG-AS1 expression and IFNG-AS1, in turn, drives IFNG expression or 2) that the major role of IFNG-AS1 is to drive expression of the IFNG enhancer associated lncRNAs that, in turn reinforce NF-κB binding to IFNG enhancer elements and stimulate IFNG expression. Additional studies will be required to distinguish among these hypotheses.
Conversely, the transcriptional regulator, BACH2, is highly expressed in TN cells, represses effector Th programs, and is actively silenced in response to stimulation of Th cell differentiation, and silencing is sustained in TCM and TEM cells (38). The BACH2 gene locus transcribes a dense pattern of >80 novel lncRNAs in TN cells and these novel lncRNAs are silenced in TCM and TEM cells. Inhibition of enhancer function in TN cells results in the loss of BACH2-associated lncRNAs and BACH2 mRNA, linking these novel lncRNAs to regulation of BACH2. Thus, via generation of these novel associated lncRNAs, enhancers may positively regulate expression of protein-coding genes to allow commitment along various cell differentiation paths.
The general argument has been made that genes critically important for cell identity and function, such as BACH2 and IFNG, are marked by SEs and TEs and, by extension, these novel enhancer-associated lncRNAs (21). The converse argument, that loci marked by high densities of SEs, TEs, and enhancer-associated lncRNAs may identify genes critical to maintain stability and function of a given cell lineage, such as TN, TCM, and TEM lineages, may also be correct. Therefore, it may be possible to identify new protein-coding genes critical for TN, TCM, and TEM stability and function by ranking them according to SE, TE, and novel lncRNA load as we have done in this study. For example, KLF12 is marked by high levels of lncRNAs and SEs and TEs in TM cells yet the role of this transcription factor in TM cell function and stability is poorly understood.
We understand a number of basic mechanisms that annotated lncRNAs employ to achieve their diverse regulatory functions. In contrast, our understanding of biologic functions of enhancer-associated RNAs is more limited (17, 18, 55). One function of RNAs transcribed from TEs and SEs is that they bind the same transcription factors that bind the DNA elements from which the RNA is transcribed, and these eRNAs help tether transcription factors to these regulatory DNA elements, which in turn stabilizes gene-expression programs and reinforces cell phenotypes (46). Our studies are consistent with this type of mechanism. Other studies suggest that eRNAs may contribute to chromosome conformation and nuclear localization (56). Additional studies will be required to fully understand the diverse functions of this class of RNAs, not only in immunological memory, but also how they contribute to biologic complexity in general. Nevertheless, our results support the notion that both annotated lncRNAs and enhancer-associated lncRNAs play critical roles in maintaining the identity of these key T cell lineages, enabling these cells to carry out their unique adaptive immune functions.
A major difference between TN and TM cells is that TN cells require TCR signaling and additional stimuli, such as cytokine receptor stimulation, to induce or activate a number of transcription factors that initiate differentiation paths, making these cells competent to produce key inflammatory cytokines, chemokines, and other proteins required by the adaptive immune response, whereas TM cells only require TCR signaling. This is also revealed by the transcription factor requirements to induce these genes. For example, Stat4, T-bet, NF-κB, and other transcription factors play key roles in the induction of IFNG by developing Th1 cells whereas rapid IFNG expression by TM cells is relatively Stat4 and T-bet independent but still NF-κB dependent, and TCR signaling alone is sufficient to activate NF-κB (6, 57). We propose that expression of enhancer-associated lncRNAs enables NF-κB recruitment to the IFNG locus by T-bet-/Stat4-independent mechanisms driving rapid expression of IFNG by TM cells, thus reducing transcription factor requirements for expression of this key gene of the memory response. This underlying mechanism may also contribute to sustained and inducible expression of other key protein-coding genes, allowing TM cells to carry out their unique functions providing lifelong immunity against pathogens, a key element of the adaptive immune response.
Footnotes
This work was supported by grants from the National Institutes of Health (R01 AI044924, DK007383, and DK20593), the National Science Foundation (DGE0909667), and by Clinical and Translational Science Award UL1TR000445. Vanderbilt Technologies for Advanced Genomics’ core facility was supported in part by grants from the National Institutes of Health (P30 CA68485, P30 EY08126 and G20 RR030956).
The online version of this article contains supplemental material.
Abbreviations used in this article:
- BRD
bromodomain-containing protein
- ChIP
chromatin immunoprecipitation
- eRNA
enhancer RNA
- FDR
false discovery rate
- FPKM
fragments per kilobase per million reads
- H3K27Ac
H3K27-acetylation
- lncRNA
long noncoding RNA
- MSigDB
Molecular Signature Database
- ncRNA
noncoding RNA
- RNA-seq
RNA sequencing
- SE
superenhancer
- siRNA
small interfering RNA
- TCM
central memory T
- TE
typical enhancer
- TEM
effector memory T
- TM
memory T
- TN
naive T.
References
Disclosures
The authors have no financial conflicts of interest.