The precise annotation of repetitive regions is central to the accurate assessment of their activities. Until recently, this has relied upon the use of manually curated consensus sequences (Bao et al., 2015) with BLASTn-based search methods to define regions of interest. In place of these flattened representations, hidden Markov models (HMMs) can now also be used to represent repeat families, better representing the full range and variability of their sequence space (Hubley et al., 2016). This profile-based masking improves both accuracy and sensitivity, and annotates an additional 5.5 and 5.1% of the mouse and human genomes, respectively (Hubley et al., 2016). Using this method, the mouse and human genomes (GRCm38.78 and GRCh38.78, respectively) were masked using RepeatMasker¹ configured with nhmmer (Wheeler and Eddy, 2013) in sensitive mode using the Dfam 2.0 library (v150923). RepeatMasker annotates LTR and internal regions separately, complicating the summation of reads spanning these divides. Tabular outputs were, therefore, parsed to merge adjacent annotations for the same element and to produce gene transfer format (GTF) files compatible with popular read-counting programs. GTF files for both genomes are freely available upon request.

The expression data used in this study have been previously described and are publicly available. Ethical review, experimental and methodological details relating to study design and data acquisition can be found in the original reports. The following accessions were used: E-MTAB-2499; GSE61608; GSE60927; GSE68769; GSE65422; GSE60424; GSE72420 and GSE62241, which are a mixture of single-end and paired-end Illumina RNA-seq reads. Adapter contamination, assessed with FastQC², was removed using Trimmomatic (Bolger et al., 2014), with additional quality trimming (Q20) and subsequent length filtering (both reads of a pair ≥ 35 nts). The resulting read pairs were aligned with HISAT2 (Kim et al., 2015) and primary mappings counted with featureCounts (Subread, Liao et al., 2014) using standard GTFs for annotated genes and the curated RepeatMasker GTFs for repeat regions. For accuracy and to prevent ambiguity, only reads that could be uniquely assigned to a single feature were counted. This may underestimate total expression in certain situations, but ensures confident count allocation to individual features. Features with no assigned reads across all samples within an experiment were discarded. Those remaining were normalized to account for variable sequencing depth between samples using DESeq2 (Love et al., 2014). In comparison to the use of normalization to transcripts-per million (TPM), for example, normalized read counts do not facilitate comparison of individual feature expression levels between experiments, but are nevertheless preferable for the assessment of repetitive element expression. Methods normalizing expression to TPM or reads per kilobase million, RPKM, require the accurate knowledge of transcript lengths, which cannot simply be determined for repetitive elements and are, in fact, often variable between treatments and systems. Normalized read counts were subsequently imported into Qlucore Omics Explorer (Qlucore, Lund, Sweden) for all downstream analysis and visualization. This included all statistical comparisons, calculation of fold-changes in transcript abundance, computation of Z-Scores (the number of standard deviations from the mean of each variable for each data point), and plotting either Z-Scores or log₂ fold-changes in heat-map form.

Induction of endogenous MLVs in LPS-stimulated B cells provided one of the earliest examples of ERE modulation upon immune cell activation (Stoye and Moroni, 1983). We therefore focused on murine B cells to examine the transcriptional response of LTR and non-LTR elements to B cell stimulation. To this end, we analyzed RNA-seq data (E-MTAB-2499) from mature B cells, isolated from the spleens of C57BL/6 (B6) mice and stimulated in vitro with LPS, a-IgM antibodies, or a combination of CD40 ligand (CD40L) and IL-4 (Hartweger et al., 2014).

As expected, analysis of this dataset highlighted a strong modulation of a great number of non-viral gene transcripts, with just over half of responding genes (53.6%) upregulated upon stimulation, relative to unstimulated cells (Figure 1A). In contrast, under the same conditions of stimulation, the majority of LTR element and LINE transcripts (85.8 and 89.3%, respectively) were proportionally downregulated (Figure 1A). This apparent reduction in ERE transcription is likely the effect of the increase in overall gene transcription in response to stimulation. Closer inspection of the top 31 LTR elements that were induced in these B cell upon stimulation, revealed several different groups of LTR elements (Figure 1B). However, notable was the over-representation of xenotropic endogenous MLVs (Supplementary Table 1). These included two closely integrated proviruses on chromosome 1, Xmv41 (LTR/ERV1|MuLV-int∼RLTR4_Mm|1|171481146|171489815) and Xmv43 (LTR/ERV1|MuLV-int∼RLTR4_Mm|1|170941521|170950177), which have been previously shown to be LPS-responsive (Young et al., 2014), as well as a previously unlocalized provirus, Xmv45 on chromosome 5 (LTR/ERV1|MuLV-int∼RLTR4_Mm|5|23700579|23709245) (Supplementary Table 1). Four additional xenotropic proviruses, previously uncharacterized due to their location of the Y chromosome (Frankel et al., 1989), were also significantly upregulated (Figure 1B and Supplementary Table 1). However, these were highly homologous with Xmv41, Xmv43 and Xmv45 (Supplementary Figure 1), making it difficult to discern whether the Y-linked proviruses are genuinely expressed or whether they report expression of Xmv41, Xmv43 or Xmv45.

As an independent confirmation of the observed pattern of LTR element transcriptional activation, we analyzed a second set of RNA-seq data (GSE61608) from mature B cells, isolated from spleens of B6 mice and stimulated in vitro with LPS or a-IgM antibodies (Fowler et al., 2015). Again, 12 out of the top 31 LTR elements identified in the previous set, were also significantly induced at the earlier time-point of 2 h in this set and, notably, these included Xmv45 (Figure 2A). Furthermore, Xmv45 was also significantly induced in a third set of RNA-seq data (GSE60927) by longer in vitro stimulation of B6-derived B cells with LPS or a combination of CD40L, IL-4 and IL-5 (Shi et al., 2015) (Figure 2B). More importantly, Xmv45 seems also to be transcriptionally induced in vivo, as splenic plasma cells assessed directly ex vivo, which represent a state of recent B cell activation, showed elevated Xmv45 transcription and clustered closely with in vitro LPS-stimulated B cells (Figure 2B). Lastly, further support for transcriptional induction of Xmv45 in vivo was provided by analysis of RNA-seq data (GSE68769) from B cells isolated from the lymph nodes of mice responding to Influenza A virus vaccination. Indeed, Xmv45 was one of the few LTR elements that were significantly induced over the course of vaccination, despite the fact that Influenza A-specific B cells should constitute only a small fraction of total lymph node B cells that were analyzed (Figure 2C).

We next explored whether the consistency with which Xmv45 transcriptional induction was detected in multiple datasets was explained by the degree of this induction. Indeed, transcription of Xmv45 in stimulated B cells was much higher than any other induced LTR element (Figure 2D) and dwarfed transcription of other MLV proviruses that are either weakly (Emv2) or strongly (Xmv41 and Xmv43) inducible by LPS stimulation (Figure 2E; please note that Xmv45 expression is plotted on a scale that is 20-times higher than the rest). Together, these data suggest that transcription of a small selection of LTR elements, exemplified by Xmv45, is consistently induced in murine B cells by a multitude of in vitro and in vivo stimuli and validate the capacity of our analysis to detect this induction in multiple datasets.

We next explored whether ERE transcriptional modulation as observed during physiological B cell activation overlapped with modulation that may occur following B cell transformation. For this purpose, we compared RNA-seq data (GSE65422) from non-transformed B cells (resting splenic B cells and germinal center B cells analyzed directly ex vivo; and B cells activated in vitro with a-CD40 and a-IgM antibodies) with B cells resembling diffuse large B cell lymphoma (DLBCL) (Zhang et al., 2015). The latter were obtained from mice that develop spontaneous B cell lymphomas as a result of deregulated expression of BCL6 under the immunoglobulin (Ig) heavy chain Iμ promoter and of deregulated activation of the alternative NF-κB pathway by expression of the NF-κB inducing kinase (NIK) under the ROSA26 promoter (Zhang et al., 2015). Both these genetic alterations were restricted to the germinal center lineage by conditional mutagenesis, using the Cγ1-cre transgene (Zhang et al., 2015).

Comparison of a-CD40 and a-IgM in vitro activated B cells with resting B cells in this dataset (Figure 3A), revealed a picture comparable with that from the previous dataset (Figure 1A), with transcriptional activation favoring gene induction, and to a lesser extent LTR element and LINE transcription. In contrast, LTR elements and LINEs were dominating the transcriptional differences between resting B cell and B cell lymphomas, with minimal overlap between B cell lymphomas and activated B cells (Figure 3B). Of note, whereas Xmv45 was still the most induced provirus upon in vitro B cell activation, B cell lymphomas where characterized by significantly elevated expression of Emv2 (Figure 3C). In fact, expression of Emv2 in lymphomas was ∼50-times higher than of Xmv45 (Figure 3C; please note that Emv2 expression is plotted on a scale that is 10-times higher than that of Xmv45). The elevated ecotropic MLV expression likely reflects restoration of Emv2 infectivity, which has been previously observed in cancer cell lines (Li et al., 1999) and immunodeficient animals (Young et al., 2012). In contrast to Xmv45, Emv2 transcription in immune cells is only weakly inducible by LPS, but strongly inducible by the epigenetic derepression through BrdU treatment (Young et al., 2014), suggesting that the primary cause of its upregulation in B cell lymphomas is loss of epigenetic repression, followed by restoration of infectivity.

Consistent with different mechanistic origins of LTR element modulation during B cell activation and B cell transformation, approximately one-third of LTR elements that were transcriptionally induced in B cell lymphomas were also induced either in germinal center B cells or in in vitro activated B cells (in equal proportions between the two), whereas the majority (two-thirds) were unique to B cell lymphomas (Figure 3D).

Given the evolutionary divergence between EREs in different host species, we next asked whether the specificity with which EREs are modulated in murine B cells in distinct conditions, also characterized ERE modulation in B cells from a different host, namely the human. We started by investigating gene and ERE transcriptional modulation in RNA-seq data (GSE60424), generated from peripheral blood B cells, isolated from healthy individuals and those with infectious, degenerative or autoimmune diseases, including Sepsis, Amyotrophic Lateral Sclerosis (ALS), Type 1 Diabetes (T1D) and MS (Linsley et al., 2014).

As might be expected by its acute and severe nature, Sepsis accounted for the majority of the 2,159 genes that were differentially regulated between the studied conditions, with a smaller, but clearly evident signature in MS patients shortly after the first treatment with IFNβ (Figure 4). A comparable number of LTR elements were also differentially expressed between the conditions (Figure 4). Interestingly, transcription of LTR elements appeared more distinct between the conditions, with the exception of T1D, than overall gene expression, with a particularly strong signature in the IFNβ-treated subset of MS patients (Figure 4). Moreover, B cells from these individuals differentially expressed more than twice the number of LINEs than of genes, with clusters of LINEs clearly distinguishing the different conditions, again with a very strong signature evident in the IFNβ-treated subset of MS patients (Figure 4).

To probe further the IFNβ-responsiveness of EREs, we first examined whether the pattern observed in purified B cells was also present when the entire complement of blood cell types was analyzed. We focused on LTR elements as they include phylogenetically more diverse groups than LINEs. Indeed, a sizeable set of the IFNβ-inducible LTR elements upregulated in B cells from IFNβ-treated MS patients was also detectable and highly induced in the same patient group, when RNA-seq data from whole-blood was analyzed (131 of 779 elements, Figure 5A).

We next examined the overlap between LTR elements that are induced by IFNβ treatment of MS patients and those that might be naturally induced in a setting of elevated levels of endogenously produced type I IFN (IFN I). To this end, we analyzed RNA-seq data (GSE72420) from whole-blood samples obtained from Systemic Lupus Erythematosus (SLE) patients (Hung et al., 2015), an autoimmune disease with a strong IFN I signature (Obermoser and Pascual, 2010). The intersection of LTR elements that were upregulated in purified B cells in response to IFNβ treatment of MS patients and those induced in SLE patients as a group, identified 219 common LTR elements (Figure 5B). Importantly, significantly elevated expression of these LTR elements was present in the vast majority (88%; 66 of 75) of SLE patients with a high Interferon Signature Metric (ISM) score, but not in any of the SLE patients with a low ISM score (0/24) or any healthy individuals (0/18) (Figure 5B). These results indicated that an overlapping set of LTR elements (Supplementary Table 2) were responsive to IFN I both in IFNβ-treated MS patients and in SLE patients with a high ISM score.

To explore whether IFN I was preferentially inducing certain LTR groups, we compared the composition of all LTR elements expressed in MS or SLE patients with that of the IFN I-inducible LTR elements shared between the two conditions (Supplementary Table 2). This comparison uncovered significant enrichment for the ERV1 groups as a whole, with members of the LTR48, HERV4, MER41D and HERVFH19 subgroups being frequently responsive to IFN I (Figure 5C). Thus, transcriptional induction of LTR elements by exogenously administered or endogenously produced IFN I displays a certain degree of specificity.

Lastly, we investigated how human B cell transformation might influence ERE transcriptional behavior. To achieve this, we compared ERE transcription in RNA-seq data (GSE62241) from follicular B cell lymphoma and from non-transformed B cells (Koues et al., 2015). The groups included B cells purified from follicular lymphoma biopsies; centrocytes, the non-cycling fraction of germinal center B cells, isolated ex vivo from tonsillar tissues; and activated B cells, isolated from peripheral blood samples and stimulated in vitro using a combination of IL-4, a-CD40, a-IgM, and a-IgD (Koues et al., 2015). Consistent with their original description (Koues et al., 2015), follicular lymphoma B cells differentially expressed a substantial number of genes, in comparison with non-transformed activated B cells or centrocytes (Figure 6A). Notably, a much larger number of EREs were dysregulated in follicular lymphoma B cells, with twice as many LTR elements and four-times as many LINEs upregulated in follicular lymphoma B cells as genes (Figure 6A). More importantly, investigation of the overlap between LTR elements upregulated in follicular lymphoma B cells with those induced in B cells from diseases other than cancer uncovered distinguishable, non-overlapping patterns, with the majority of induced LTR elements specific to B cell lymphoma (Figure 6B). Together, these results suggested that distinct LTR elements are transcriptionally activated in cancer and in other diseases.