Long noncoding RNAs in B-cell development and activation

Long noncoding RNAs (lncRNAs) have emerging roles in innate and adaptive immunity. For example, lnc-DC is required for normal dendritic cell differentiation and function,¹ IL1β-eRNA and IL1β-RBT46 are required for lipopolysaccharide-induced pro-inflammatory responses in monocytes,²  and NRAV modulates cellular responses to viral infections.³  In T cells, an intronic lncRNA NRON abrogates the nuclear transport of nuclear factor of activated T cells, and hence modulates expression of interleukin-2.⁴  In B-cell lymphomas, the lncRNA Fas-AS1 modulates expression of soluble Fas receptor messenger RNA, an important regulator of apoptosis.⁵  Thus, lncRNAs have the potential to influence both normal and pathological immune cell development and function.

LncRNAs may operate via a variety of molecular mechanisms.⁶  For example, enhancer-associated lncRNAs (eRNAs) act in cis and originate from transcribed extragenic or intragenic enhancer regions, whereas promoter-associated lncRNAs (pRNAs) can act in trans and originate from canonical promoter-derived transcriptional activity.^7,8  These 2 broad lncRNA categories are distinguished by the ratio of mono- vs tri-methylation of histone 3 lysine 4 (H3K4me1/H3K4me3).⁸  Compared with pRNAs, eRNAs tend to exhibit more restricted expression and their RNA sequences show less constraint.⁸ 

Advances in sequencing technology have enabled the identification of large numbers of putative lncRNA loci.^9,10  However, the proportion of lncRNAs with clearly defined function is small,^11,12  caused in part by poor annotation of lncRNAs expressed in a tissue of interest, making it difficult to select candidate lncRNAs for targeted studies. This is a consequence of the expression patterns of lncRNAs, which are often restricted to 1 or very few tissues or cell types.⁹  Recent studies have addressed this limitation by surveying lncRNA expression in several organisms and tissues,^13-17  including murine T cells¹⁸ ; however, there have been no comparable attempts to use sequencing technologies to describe the murine B-cell lncRNA repertoire.

To facilitate the study of lncRNA biology in B cells, we describe a catalog of 4516 de novo assembled high-confidence lncRNAs expressed in 11 mouse B-cell populations. We identify human lncRNAs that may be orthologs of the mouse genes. Furthermore, we classify subsets of eRNAs and pRNAs and perform an unsupervised clustering analysis to associate lncRNAs with messenger RNAs at key stages of B-cell development. Finally, we use the lncRNA catalog to show that PAX5, a transcription factor required to specify the B-cell lineage,¹⁹  binds to and regulates expression of lncRNA loci in both pro-B and mature B cells as well as in acute lymphoblastic leukemia.

All RNA-seq and chromatin immunoprecipitation (ChIP)-seq experiments were performed with female C57BL/6JNimr mice aged 7-9 weeks, except for RNA-seq of plasmablasts and plasma cells, which were obtained from 12- to 14-week-old Blimp1-GFP mice.²⁰ 

Gating strategies for cell sorting are shown in supplemental Figure 1, available on the Blood Web site.

Sorted populations of cells were resuspended in Trizol (Life Technologies), and RNA was purified using the RNeasy Mini Kit (Qiagen). RNA quality was assessed using the 2100 Expert Agilent Bioanalyzer. For all samples except plasmablasts and plasma cells, stranded polyA-enriched libraries were made using the Stranded TruSeq RNA Sample Preparation Kit (Illumina) and sequenced on the HiSeq 2500 (Illumina), collecting 100 base paired-end reads. For plasma cells and plasmablasts, unstranded non-ribosomal RNA–enriched libraries were made using the SMARTer Ultra Low Input RNA Kit for Sequencing v3 (Clontech) and sequenced collecting 50 base paired-end reads.

ChIP immunoprecipitation-sequencing was performed in triplicate for all stages of B-cell development, except plasmablasts and plasma cells, as described previously.²¹  For details, see the supplemental Methods.

RNA-seq reads were aligned to the C57BL/6J mouse reference genome (mm10, GRCm38) using STAR²²  v2.3.0e. Transcript assembly was performed separately for all samples except plasmablasts and plasma cells using Cufflinks²³  v2.2.0. Nonuniquely mapping reads were retained during assembly and potential transcripts (transfrags) were discarded when they contained fewer than 5 successfully mapped reads. Individual assemblies were subsequently compared and transfrags were discarded if they were not assembled in at least 2 samples.

The lncRNA discovery pipeline is depicted in supplemental Figure 2A. Following assembly, transcripts <200 bp in length were discarded. Remaining transcripts were filtered against existing databases (Ensembl v72 and National Center for Biotechnology Information RefSeq) and discarded if they intersected (≥1 bp, on the same strand) intervals annotated as anything other than noncoding RNA. Transcripts were also discarded if they intersected nuclear mitochondrial DNA or pseudogenes predicted using Exonerate,²⁴  if they were in close proximity downstream of a protein-coding gene, or if they were classified as coding by both Coding Potential Calculator²⁵  and PhyloCSF.²⁶  Further details, including a description of lncRNA nomenclature, are provided in supplemental Methods. LncRNAs classed as intergenic (>5 kb from a protein-coding gene) were identified as eRNAs or pRNAs using H3K4me1 and H3K4me3 ChIP-seq data in conjunction with a stringent classification pipeline (supplemental Figure 2B), in which successfully classified loci were required to intersect a called chromatin peak and to show 4-fold greater coverage of reads arising from their characteristic chromatin mark.

Genomic data has been visualized as a UCSC public hub, accessible at https://www.cgat.org/downloads/public/projects/proj010/UCSC_track_hub/hub.txt.

To identify and characterize lncRNAs expressed during B-cell development and activation, we used fluorescence-activated cell sorting to obtain 5 replicate cell samples for each of 8 purified B-cell populations derived from the bone marrow, spleen, and peritoneal cavity of adult female C57BL/6JNimr mice (Figure 1A; supplemental Figure 1A-D). These populations represent the major subsets of developing B cells in the bone marrow (pro-B cells, pre-B cells, and immature B cells) and subsets of mature B cells in the bone marrow, spleen (follicular and marginal zone B cells), and peritoneal cavity (B1a cells). In addition, we also sorted germinal center (GC) B cells, an activated B-cell subset. From each of these 40 samples, we extracted and sequenced the polyA⁺ RNA fraction in a strand-specific manner to an average depth of 42 million uniquely mapped reads, using a 2 × 100 base paired-end sequencing protocol.

To identify lncRNAs, we used a de novo transcript assembly pipeline (supplemental Figure 2A). Stringency of selection of potential lncRNA loci was assured by discarding transcripts whose assembly was not replicated in at least 2 of the 40 sample libraries, by discarding loci represented in public databases of protein-coding annotations (Ensembl and RefSeq), and by discarding loci containing transcripts with predicted coding potential (supplemental Figure 3A-C). Using this process, we identified 4516 lncRNA loci that were expressed in at least 1 of the 8 B-cell populations, and comprised, on average, 1.8 transcripts per locus (supplemental Table 1). To validate the assembly process, we used polymerase chain reaction to check for expression of predicted lncRNAs. Of 53 lncRNAs with expression values ranging between 0.9 and 667 fragments per kilobase per million (FPKM), we were able to successfully validate 47 (89%), demonstrating that the assembly pipeline had generated reliable transcripts (supplemental Figure 4).

Of the predicted lncRNA loci, 3025 (67%) showed no evidence of splicing (single-exon loci), whereas the remaining 1491 showed evidence of at least 1 spliced transcript (multiexon loci). Approximately half of all single and multiexon loci were intergenic, with the remainder being either flanking (<5 kb from a protein-coding gene) or overlapping a protein-coding gene on the antisense strand (Figure 1B; supplemental Table 2). The single-exon and multiexon lncRNAs showed similar coding potential, expression levels and mean exon sizes, supporting the view that they belong to the same gene class (supplemental Figure 5A-D).

We compared the 2349 intergenic lncRNA loci identified in this study with intergenic lncRNAs reported in mouse T-cell subsets,¹⁸  and with the latest Ensembl (v74) long intergenic noncoding RNA annotations. A total of 1829 (78%) of the intergenic lncRNAs in our B-cell catalog did not overlap on the relevant strand with either the T cell or Ensembl long intergenic noncoding RNAs (Figure 1C), reinforcing the notion that many lncRNAs are tissue- or cell type–specific in their expression.^9,12,18 

The highly tissue-specific expression of B-cell lncRNAs was further confirmed by comparison of our catalog with the Functional Annotation of the Mammalian Gene 5 (FANTOM5) study of transcriptional start sites (TSS) in 128 mouse primary cell types,²⁷  which included a single B-cell sample. Because FANTOM annotation required TSS to be detected in 2 or more samples, we did not expect this dataset to include TSS for lncRNAs expressed only in B cells. Accordingly, only 299 (13%) of our intergenic lncRNAs were within 500 bp of a FANTOM-annotated TSS (Figure 1D-E). Notably, a greater overlap with FANTOM-annotated TSS was observed for multiexon transcripts than for single-exon transcripts. The same trend was observed in comparisons with previous lncRNA annotations (supplemental Figure 5E-F). We therefore conclude that expression at an lncRNA locus is more likely to be recapitulated in different cell types if the locus contains evidence of splicing.

A previous study had shown that many lncRNAs overlap with transposable elements (TEs).²⁸  In agreement with this, we find that approximately 40% of the lncRNAs identified here overlap with TEs, and a comparison of these with the genomic content of TEs shows an enrichment of endogenous retroviruses and a relative depletion of LINE elements, similar to the previous report (Figure 1F; supplemental Table 3).

In summary, we have generated a novel, high-confidence catalog of lncRNAs expressed across key stages of B-cell development. To facilitate the use of this catalog, we have visualized our genomic data as a UCSC public hub.

To identify intergenic lncRNAs within our catalog that may be either eRNAs or pRNAs, we determined the genome-wide abundance of H3K4me1 and H3K4me3 chromatin marks known to distinguish these 2 lncRNA subtypes⁸  (Figure 2A-B). We quantified the relative abundance of these marks at the proposed TSS of 2349 intergenic lncRNAs across all 8 B-cell populations, and identified eRNAs and pRNAs using stringent selection criteria (supplemental Figure 2B).

A total of 702 eRNAs and 192 pRNAs were predicted in 1 or more B-cell populations (Figure 2C; supplemental Figure 6A-C; supplemental Table 2) and were evenly distributed between single-exon (320 eRNAs, 96 pRNAs) and multiexon loci (382 eRNAs, 96 pRNAs, Figure 2D). H3K4me1 and H3K4me3 marks at all lncRNA loci remained broadly consistent across B-cell populations (Figure 2E). Furthermore, lncRNA classified as eRNAs or pRNAs in 1 B-cell population displayed a comparable chromatin signature in other cell populations (supplemental Figure 6D-E), suggesting that these classifications represent 2 largely mutually exclusive lncRNA classes.

In agreement with earlier studies,⁸  B-cell eRNAs were distinguishable from pRNAs by a significantly higher correlation of expression with their proximal protein-coding gene (Figure 2F), by lower expression (Figure 2G), and by expression in fewer cell subsets (Figure 2H). Thus, lncRNAs with a high H3K4me1:H3K4me3 ratio at their TSS show characteristics typical of eRNAs and their expression tends to be associated with the expression of adjacent protein-coding genes.

Previous studies in the T-cell lineage had shown that intergenic lncRNA expression is highly restricted even between closely related T-cell subsets.^17,18  If this were also the case for B cells, then B-cell subsets should be distinguishable on the basis of lncRNA expression. Indeed, using both hierarchical clustering and principal component analysis on either coding gene or lncRNA expression, all 8 B-cell populations were distinguishable and replicate samples clustered closely (Figure 3A-B; supplemental Figure 7A-D). For both types of genes, follicular and marginal zone B cells were the most closely related, whereas B1a cells and GC B cells showed the greatest distinction from other cell types.

A greater proportion of B-cell lncRNAs was differentially expressed between cell populations than protein-coding genes (supplemental Figure 8A). Although cell type specificity was, as expected, dependent on expression level, even at equivalent expression levels, lncRNAs displayed greater cell-type restriction than coding genes (Figure 3C; supplemental Figure 8B-C), suggesting that lncRNAs have very restricted spatiotemporal roles during B-cell development.

The trend for lncRNAs to be more restricted in expression than protein-coding genes was observed across all stages of B-cell development. At a defined expression threshold (FPKM >1), a greater proportion of protein-coding genes than lncRNAs were expressed in all bone marrow–derived B-cell populations, in all splenic and peritoneal cavity mature B-cell populations, and in both naïve and activated B-cell populations (Figure 3D-E).

To extend this analysis further, we examined expression of the lncRNAs in antibody-secreting cells. We performed RNA-seq on splenic plasmablasts and plasma cells and on bone marrow plasma cells (supplemental Figure 1E-F; supplemental Figure 7E) and determined expression of lncRNAs and coding genes in these 3 additional subsets (supplemental Table 2). Principal component analysis showed that the 3 subsets clustered closely together (supplemental Figure 7F-G) and lncRNAs were again more cell type–specific (Figure 3D-E).

One genome-scale approach to generate hypotheses regarding lncRNA function is to exploit a “guilt by association” approach by correlating expression of lncRNAs with coding genes of known biological function.^9,29  We predicted the biological processes to which lncRNAs might contribute by clustering protein-coding loci expressed in at least 1 B-cell type (median FPKM >1) into 1 of 17 clusters, based on correlation of expression across B-cell development³⁰  (supplemental Figure 9A-B; supplemental Table 2). We subsequently correlated individual lncRNAs with the representative gene expression profile (eigengene) for each protein-coding cluster (Figure 4A), and defined an association between an lncRNA and a cluster of protein-coding genes at a Pearson correlation coefficient of |ρ| > 0.8 (Figure 4B). At this threshold, 987 (22%) lncRNAs associated with protein-coding gene clusters, suggesting that these lncRNAs contribute to the same biological processes as the protein-coding genes within the respective cluster. Although some clusters were characteristic of individual cell stages (Figure 4A), Gene Ontology enrichment analysis³¹  identified others that were associated with cellular processes (supplemental Table 4). For example, the lightpink4 cluster contains genes up- or downregulated in marginal zone B cells, whereas the floralwhite cluster is enriched for genes associated with the cell cycle, and its eigengene is highest in pro-B cells, pre-B cells, and GC B cells, the 3 subsets with the highest rates of cell division.

We previously showed an increased correlation of expression between eRNAs and their proximal protein-coding gene compared with pRNAs (Figure 2F); the association of lncRNAs with clusters of protein-coding genes allowed us to test this hypothesis further. Permutation tests demonstrated that eRNAs showed strongest correlation with the cluster containing their proximal protein-coding gene (P = .04), whereas no such trend was observed for pRNAs (P = .51).

Combining information about genomic location, chromatin state, and expression allowed us to identify eRNAs that may function in cis with genes of distinctive biological function. For example, Figure 4C shows a single cluster of protein-coding genes that exhibit large changes of expression within marginal zone B cells. Of the lncRNAs correlated with this cluster, we show 1 (LNCGme02323) that is proximal to cluster member Zc3h12c, encoding a zinc finger protein implicated in the regulation of pro-inflammatory activation in macrophages.³²  Using our catalog, LNCGme02323 is identifiable as an eRNA that is coexpressed with Zc3h12c in marginal zone B cells and that lies ∼67 kb upstream of the coding gene TSS (Figure 4D-E), suggesting that this lncRNA may be transcribed from an enhancer controlling the expression of Zc3h12c. In total, we were able to identify 126 eRNA loci whose expression was most strongly correlated with that of a weighted gene coexpression network (WGCNA) module containing an adjacent protein-coding gene, including 48 with a correlation of |ρ| > 0.8 (supplemental Table 5).

To extend the utility of this catalog of mouse lncRNAs, we sought to identify potential orthologous human lncRNAs by comparing the mouse genes with lncRNAs previously identified in human B cells.^17,33-35  Based on pairwise genomic alignments, we identified 185 mouse lncRNAs that have a human syntenic ortholog (supplemental Table 6). Taking advantage of the fact that eRNAs are a distinct class of lncRNA throughout B-cell development (supplemental Figure 6D), we compared mouse eRNAs with human eRNAs reported in CD19⁺ B cells.³⁴  Pairing eRNAs adjacent to orthologous human and mouse coding genes identified a further 228 mouse eRNAs that have 1 or more human eRNA counterparts (supplemental Table 6). Such conserved location indicates that the human/mouse pairs may be marking enhancers of orthologous coding genes.

The transcription factor PAX5 is expressed throughout B-cell development and is crucial for commitment to the B-cell lineage as well as for maintaining B-cell homeostasis.^19,36  PAX5 mediates these functions through the regulation of protein-coding gene expression, but it may also regulate lncRNA expression. Using previously published PAX5 binding sites,³⁶  we identified 784 and 717 lncRNAs in pro-B cells and mature B cells, respectively, whose loci either overlapped or were <1 kb downstream of a PAX5 binding site that had no comparable association with a protein-coding gene. This represented a statistically significant enrichment compared with a null model of PAX5 binding (pro-B cells: 6.1-fold, P < 10⁻⁴; mature B cells: 10.5-fold, P < 10⁻⁴). Using published RNA-seq data for wild-type or Pax5-deficient pro-B cells or mature B cells,³⁶  we further identified 28 and 87 of these lncRNA loci, respectively, for which there is evidence that PAX5 both binds and regulates their expression (supplemental Figure 10A-D; supplemental Table 7).

Disruption of PAX5 binding during early B-cell development is associated with development of B-progenitor acute lymphoblastic leukemia (B-ALL).³⁷  For this reason, we extended our analysis of PAX5 binding and lncRNA expression using RNA-seq data from a mouse model of B-ALL in which Pax5 expression can be regulated by doxycycline.³⁸  Induction of Pax5 expression within leukemic cells causes the tumor cells to differentiate and the leukemia to regress. Comparison of B-ALL cells with and without doxycycline-induced Pax5 expression revealed 331 differentially expressed lncRNAs (q < 0.05), of which a large majority (91%) were upregulated as a result of Pax5 expression (supplemental Table 7). Furthermore, most (184) of these lncRNAs originated from loci bound by PAX5 in pro-B cells or mature B cells (Figure 5A-B). As an example, we show 2 previously annotated lncRNA loci (LNCGme00432 and LNCGme00344) with a novel locus (LNCGme00345), all of which are eRNAs and are bound by PAX5, differentially expressed in the absence of PAX5, and lie downstream of the B-cell lymphoma 11a gene (Bcl11a), whose expression is also PAX5-dependent in B-ALL cells (Figure 5C). We note that Bcl11a does not have PAX5 binding sites near its TSS, and yet its expression is PAX5-dependent. Taking these observations together, we hypothesize that LNCGme00432, LNCGme00344, and LNCGme00345 may be transcribed from PAX5-induced enhancers that in turn regulate expression of Bcl11a and thus explain how expression of Bcl11a is indirectly PAX5-dependent. Because BCL11A has been proposed to activate expression of the Rag1 and Rag2 genes and hence to promote recombination of V, D, and J immunoglobulin genes,³⁹  it may contribute to PAX5-dependent differentiation of B-ALL cells and subsequent regression of the leukemia. Furthermore, BCL11A is essential for B-cell development^40,41 ; thus, these lncRNAs may also mediate the effect of PAX5 on expression of Bcl11a during normal B-cell development. Interestingly, LNCGme00344 and LNCGme00345 have human eRNA orthologs (BMThy_chr2_0943 and BMThy_chr2_0945; supplemental Table 6), suggesting that this transcriptional regulation may be conserved in human B cells. In total, we report 199 lncRNAs (including 73 eRNAs) that are bound by PAX5, show PAX5-dependent expression, and are neighboring a protein-coding gene that also shows PAX5-dependent expression (supplemental Table 7).

To aid in the identification of lncRNAs active within the immune system, we present a comprehensive catalog of 4516 lncRNAs, including 2349 intergenic loci, expressed across 11 murine B-cell populations. Translating genome-scale lncRNA identification into in vivo functional studies remains a major challenge that necessitates the use of model organisms. Previous work has identified lncRNAs expressed in a subset of these populations in humans^17,33-35 ; however, given the high tissue specificity and rapid evolution of lncRNAs, a murine catalog is therefore a prerequisite to successful lncRNA studies in the mouse. In addition, we defined patterns of expression, chromatin modification, and transcription factor binding that provide insight into the regulation and function of these transcripts and we identified a subset of lncRNAs that have human orthologs and thus are of direct relevance to human research.

The majority (67%) of lncRNA loci present in our catalog contained no evidence of splicing. Although single-exon transcripts had to be replicated in 2 or more samples to be included in this catalog, their loci were less likely to be in published datasets. However, the high rate of polymerase chain reaction validation of these single-exon lncRNAs (85%) and their comparable intersection with clearly defined chromatin peaks confirms the reliability of the transcript annotation, and the demonstrated functions of single-exon lncRNAs such as NEAT1⁴²  and Paupar⁴³  highlight the importance of including these transcripts in the catalog.

LncRNA genomic locations have the potential to provide insight into function. A substantial proportion of multiexon (48%) and single-exon (48%) loci were antisense to or flanking (<5 kb) protein-coding genes. Previous studies have classified as many as 60% of detected lncRNAs as antisense transcripts arising from bidirectional transcription at protein-coding promoters⁴⁴  that could regulate expression of genes encoding transcriptional regulators.⁴⁵  These are included within the lncRNAs classified as flanking in this study. Full-length antisense lncRNAs have also been implicated in a variety of cellular functions,^46,47  most notably FAS-AS1, which regulates alternative splicing in its antisense coding gene Fas.^5,48 

Half of all lncRNAs detected in B cells (52%) were intergenic (>5 kb from a protein-coding gene). The use of H3K4me1 and H3K4me3 chromatin marks to distinguish eRNAs from pRNAs has been previously demonstrated,⁸  and we were able to corroborate many of the features of these transcripts. The consistent classification of these loci across B-cell development strongly suggests that these markers identify 2 distinct and mutually exclusive classes of lncRNA.

Comparison with lncRNAs expressed in human B cells identified 185 mouse lncRNAs that had a human ortholog and a further 228 mouse eRNAs that have human eRNA counterparts. This offers the possibility of using the power of mouse genetics to study the in vivo function of these human lncRNAs by genetic modification of the corresponding mouse genes.

The use of expression correlation as a means of associating lncRNA function with that of protein-coding genes has been successfully demonstrated in previous studies.^17,18  By clustering protein-coding genes on the basis of coexpression, we were able to identify discrete groups whose patterns of expression transcend individual B-cell populations and whose members are associated with, for example, cell-cycle function. Post hoc association of lncRNAs with these clusters implicated a substantial proportion (22%) of the lncRNAs with these specific processes. In addition, this analysis supported the hypothesis that eRNAs may function to regulate expression of adjacent protein-coding genes.⁷ 

Using external data, we identified loci regulated by PAX5, a crucial mediator of B-cell development, whose loss contributes critically to the development of B-ALL.³⁸  After incorporating our classification of lncRNAs based on chromatin state, we identified 73 eRNAs situated proximal to protein-coding genes that also show Pax5-dependent expression in B-ALL cells. The roles of such transcripts remain to be investigated; however, their discovery demonstrates the utility of our catalog when used in conjunction with external datasets to probe the complexities of B-cell development and disease.

In conclusion, RNA-seq followed by de novo transcript assembly has the potential to identify thousands of novel long noncoding transcripts. Although a growing number of these transcripts have been assigned to diverse cellular processes, as yet, the majority have no known function. Given the large numbers of these transcripts, a major challenge when seeking to identify and characterize lncRNA function lies in the prioritization of plausible candidates for in vivo functional studies. We anticipate that this catalog will serve as a valuable resource for such future research—for example, by allowing systematic genetic screens of high-confidence lncRNAs to identify those with key functions in B-cell development and activation.

The authors thank Stephen Nutt for Blimp1-GFP mice.

This work was supported by the Flow Cytometry and Advanced Sequencing Facilities of the Medical Research Council (MRC) National Institute for Medical Research (now the Francis Crick Institute). The MRC Computational Genomics Analysis and Training Centre and V.L.J.T. were supported by the Medical Research Council (MRC program number U117527252). V.L.J.T. was supported by the Francis Crick Institute, which receives its funding from the MRC, Cancer Research UK, and the Wellcome Trust.

Contribution: T.F.B., J.S.J., and J.M. performed research; T.F.B., J.S.J., J.M., and A.H. analyzed data; and T.F.B., J.S.J., A.H., C.P.P., and V.L.J.T. designed research and wrote the paper.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Victor L. J. Tybulewicz, The Francis Crick Institute, The Ridgeway, London NW7 1AA, United Kingdom; e-mail: victor.t@crick.ac.uk.