LTR-retrotransposon transcriptome modulation in response to endotoxin-induced stress in PBMCs

Background: Human Endogenous Retroviruses (HERVs) and Mammalian apparent LTR-retrotransposons (MaLRs) represent the 8% of our genome and are distributed among our 46 chromosomes. These LTR-retrotransposons are thought to be essentially silent except in cancer, autoimmunity and placental development. Their Long Terminal Repeats (LTRs) constitute putative promoter or polyA regulatory sequences. In this study, we used a recently described high-density microarray which can be used to study HERV/MaLR transcriptome including 353,994 HERV/MaLR loci and 1559 immunity-related genes. Results: We described, for the first time, the HERV transcriptome in peripheral blood mononuclear cells (PBMCs) using a cellular model mimicking inflammatory response and monocyte anergy observed after septic shock. About 5.6% of the HERV/MaLR repertoire is transcribed in PBMCs. Roughly one-tenth [5.7–13.1%] of LTRs exhibit a putative constitutive promoter or polyA function while one-quarter [19.5–27.6%] may shift from silent to active. Evidence was given that some HERVs/MaLRs and genes may share similar regulation control under lipopolysaccharide (LPS) stimulation conditions. Stimulus-dependent response confirms that HERV expression is tightly regulated in PBMCs. Altogether, these observations make it possible to integrate 62 HERVs/MaLRs and 26 genes in 11 canonical pathways and suggest a link between HERV expression and immune response. The transcriptional modulation of HERVs located close to genes such as OAS2/3 and IFI44/IFI44L or at a great distance from genes was discussed. Conclusion: This microarray-based approach revealed the expression of about 47,466 distinct HERV loci and identified 951 putative promoter LTRs and 744 putative polyA LTRs in PBMCs. HERV/MaLR expression was shown to be tightly modulated under several stimuli including high-dose and low-dose LPS and Interferon-γ (IFN-γ). HERV incorporation at the crossroads of immune response pathways paves the way for further functional studies and analyses of the HERV transcriptome in altered immune responses in vivo such as in sepsis.

HERV transcriptome; PBMCs; LPS; Endotoxin tolerance; Signalling pathways; Sepsis

Marine Mommert, Olivier Tabone and François Mallet contributed equally to this work.

Retrovirus-like sequences represent the 8% of the human genome [[1] ]. They consist of some 200,000 Human Endogenous Retroviruses (HERVs) and 240,000 Mammalian Apparent LTR retrotransposons (MaLRs). HERVs are remnants of ancestral and independent retroviral infections within the germ line. The parental integrated retroviral DNA or provirus is flanked by two Long Terminal Repeats (LTR). The 5'LTR contains the promoter and enhancer signals initiating transcription, while the 3'LTR contains the polyadenylation signal terminating transcription. Between the two LTRs lie at least three genes coding for the structural proteins (gag), the enzymatic proteins (pro-pol), and the envelope glycoprotein (env). MaLR structure is similar except for the absence of an env gene. Since endogenization, proviruses have propagated all along the genome by reinfection and retrotransposition events. Due to the general absence of selection pressure, most of the elements contain substitutions, insertions and deletions. However, a few hundred large open reading frames (ORF) remain [[2] ], including env ORFs, of which Syncytins support essential functions in placental development [[3] ]. In addition, numerous HERV elements consist of solitary LTRs, resulting from the loss of coding genes by recombination between the two flanking LTRs. All these mechanisms lead to complex multicopy groups, (reviewed in [[4] –[6] ]). Each group consists of heterogeneous elements, all defective for replication, and thus engaged in a vertical mode of transmission exclusively. While bioinformatics approaches have identified 103 HERV groups and 1 MaLR group [[1] ], only 40 HERV groups have been characterised in wet-lab studies [[7] –[9] ]. Based on the homology of pol sequences between HERVs and exogenous retroviruses, these well-defined groups can be classified as gamma-, beta-, spuma- and epsilon- retroviruses [[10] ].

As repeated elements and due to their organisation into groups, (H)ERVs may be involved in genomic plasticity during evolution, being preferential recombination sites within or between chromosomes [[11] ]. Under physiological conditions, HERV elements are subject to strong epigenetic controls either in terms of methylation or histone code [[12] ]. Nevertheless, in various diverse situations, HERV elements have been shown to be transcribed. HERV transcription has been observed in organ-specific (e.g. brain for multiple sclerosis) and systemic (e.g. lupus erythematosus) autoimmune diseases, and HERV-driven mechanisms involving molecular mimicry and immune dysregulation have been proposed [[13] –[15] ]. HERV expression has also been researched in cancer with regard to the oncogenic properties of infectious retroviruses and epigenetic changes observed in cancer [[12] , [16] ]. The latter highlighted LTR-driven transactivation of cellular proto-oncogenes or the expression of HERV-encoded Env, NP9 or Rec candidate oncogenes. HERVs also contribute to the physiopathology of their host at multiple levels. In brief, (i) solo or proviral LTR can modulate the expression of adjacent cellular genes, in addition to their autonomous function in controlling retroviral expression [[17] –[20] ], (ii) the expression of HERV proteins with conventional retroviral functions can influence the host's physiological or pathological states, like fusion for Syncytin-1 [[21] ], immunomodulation for Env HERV-H and Syncytin-2 [[22] , [23] ], RNA nuclear export for Rec [[24] ], and even viral-like particle formation derived from HML-2 [[25] ], and finally (iii) non-coding HERV-expressed sequences may also be biologically active, e.g. HERV-H loci involved in the maintenance of pluripotency in human cells [[26] ]. In all cases, deciphering an HERV biological function starts with the analysis of its expression, which is far from being simple due to both complex biological mechanisms and technical challenges. Such biological complexity is illustrated by one extensively described HERV-W element, the ERVWE1 locus, encoding placental fusogenic Syncytin-1 Env. ERVWE1 expression is driven by its own 5'LTR U3 promoter and adjacent MaLR LTR enhancer. Restricted expression outside the trophoblast is controlled at the LTRslevel by CpG methylation and/or repressive histone mark H3K9me3, and at the splicing level, at least in part, by H3K36me3 along the intron–exon boundary (reviewed in [[3] ]). Experimentally, the challenge of the individual identification of transcriptionally active HERV loci was recently addressed using NGS [[27] ] and high-density microarray [[28] –[30] ] technologies. Indeed, first [[28] ], second [[29] ] and third [[30] ] generations of custom HERV-dedicated microarrays aimed to solve the antagonism between the specificity of individual locus recognition and exhaustiveness of the HERVome. Although addressing a limited number of groups, the first two generations of HERV-dedicated microarrays confirmed that reproductive organs and solid tumours are major sites of HERV expression, and highlighted the tissue specificity/tropism of expressed HERV elements [[28] , [29] ].

Most of the literature concerning HERVs focuses on autoimmune diseases, cancers and placental physiopathology, all these contexts being associated with local or systemic modulation of the immune response. Indeed, there is a growing line of evidence that HERVs may directly shape and regulate our immune system [[31] –[33] ]. The first evidence of the expression of HERVs beta-retroviruses in PBMCs was reported in healthy volunteers 20 years ago, using pol-based pan retroviral PCR [[34] ], following Northern-blot-based seminal scrutiny [[35] ]. This observation was extended to the HERV-H and HERV-W gamma-retrovirus groups using similar PCR-based technology 20 years later [[36] ]. It was demonstrated that the level of HERV expression in the PBMC compartment is modulated in solid organ cancers, autoimmune diseases, infectious diseases, as well as immunocompromised states, e.g. HML-2 in prostate cancer [[37] ], HERV-W/MSRV in multiple sclerosis [[38] ], HERV-W in EBV-infected multiple sclerosis patients [[39] ], and HML-2 in HIV infected patients [[40] ]. In line with this, HERV expression in PBMCs is modulated by microbial components, the differentiation state of the cell and cytokines. Bacteria-derived components such as LPS increase HERV-W, HERV-K, HERV-H, and decrease HERV-E group expression in monocyte-derived macrophage (MDM) cell lines [[41] ]. Recently, taking advantage of probes targeting HERVs in commercial microarrays, the immune cell activation by microbial signals in vitro induces global modulation of endogenous retroelements [[42] ]. HERV-W, HERV-K and HERV-H RNA levels are increased during monocyte differentiation [[41] ]. HERV expression is modulated by cytokines, as observed for ERV3 in the U-937 monocytic cell line [[43] ], or also for MSRV released from PBMCs stimulated by TNF-α and IFN-γ and inhibited by IFN-α [[44] ]. Moreover, Env of HERV-W and HERV-H groups are expressed on the surface of B cells and monocytes in patients with active multiple sclerosis [[45] ]. Syncytin-1 has also been observed in other pro-inflammatory states in skin-homing non-recirculating mycosis fungoides T cells [[46] ]. Conversely, immunosuppressive properties of Syncytin-1 [[47] , [48] ] and Env HERV-H [[22] ] have been described in different contexts [[49] ]. These data support the hypothesis that HERVs are expressed in inflammatory and immunosuppressive contexts and have direct or indirect interactions with the immune host response.

Thus, the considerable improvement inthe molecular and cellular biology tools helps demonstrate that many HERVs are not silent in many contexts. The use ofunique locus candidate approaches helps understand their role in pathophysiological phenomena. We recently introduced a third generation of HERV-dedicated microarray, in high-density Affymetrix format, which makes it possible to measure HERVs at the individual locus level. It targets almost complete coverage of HERVs and MaLRs LTR-retrotransposons of the GRCh38 version of the human genome and Dfam 1.1 database of repetitive DNA element sequence alignments. The chip also targets more than 1500 human genes mainly coding for immunity-related proteins [[30] ]. We investigated HERVs and MaLRs expression in PBMCs by modelling monocyte endotoxin tolerance, consisting in low-dose LPS priming of PBMCs and mimicking monocyte anergy observed in sepsis patients [[50] –[52] ]. Symmetrically, PBMCs were stimulated by single high-dose LPS mimicking gram-negative bacterial infection and illustrating an inflammatory context. This article i) provides an overview of the HERV transcriptome in PBMCs associated with functional LTR characterisation; ii) shows the stimulus-dependent co-modulation of HERVs/MaLRs and genes, including tolerisable reversible phenotypes; and iii) demonstrates an integrative view of HERVs/MaLRs and genes in canonical immunity pathways, illustrating the multifaceted nature of interactions between LTR retrotransposons and genes.

In order to present an overview of the HERV/MaLR transcriptome in PBMCs, we used various LPS challenges to mimic healthy, inflammatory and immunocompromised states using the previously described endotoxin tolerance model [[50] ]. The transcriptome was scrutinised using a custom Affymetrix HERV-V3 microarray [[30] ] which can discriminate 174,852 HERV elements, 179,142 MaLR elements and a set of 1559 genes. In addition to the gene probesets designed in U133plus2 and HTA Affymetrix formats, the chip contains a set of highly informative probesets corresponding to accurately annotated HERV loci hereinafter referred to as "HERV_prototypes", and two sets of probesets corresponding to roughly annotated HERVs/MaLRs elements called hereafter "HERV_Dfam" and "MaLR_Dfam" (summarised in Additional file 2: Table S1). In brief, the "HERV prototypes" repertoire was retrieved from selected prototype loci. These are sets of loci which maintain the largest open reading frames for gag pol env genes within a proviral structure flanked by two complete LTR sequences. From these elements, a repeatmasker-based alignment procedure retrieved 29,271 loci divided into 42 groups. The "HERV_Dfam" and "MaLR_Dfam" repertoires were retrieved from Dfam, a database of repetitive elements detected by RepBase consensus and based on Hidden Markov Models (HMM) [[53] ] ("Proto versus Dfam" tab, Additional file 2: Table S1). Redundancy between repertoires was removed from the HERV/MaLR repertoire. A summary of the absolute counts and the relative abundance of transcriptionally active elements of the HERV/MaLR transcriptome is given in Fig. 1, at the probeset level. Notably, although the hybridisation assessment quality was equivalent across repeated elements and gene repertoires [[30] ], a higher proportion of gene probesets were transcriptionally active, i.e. 52% (42,560 probesets). Overall, 5.6% of targeted HERVs/MaLRs (71,063 probesets) were transcriptionally active in PBMCs. More precisely, 9.4% of the well-described "HERV_prototypes" repertoire, and 5.5% of HERV_Dfam and MaLR_Dfam were expressed. Among the 9.4% expressed prototype elements, 6.1, 1.7 and 1.6% belonged to gamma-, beta- and spuma/epsilon-like retrovirus classes, respectively. On moving from classes towards groups, all well-defined HERVs groups had expressed loci in PBMCs. Notably, within gamma-retroviruses, the largest HERV-H group is that with the highest proportion of active probesets (940 expressed probesets, 17% of the whole group). The lesser known PRIMA-41 group is the third most represented among gamma-retroviridae (581 probesets, 5% of the whole group). HML-8 (713 probesets) and HML-1 (178 probesets) are the groups with the highest proportion of expressed probesets (both 12%, of the whole corresponding group) among beta-retroviridae. The HML2 group, which is considered to be highly active, only showed 92 expressed probesets (9% of the whole group). Interestingly, although represented by a much smaller subgroup, centromeric HML2 elements seemed to be more expressed than the non-centromeric HML-2 ones. Finally, in spuma-retroviridae, HERV-L group provides the largest amount of active probesets (1197 probesets, 9% of the whole group). It should be noted that we did not observe significant enrichment or depletion of a specific repertoire, class or HERV group, according to healthy, inflammatory or immunocompromised-like experimental conditions (data not shown).

Beyond the basic observation of transcriptional expression, it would be informative to know whether HERVs/MaLRs exhibit any transcription bias due to their structure. The HERV/MaLR transcriptomic activity observed collectively in stimulated and unstimulated PBMCs is depicted at the LTR/gag/pol/env region levels for each group/repertoire (Additional file 2: Table S1). Figure 2a provides a simplified comparative view of the "HERV_prototype" repertoire addressed by the microarray and its related transcriptome. We observed some differences between chip and transcriptome. Active LTRs represent 77.7% of the HERV transcriptome, whereas they represent 71.5% of the chip. This means that LTRs are more represented in active elements than internal HERV/MaLR regions. More specifically, solo LTRs are more abundant than proviral LTRs (5'LTRs and 3'LTRs). Nevertheless, it seems that proviral LTRs are over-represented in the transcriptome (21.4% versus 17.7%), whereas solo LTRs are not (56.3% versus 53.8%). Notably, internal proviral genes are under-represented in the transcriptome (22.3% versus 28.5%). Surprisingly, if we compare the "HERV_prototypes", "HERV_Dfam" and "MaLR_Dfam" repertoire transcriptional activities, the HERV prototypes (14.3%) had higher numbers of expressed loci than Dfam elements (5.2% for HERVs and 4.4% for MaLRs). Attributable LTRs (solo, 5' or 3') are defined as LTRs bearing U3 and U5 regions. For these LTRs, we are able to attribute promoter or polyA functions. In the prototype repertoire, 6404 (31%) were attributable LTRs. According to cut-offs and fold change criteria (see material and methods), a function could be attributed to these LTRs in each of the 45 samples (Fig. 2b). Promoter (Pr) activity was assigned to 15.2% of LTRs and polyadenylation (pA) signal was observed in 11.8%. Most of the LTRs were silent (70.5%) and a minority (2.5%) were classified as readthrough (RdT).

Finally, we wanted to know whether a single LTR had the potential to change status (Fig. 2c). As expected, 56.6% of LTRs (3076 LTRs) were systematically silent in all samples. Interestingly, 27.6% of LTRs shifted from a silent status to a promoter (902 Silent/Pr LTRs) or polyadenylation (598 Silent/pA LTRs) function. Although poorly represented (2%), the same shift is observed for RdT LTRs. Some LTRs are exclusively RdT (20 LTRs) and some shifted from RdT status to a promoter (37 RdT/Pr LTRs) or polyadenylation (51 RdT/pA LTRs) function. A significant proportion of LTRs (13.1%) were in promoter (402 Pr LTRs) or polyA (311 pA LTRs) status in all samples (candidate loci validation, Additional file 8: Figure S6d). Almost no LTR (0.6%) shifted through at least three different features including both promoter and polyA functions. Consequently this observation confirms that the shift from promoter to polyA function is an extremely rare if significant event. The same analysis was performed using the complete HERV/MaLR dataset (Additional file 3: Figure S2). Although results should be considered with caution due to the imprecise annotation of HERV_Dfam and MaLR_Dfam LTRs (Additional file 2: Table S1), the overall trends were similar to those observed at the "HERV_prototypes" repertoire level.

The genomic environment encompassing the functional and silent LTRs is depicted in Additional file 4: Figure S3. As previously observed [[29] ], the gene density ratio is almost 1.2 times higher for promoter LTRs than for silent LTRs. Meanwhile, the proportion of intragenic LTRs is biased towards the antisense representation, with about two-thirds of LTRs being antisense to the gene in which they are located, regardless of their functional category. About 80% of intragenic LTRs overlap with introns. Although no major bias in the genomic environment of intergenic LTRs (Additional file 4: Figure S3b) could be associated with their function, some trends are observed. A plateau is observed with constitutive promoter LTRs, reflecting a slightly lower occurrence of genes in the sense orientation up to 10 kb upstream of these LTRs. Symmetrically, sense gene occurrence apparently rises faster than for antisense genes in the downstream zone of silent LTRs compared to promoter LTRs.

To gain insight into the modulation of genes and HERV/MaLR expression associated with PBMC stimulation conditions, we performed both a hierarchical clustering analysis and a supervised statistical analysis in pairwise conditions. For the record, PBMCs from 5 healthy volunteers were cultured in triplicate, (i) without any additional stimulation (NS for non-stimulated), (ii) with a high concentration of LPS to mimic inflammation, (iii) and primed with a low-dose LPS and latterly boosted with a high-dose LPS defining the so-called endotoxin tolerance model (ET) that mimics monocyte anergy (Fig. 3a). We checked our model by means of TNF-α pro-inflammatory and IL-10 anti-inflammatory cytokine quantification both in supernatants (protein level) and cellular extracts (mRNA level). The combined profiles of TNF-α as a "tolerisable gene" which exhibits a lower response, and as IL10, a "non tolerisable gene", whose expression was increased or unaltered, validate the endotoxin tolerance model (Additional file 5: Figure S4). Figure 3b shows a heatmap from hierarchical clustering of the 1% most variant probesets (genes and HERVs/MaLRs) across samples. We observed 3 main expression profiles: (i) non tolerisable probesets with low expression for the NS condition and high expression for the ET and LPS conditions (top, Fig. 3b), (ii) tolerisable probesets with low expression for NS and ET, and high expression for the LPS conditions (bottom, Fig. 3b), and (iii) down-modulated probesets with high expression for the NS condition, and low expression for the LPS and ET conditions (middle, Fig. 3b). Hence, as stimulation conditions appeared to be strong drivers of PBMC transcriptome modulation, we performed a differential expression analysis. Differentially expressed genes are hereinafter referred to as DEG and differentially expressed HERVs/MaLRs loci referred to as DEL. Volcano plots are depicted at the probeset level for both genes and HERVs/MaLRs. They illustrated that similar amounts of elements appeared to be over- or under-expressed in an inflammatory context (LPS), whereas more gene and HERV/MaLR probesets were down-modulated under tolerance condition (ET) (Fig. 3c, top row).

Overall, differential expression analysis of LPS versus NS conditions identified 785 up-regulated probesets and 847 down-regulated probesets (corresponding to 243 distinct HERV/MaLRloci), while analysis of ET versus LPS conditions merely identified 38 up-regulated probesets and 677 down-regulated probesets (corresponding to 105 distinct HERV/MaLR loci) (adjusted p-value < 0.05, |log2FC| > 1) (Fig. 3c, middle, and Additional file 6: Table S2). Following LPS stimulation, out of the 111 differentially expressed HERV/MaLR elements, 38 belong to HERV_Dfam, 51 to MaLR_Dfam, and 22 to HERV_prototypes. Although DEL analysis does not indicate any group or chromosomal enrichment, some characteristic points could be observed. ERV9, ERV-E4.1, HERV-FRD, HERV-H, HERV-I, HML8, PABL-A and PRIMA4 groups were exclusively modulated under inflammatory condition (LPS vs NS). HERV-Fb, HERV-L, HERV-T and PRIMA41 groups were both modulated for LPS vs NS and ET vs LPS (see below). Among the DELs, 64% consist of solo LTRs and 36% of complete or partial proviruses. After LPS stimulation, PTGS2, IL1B, IL12B, IL6, IL1A and TNF-α pro-inflammatory cytokine as well as CCL20 and PTX3, were the most up-modulated genes. MERTK, CLEC7A, CD36, TIMP2 and CCL13 were the most down-modulated genes. Moreover, we observed, with a pathway analysis, an enrichment of the Death Receptor Signalling and NF-ҡB Signalling pathways (Fig. 3c, last row). Notably, these results highlighted significant inactivation of "LPS-stimulated MAPK signalling". This may be due to negative feedback resulting from the high production of TNF-α following LPS stimulation.

In a tolerance context, among the 47 differentially-expressed HERV/MaLR elements, 20 belonged to HERV_Dfam, 20 to MaLR_Dfam and 7 to HERV_prototypes. As observed for LPS vs NS, no group is enriched or depleted among the most DELs. Notably, 11 elements out of 47 (23.4%) are located on chromosome 12. The HERV-HS49C23, HERV-F and HML6 elements were exclusively modulated in the ET vs LPS condition. Some groups had differentially expressed loci in LPS vs NS and other differentially expressed loci in ET vs LPS, such as HERV-Fb, HERV-L and HERV-T groups. Among the 47 DELs, 86% consist of solo LTRs and 14% of complete or partial proviruses. In a tolerance context, IFNγ, NEFH, MMP10, SERPINB2 and THBD were the most up-modulated genes. MX1, OAS1, CCL15, OAS3, EIF2AK2 and TNF-α were the most down-modulated genes. Consistently, there were 3 inhibited pathways: "role of pattern recognition receptors in recognition of bacteria and viruses", "NF-κB signalling" and "TLR signalling"(Fig. 3c, last row). Conversely, "LXR/RXR signalling" was highly activated, putatively reflecting LXR-induced inactivation of the NF-ҡB signalling pathway leading to the anti-inflammatory macrophage phenotype in atherosclerosis (Fig. 3c, last row; [[54] ]).

Differential expression analysis allowed us to identify modulated HERVs/MaLRs and genes involved in inflammatory and tolerance contexts. Before considering functional linking of HERVs/MaLRs and genes, we validated the results obtained with microarrays through the use of the RT-qPCR reference method.

To confirm modulation of expression observed with the HERV-V3 microarray, we selected the 44 HERV/MaLR probesets for RT-qPCR validation,. Of these, 31 were the most differentially expressed probesets in Dfam repertoires and 13 the most differentially expressed probesets among the prototype repertoire. HERV/MaLR locus-specific RT-qPCR systems were meticulously designed and validated to secure locus specificity (see Additional file 7: Figure S5), leading to 32 primer pairs out of 44 of selected candidates. We confirmed HERV/MaLR modulation on the samples used for microarray analysis, and then on an independent cohort of 6 healthy volunteers. PCR products could be obtained on 23 out of 32 primer pairs (depicted in Additional file 7: Figure S5). Overall, 87% of the detectable elements had concordant profiles with HERV-V3 microarray data. Twenty loci exhibited similar expression profiles in microarray and RT-qPCR experiments and 3 exhibited conflicting profiles (Additional file 8: Figure S6).

Notably, a majority of HERV/MaLR elements exhibited similar patterns as "tolerisable" or "non tolerisable" genes, i.e. divergently modulated following LPS and ET treatments. Figure 4a illustrates the comparable "tolerisable" behaviour of TNF-α, 121601901-HERV0116uL, and 08114670-MALR1129uL HERV loci. Figure 4b depicts the similar "non tolerisable" behaviour of IL10, 070278702-MALR1045uL, and 043166701-MALR1020uL MaLR loci. These two phenotypes were observed regardless of their distance from genes. Nevertheless, 121601901-HERV0116uL, 070278702-MALR1045uL, and 043166701-MALR1020uL are located within OAS3, ITGB8 and MIR3945HG genes, respectively. Conversely, 08114670-MALR1129uL is at a distance of more than 100 kb from the closest gene. Interestingly, as known for TNF-α [[51] ], the tolerisable HERV phenotype was reversed by IFN-γ (Fig. 4a and b, column c). Hence, as HERVs/MaLRs and genes seemed to share similar control of expression following stimulation, we attempted to integrate HERVs/MaLRs in gene pathways based on their common transcriptional behaviour.

Differential expression analysis identified genes and HERV/MaLR loci for which expression was modulated in inflammatory or tolerance contexts. Some HERVs/MaLRs and genes had strongly correlated expression profiles. We sought to identify which HERV elements may be co-expressed with genes which were integrated in regulatory networks. First, using Ingenuity Pathway Analysis (IPA), we identified 13 activated or repressed canonical pathways in inflammatory or tolerance contexts. In 11 out of the 13 pathways, 26 genes belonged exclusively to one pathway. We then identified 72 probesets corresponding to 62 HERV/MaLR loci strongly correlated with the 26 genes (correlation ≥0.8) (Additional file 9: Figure S7a-c, Additional file 10: Table S3). Among those 62 HERV and MaLR loci, 7, 26 and 29 belonged to "HERV_prototypes", HERV_Dfam and MaLR_Dfam repertoires, respectively. This allowed us to build a global network integrating HERV and MaLR loci within gene pathways (see Additional file 9: Figure S7d). Twenty-six out of the 62 HERV/MaLR elements belonged exclusively to one pathway, and 36 HERV/MaLR loci are associated with 2 to 7 pathways. Eleven HERV or MaLR loci were identified in the vicinity of genes (≤ 40 kb), and they were integrated into 1 to 4 pathways. Conversely, although spread on various chromosomes, a bundle of 8 HERV/MaLR elements contributed to both the LXR/RXR activation pathway activated under ET condition and the NF-κB signalling pathway activated during inflammation, via PTGS2 and FLT1 genes (see Additional file 9: Figure S7e).

A gene centred view of the "Role of pattern recognition receptors in recognition of bacteria and viruses" (PRR) pathway is depicted in Fig. 5a. It includes 32 HERV/MaLR loci and 7 genes exclusive to the pathway. Ten out of the 32 HERV/MaLR loci belong exclusively to this pathway. Most of the genes are co-expressed with several retroviral elements, 18 for IFIH1, 18 for IRF7, 17 for OAS3, 15 for OAS2, 7 for PTX3 and 3 for C5AR1, widespread on distinct chromosomes. Interestingly, 4 tolerisable HERV/MaLR loci are present on chromosome 12 within a 39 kb region which overlaps with OAS2 and OAS3 tolerisable genes. Moreover, 4 HERV/MaLR loci are on chromosome 1 in a 32 kb region overlapping with IFI44L and IFI44 genes. RT-PCR amplification of IFI44L, IFI44 genes and 011052301-HERV0472uL HERV loci showed a marked decrease for the ET condition (data not shown). Notably, most of these HERV/MaLR elements expressed belonged to the 3'UTR of several gene transcripts: 121601802-HERV0492uL and 121601901-HERV0116uL for OAS3–201, 121602201-HS49sLRp for OAS2–203, 121602901-MALR1023uL for OAS2–202, and 011052301-HERV0472uL and 011052401-HERV0462uL for IFI44L-201.

Alternatively, an HERV/MaLR-centred view is proposed in Fig. 5b and highlights 18 loci which belong to at least 4 different networks. Eleven out of 15 HERV/MaLR tolerisable loci were shared by 4 pathways (PRR, RIG-1, Interferon and IRF pathways). These 11 elements were co-localised with or close to conventional genes (≤ 10 kb). Non-tolerisable elements 070278702-MaLR1045uL and 043166701-MaLR1020uL and tolerisable element 081146702-MALR1129uL are at the crossroad of several pathways, namely "LXR/RXR activation", "TREM1 signalling", and "Dendritic cell maturation" Altogether, this interplay between HERVs/MaLRs, genes and pathways suggested complex and heterogeneous transcriptional regulation mechanisms.

We used the HERV-V3 chip to provide a first overview of the HERV transcriptome in PBMCs. Healthy, inflammatory and immunocompromised/monocyte anergy states were simulated using an endotoxin tolerance model [[50] ]. We observed that about 5.6% of LTR retrotransposons were transcriptionally active in PBMCs, reaching 9.4% in the well-annotated "HERV_prototypes" repertoire. Previous works estimated the extent of the HERV transcriptome, compared to the HERV genome, to be between 7 and 30% [[29] , [55] –[57] ]. These differences could be related to (i) disease conditions [[29] , [55] , [57] ], (ii) tissue specificity [[29] , [55] ], and (iii) technology [[29] , [55] –[57] ]. Although various intensity thresholds (background/expressed) may also account for such differences, we confirmed in our study that a significant amount of HERV and MaLR sequences thought to be silent are actually expressed. HERV or MaLR regions might be either embarked by conventional genes, or self-induced. Indeed, an LTR could possess promoter activity or harbour polyadenylation signals used for proviral gene expression or adjacent non-retroviral genomic sequences. The vast majority of LTRs were silent (70.5%), and the minority were transcribed via readthrough mechanism (2.5%), regardless of the condition. The large amount of silent LTRs is in line with the accumulation of inactivating mutations [[58] ] as well as silencing by epigenetic mechanisms [[59] , [60] ]. A relatively balanced state was observed between putative promoter (15.2%) and polyadenylation (11.8%) functions. The same overall trends were previously observed in cancer tissues [[29] ], although the amount of silent LTRs was higher in our study (70.5% versus 47%), probably due to the largely increased number of targeted groups. Altogether, integrating both more precise data derived from the "HERV_prototype" repertoire and results from the complete HERV/MaLR dataset, about one tenth of LTRs exhibit a constitutive promoter or polyA function while roughly one quarter of LTRs may shift between silent and promoter or polyA functions. Only a few thousandths of LTRs can shift between promoter and polyA, as initially described in cancer [[29] ], strengthening what we called "operational determinism", i.e. an LTR is predetermined to act as a promoter or a polyA site. Although all "prototype" groups are active in PBMCs, we observed a higher proportion of gamma-retroviruses, including notably the super spreader HERV-H group and the HERV-W group containing Syncytin 1. This is consistent with group-based PCR approaches in PBMCs [[36] ], and in MDM cell lines [[41] ], as well as with the detection/modulation of expression of MSRV/HERV-W and HERV-H loci in PBMCs of healthy donors and multiple sclerosis patients [[15] , [45] , [61] , [62] ]. Concerning beta-retroviruses, HERV-K/HML elements were expressed in PBMCs as previously observed in healthy subjects [[36] , [41] ], in prostate cancer [[37] ], and Henoch-Schönlein purpura [[63] ] patients.Notably, centromeric HML2 elements seemed to be relatively more represented than the other HML-2 elements, as observed in the blood of HIV-infected patients [[64] ]. Nevertheless, this exhaustive HERV-dedicated microarray allowed us to detect expression of poorly characterised groups such as PRIMA41 and MER52A which merit further investigations.

After looking at the HERV transcriptome landscape in PBMCs, we analysed whether individual HERVs/MaLRs could be finely regulated upon stimulation. Hierarchical clustering aggregated samples according to their stimulatory status and DE analysis showed that many HERVs/MaLRs and genes were modulated between conditions. We observed a similar amount of up-regulated and down-regulated elements in the LPS condition compared to NS, as previously shown using group-based RT-PCR systems [[36] , [41] ]. Such dichotomic activity of HERVs/MaLRs was found in the highly inflammatory context of burn patients [[65] , [66] ]. In mice, LPS stress induced primary lymphoid cell-specific production of MLV-ERV virions [[67] ]. More generally, microbes modulated ERV transcription in mice [[42] ]. In a tolerance context, HERVs/MaLRs and genes tended to be down-regulated. Altogether, these observations suggest various levels of HERV/MaLR control, including a fine-tuned control similar to conventional genes.

The NF-ҡB signalling pathway was both up-regulated following LPS stimulation and down-regulated following tolerance induction (LPSvsNS, Z score: 2.121; ETvsLPS, Z score: − 2.646). This is highlighted by the pro-inflammatory cytokines, IL12B, IL6, IL1A, TNF-α modulation, as well as IL10 anti-inflammatory variation, as previously described in such contexts [[50] , [68] –[70] ]. As observed for TNF-α and IL10 genes, HERV/MaLR elements exhibit a dichotomy of tolerisable versus non tolerisable phenotypes, although some elements were found to be down-regulated in both conditions. As illustrated by 121601901-HERV0116uL and 081146702-MALR1129uL loci, HERV/MaLR tolerisation was surprisingly reversible upon IFN-γ addition, as previously described for TNF-α [[51] , [52] ]. Again, this demonstrates that HERVs/MaLRs and genes share similar regulation control following stimulation inducing inflammation or anergy in PBMCs. These results confirm the existence of tight HERV expression regulation control, as previously suggested by the tropism-related behaviour of HERV elements in solid tissues and in particular in reproductive tissues [[3] ] and cancer [[29] ].

The comparable responses of HERVs/MaLRs and genes to high-dose, low-dose LPS and IFN-γ-dependent reversibility of tolerance may be due to either autonomous HERV expression potentially driving gene expression, or gene expression potentially embedding HERV expression. Pathway and co-expression analyses allowed us to include 62 HERV/MaLR loci and 26 genes which are exclusive to 1 pathway into regulatory networks. Eleven canonical pathways were activated or repressed in inflammatory or tolerance contexts. Among the integrated HERV/MaLR elements, ten loci mapped to the 3'UTR of OAS2, OAS3, IFI44L, IFI44, C5AR1 and C5AR2 genes. These retroviral elements may contribute to the post-transcriptional control of these transcripts, including polyA signalling, nucleocytoplasmic transport, translation efficiency, localisation and stabilisation of mRNA [[71] ]. The complex interplay of HERVs and genes can be illustrated by the IFI44L to IFI44 gene region on chromosome 1. The 2 genes and the 4 HERV/MaLR loci located between them are similarly down-regulated under ET condition. Such similar transcriptional expression of IFI44L and IFI44 was previously described in purified CD14 monocytes of patients with Sjogren's syndrome [[72] , [73] ]. Interestingly, dissociated IFI44L versus IFI44 expression was observed in high versus low IFN-γ producers following Leishmania braziliensis stimulation in PBMCs [[74] ]. In addition, among the 4 HERV/MaLR elements, the 011052702-MALR1044uL locus corresponded to an annotated CTCF binding region which defines the boundary between active and heterochromatic DNA. The 011052702-MALR1044uL locus may therefore regulates IFI44L and IFI44 expressions in some situations.

The co-expression of genes with a large number of HERV/MaLR elements scattered among different chromosomes, e.g. 18 HERV/MaLR loci located on 10 distinct chromosomes linked with IFIH1, 6 HERVs/MaLRs located on 6 distinct chromosomes linked with PTX3, suggested that they are part of gene networks regulated by shared transcription factors, as previously proposed [[33] ]. However, we did not observe any tissue-specific transcriptional factor binding sites (TFBS), but an enrichment of the AP-1 binding site with promoter LTR was observed as compared to silent LTR (p-value: 1.14 10− 4, data not shown). The similar modulation of MALR1126 LTR promoters belonging to LXR/RXR and NF-κB networks and located on chromosome 5 (050017102-MALR1126uL) and 10 (100175702-MALR1126uL) may reflect such shared regulation. The integration of LTRs within several pathways is a first suggestion that HERVs/MaLRs and genes are similarly regulated. The HERVs/MaLRs integrated in PRR, RIG1, IFN and IRF pathways and which are mapped into gene transcripts, are mostly tolerisable and carry H3K36me3, a mark of actively transcribed regions in normal haematopoietic cells. Most HERV/MaLR elements, either tolerisable or non-tolerisable and integrated into LXR/RXR, TREM1 and dendritic cell maturation or LXR/RXR and NF-ҡB pathways, are at a distance of more than 25 kb from the closest gene and carry mainly H3K27me3, and occasionally H3K9me3, mark of repressed regions in normal haematopoietic cells. Notably, the IFI44-associated 011052702-MALR1044uL and SLC30A4-associated 150312301-HERV0498uL loci appeared to be decorated with H3K9me3 and H3K27me3 repressive histone marks, respectively. As a general trend, LTRs were screened for particular histone modification signals by overlap with Encode peaks for H3K27ac, H3K27me3, H3K36me3, H3K4me1, H3K4me3 and H3K9me3 in different immune cells; we found enrichment for the H3K9me3 mark for CD14+ monocytes. The LPS stimulations probably modified the local chromatin configuration leading to global modulation of expression, of the gene and HERV/MaLR [[75] ]. Intriguingly, it appears that CD14-positive monocytes are at least twice more enriched in H3K9me3 histone mark than CD4- and CD8-positive T cells, B cells and natural killer cells (data not shown). Taken together, the variability of histone marks, the tolerisable/non-tolerisable HERV/MaLR phenotypes and vicinity of genes suggest a common control of expression between HERV/MaLR and genes. In addition to a contribution to cis-regulation, HERV/MaLR presence at the crossroads of different regulatory networks and conservation in the human population (data not shown), may suggest a role in trans-regulation linking enhancer and promoter regions [[33] , [76] ]. A better understanding of the causal relationship between HERV/MaLR, genes and regulatory pathways would merit further investigations.

This microarray-based approach revealed the expression of about 47,466 distinct HERV loci and identified 951 putative promoter LTRs and 744 putative polyA LTRs in PBMCs. HERV/MaLR expression was shown to be tightly modulated following several stimuli including high-dose and low-dose LPS as well as IFN-γ. This allowed us to propose an integrative view of HERVs/MaLRs and genes in global functional pathways. Further systematic analyses will be required to gain insight on the modulation of expression of HERV/MaLR loci in different haematopoietic cell types, including monocytes, B and T cells, as well as neutrophils and NK cells. This may help decipher the multiple levels of HERV functions in haematopoietic cells, as locally illustrated by the surface or intracellular envelope on monocytes [[45] ] or glial cells [[77] ], or by non-coding elements involved in the control of cell differentiation [[78] ]. From an in vivo point of view, this approach paves the way for systematic deciphering of modulated retroviral elements associated with autoimmune diseases such as systemic lupus erythematosus, inflammatory diseases such as type 1 diabetes [[79] ] inherited autoimmune and auto-inflammatory disorders such as type 1 interferonopathies (reviewed in [[80] , [81] ]), and virus- or drug-induced immunocompromised states [[82] ], as well as resulting from a compensatory response to hyperinflammation such as in sepsis [[83] ]. Notably, it would be of interest to investigate whether the altered histone methylation recently observed for genes in LPS-induced tolerance and in septic patients [[75] , [84] ] may affect HERV expression and contribute to sepsis.

Citrated pouches or heparinised tubes blood were obtained from EFS (Etablissement Français du Sang) and used immediately. According to EFS standardised procedures for blood donation and to provisions of the articles R.1243–49 and following ones of the French Public Health Code, a written non-opposition to the use of his donation for research purposes was obtained from healthy volunteers. The blood donors' personal data were anonymised before blood transfer to our research lab. We obtained the favourable notice of the Local Ethical Committee (Comité de Protection des Personnes Sud-Est II, Bâtiment Pinel, 59 Boulevard Pinel, 69,500 Bron) and the acceptance of the Ministère de la Recherche (declaration DC-2008-64) for handling and conservation of these samples. Peripheral blood mononuclear cells (PBMCs) were isolated with Unisep tube density gradient centrifugation (Eurobio) and washed with sterile PBS (phosphate buffered saline) (Eurobio). The PBMCs were adjusted to 2 × 106 cell/mL and cultured in X-Vivo 20 Medium (Lonza) at 37 °C and 5% CO2. All the experiments were perfomed in triplicate. Lipopolysaccharide was purchased from Sigma-Aldrich and was a mix of Escherichia coli O111:B4, O55:B5 and O127:B8 (Sigma). In this ex vivo endotoxin tolerance model, the PBMCs were first cultured for 15 h without (control group NS and LPS cells), or with 2 ng/ml LPS (ET cells). After washing steps, the PBMCs were incubated a second time for 6 h without (control group NS), or with 100 ng/ml LPS (LPS and ET cells) (Fig. 3a). In this model, when specified in the text, the effects of recombinant human IFN-γ to reverse tolerance effects were studied. Another incubation phase was performed for 24 h with 100 ng/mL of human IFN-γ1b (Miltenyi Biotec) or vehicle, between the two LPS incubations. At the end of the experiments, the supernatants were retrieved and stored at − 80 °C. Pro-inflammatory cytokine TNF-α and anti-inflammatory cytokine IL10 concentrations in the PBMC culture supernatants were detected using commercially-available ELISA kits from R&D System, in accordance with the supplier's recommendations. The cells were harvested, lysed in RLT buffer supplemented with β mercaptoethanol and stored at − 80 °C until further processing. The total RNA was extracted from PBMCs using RNeasy Mini kit (Qiagen) according to the manufacturer's instructions. For each RNA extraction, the residual genomic DNA was digested using the gDNA Eliminator spin column (Qiagen), and directly on RNeasy spin column using RNase-Free DNase Set (Qiagen). RNA quantity and quality were determined using Nanodrop (Thermo Scientific) Bioanalyser 2100 (Agilent) according to the manufacturer's instructions [[51] , [52] ].

HERV-V3 targets 353,994 loci-elements, represented by 4,410,200 probes. The custom HERV GeneChip can discriminate between distinct HERV elements composed of a set of highly informative probesets (located in U3, R, U5 subdomains of solo, 5' and 3' individual LTRs and gag/pol/ env regions), hereinafter referred to as 'HERV prototypes repertoire', and a set of probesets with lower-quality annotations (located in the first third and last third of the complete LTR, and every 2.5 kb in the region in between LTRs), hereinafter referred to as 'HERV/MaLR_Dfam repertoire'. The custom HERV GeneChip also contains probesets targeting LINE1, lncRNA, viruses, and the gene repertoire. The descriptions of the HERVgDB4 database and of the final contents of the HERV-V3 microarray are provided in Additional file 2: Table S1 [[30] ].

The cDNA synthesis and amplification steps were performed using 16 ng of RNA with the Ovation Pico WTA System V2 kit (Nugen) according to the manufacturer's instructions. Five micrograms of amplified purified DNA were fragmented into 50–200 bp fragments and were 3-labeled using the Encore Biotin Module kit (Nugen) according to the manufacturer's instructions. The HERV-V3 microarrays were hybridised at 50 °C for 18 h in an oven with constant stirring (60 rpm). Washing and staining were carried out according to the protocol provided by the manufacturer, using the GeneChip fluidics station 450 (Affymetrix). The arrays were finally scanned using the GeneChip scanner 3000 7G (Affymetrix) fluorometric scanner. Images (DAT files) were converted to CEL files using GCOS software (Affymetrix) [[30] ]. The experimental data generated have been filed with the National Center for Biotechnology Information (NCBI) and are available on the GEO DataSets site under access number GSE108239.

Microarray analysis pre-processing was detailed in supplementary methods (Additional file 11: Supplementary Methods). We chose a specific threshold to define LTR functions for the HERV_Dfam repertoire. We used the dichotomy of probeset signal targeting. More specifically, to retain sensitivity and robustness with regard tofunction assignation, we voluntarily and arbitrarily selected a relatively low expression level cut-off of 24.5 for positive signal attribution coupled with a significant fold change between U3 and U5 signals. Therefore, an LTR was referred to as 'promoter' (Pr) in cases where the signal of the U5-associated probeset was (i) over the threshold, and (ii) at least 3 times higher than its U3 counterpart, and as 'polyadenylation signal' (pA) if the intensity of its U3-associated probeset was (i) over the threshold and (ii) at least 3 times higher than its U5 counterpart. An LTR was assigned as 'readthrough' (RdT) if both U3 and U5 signals (i) were over the threshold, and (ii) without significant fold change between its. Finally, an LTR was assigned as silent if U3 and U5 were both under the threshold; all other remaining LTRs were classified as undetermined (112 LTRs). To visualise HERV and gene co-expression, hierarchical clustering based on correlation distance with the average method was performed on the 1% most variable probesets. Subsequently, comparisons between i) unstimulated PBMCs (NS) and PBMCs stimulated once with LPS (LPS), and ii) tolerant PBMCs re-stimulated with LPS (ET) and LPS were carried out. For all probesets, for differential expression analysis, moderated t-tests were performed (Limma, v3.22.7 [[85] ], and p-values adjusted for multiple testing using the Benjamini-Hochberg procedure [[86] ]. A probeset was considered to be statistically significantly differentially expressed when the absolute log2 Fold Change (|log2FC|) was over 1 and the adjusted-p-value under 0.05. Graphs were generated using ggplot2 (v2.2.0) or pheatmap (v1.0.8). Finally, the Ingenuity Pathways Analysis tool (IPA, Ingenuity® Systems, https://analysis.ingenuity.com) was used to assess upstream regulators, canonical pathways, disease, and functions. Details of this analysis were presented in supplementary methods (Additional file 11: Supplementary Methods). HERVs with a correlation coefficient of over 0.8 with an identified gene were selected for integrative pathway view analysis.

All authors were involved in the analysis and interpretation of data as well as drafting the manuscript or revising it critically for important intellectual content. FM, GM, FV, JT, AP and KBP made substantial contributions to the conception and design of the study, MM, OT and FM designed the experiments, MM performed the endotoxin tolerance model. MM, GO and EC designed and performed the microarray experiments, OT and MM performed the data analyses. MM, OT, JT and FM performed data interpretations, MM, AG and PF performed the primer pairs design, screening and validation, RT-qPCR experiments and statistical analysis, FM, MM, OT, MN and JNV conceived and performed genomic environment and transcription factors enrichments analysis, MM, OT and FM wrote the paper. All authors read and approved the final manuscript.

Citrated pouches or heparinised tubes blood were obtained from EFS (Etablissement Français du Sang) and used immediately. According to EFS standardised procedures for blood donation and to provisions of the articles R.1243–49 and following ones of the French Public Health Code, a written non-opposition to the use of his donation for research purposes was obtained from healthy volunteers. The blood donors' personal data were anonymised before blood transfer to our research lab. We obtained the favourable notice of the Local Ethical Committee (Comité de Protection des Personnes Sud-Est II, Bâtiment Pinel, 59 Boulevard Pinel, 69,500 Bron) and the acceptance of the Ministère de la Recherche (declaration DC-2008-64) for handling and conservation of these samples.

Not applicable.

The authors declare that they have no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The authors would like to express our thanks to Nadia Gaci and Laurence Ganée for their support in the primer design, and Jérémie Becker and Valérie Cheynet for their kind advice. Many thanks to Sophie Ablott who edited and proofread the English in this paper.

This work was supported by bioMerieux SA. MM and OT were supported by doctoral grants from bioMerieux. In addition, OT was supported by Association Nationale de la Recherche et de la Technologie (ANRT). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The microarray expression data have been filed on the NCBI Gene Expression Omnibus and are accessible via GEO accession number GSE108239.