1. Genetics and Genomics
  2. Immunology and Inflammation

Transposable elements regulate thymus development and function

  1. Jean-David Larouche
  2. Céline M Laumont
  3. Assya Trofimov
  4. Krystel Vincent
  5. Leslie Hesnard
  6. Sylvie Brochu
  7. Caroline Côté
  8. Juliette F Humeau
  9. Éric Bonneil
  10. Joel Lanoix
  11. Chantal Durette
  12. Patrick Gendron
  13. Jean-Philippe Laverdure
  14. Ellen R Richie
  15. Sébastien Lemieux
  16. Pierre Thibault
  17. Claude Perreault  Is a corresponding author
  1. Institute for Research in Immunology and Cancer, Université de Montréal, Canada
  2. Department of Medicine, Université de Montréal, Canada
  3. Deeley Research Centre, BC Cancer, Canada
  4. Department of Medical Genetics, University of British Columbia, Canada
  5. Department of Computer Science and Operations Research, Université de Montréal, Canada
  6. Fred Hutchinson Cancer Center, United States
  7. Department of Physics, University of Washington, United States
  8. Department of Epigenetics and Molecular Carcinogenesis, University of Texas M.D. Anderson Cancer Center, United States
  9. Department of Biochemistry and Molecular Medicine, Université de Montréal, Canada
  10. Department of Chemistry, Université de Montréal, Canada
https://doi.org/10.7554/eLife.91037.3
Share this article
Cite this article
  1. Jean-David Larouche
  2. Céline M Laumont
  3. Assya Trofimov
  4. Krystel Vincent
  5. Leslie Hesnard
  6. Sylvie Brochu
  7. Caroline Côté
  8. Juliette F Humeau
  9. Éric Bonneil
  10. Joel Lanoix
  11. Chantal Durette
  12. Patrick Gendron
  13. Jean-Philippe Laverdure
  14. Ellen R Richie
  15. Sébastien Lemieux
  16. Pierre Thibault
  17. Claude Perreault
(2024)
Transposable elements regulate thymus development and function
eLife 12:RP91037.
https://doi.org/10.7554/eLife.91037.3
  2. Download BibTeX
  3. Download .RIS

Abstract

Transposable elements (TEs) are repetitive sequences representing ~45% of the human and mouse genomes and are highly expressed by medullary thymic epithelial cells (mTECs). In this study, we investigated the role of TEs on T-cell development in the thymus. We performed multiomic analyses of TEs in human and mouse thymic cells to elucidate their role in T-cell development. We report that TE expression in the human thymus is high and shows extensive age- and cell lineage-related variations. TE expression correlates with multiple transcription factors in all cell types of the human thymus. Two cell types express particularly broad TE repertoires: mTECs and plasmacytoid dendritic cells (pDCs). In mTECs, transcriptomic data suggest that TEs interact with transcription factors essential for mTEC development and function (e.g., PAX1 and REL), and immunopeptidomic data showed that TEs generate MHC-I-associated peptides implicated in thymocyte education. Notably, AIRE, FEZF2, and CHD4 regulate small yet non-redundant sets of TEs in murine mTECs. Human thymic pDCs homogenously express large numbers of TEs that likely form dsRNA, which can activate innate immune receptors, potentially explaining why thymic pDCs constitutively secrete IFN ɑ/β. This study highlights the diversity of interactions between TEs and the adaptive immune system. TEs are genetic parasites, and the two thymic cell types most affected by TEs (mTEcs and pDCs) are essential to establishing central T-cell tolerance. Therefore, we propose that orchestrating TE expression in thymic cells is critical to prevent autoimmunity in vertebrates.

eLife assessment

This important study shows, based on analyses of single-cell RNA-seq data sets of thymus cells, that transposable elements (TEs) are broadly expressed in thymic stromal cells, especially in medullary thymic epithelial cells and plasamacytoid dendritic cells. The authors also show that at least some TE-derived peptides are presented by MHC-I molecules in the thymus. The study provides solid findings supporting a role of TEs in thymic T-cell selection and immune self-tolerance.

https://doi.org/10.7554/eLife.91037.3.sa0

Introduction

Self/non-self discrimination is a fundamental requirement of life (Boehm, 2012). In jawed vertebrates, the thymus is the only site where T lymphocytes can be properly educated to distinguish self from non-self (Boehm and Swann, 2014; Suo et al., 2022). This is vividly illustrated by Oncostatin M-transgenic mice, where T-cell production occurs exclusively in the lymph nodes (Terra et al., 2005). These mice harbor normal numbers of T-cell receptors (TCRs) αβ T cells but present severe autoimmunity and cannot fight infections (Blais et al., 2008). Intrathymic generation of a functional T-cell repertoire depends on choreographed interactions between the TCRs of thymocytes and peptides presented by major histocompatibility complex (MHC) molecules on various antigen-presenting cells (APCs) (Zuñiga-Pflucker et al., 1989). Positive selection depends on self-antigens presented by cortical thymic epithelial cells (cTECs) and ensures that TCRs recognize antigens in the context of the host’s MHC molecules (Breed et al., 2018; Dervović and Zúñiga-Pflücker, 2010). The establishment of central tolerance depends on two main classes of APCs located in the thymic medulla: dendritic cells (DCs) and medullary TEC (mTEC) (Lebel et al., 2020; Srinivasan et al., 2021; Cheng and Anderson, 2018). Two other APC types have a more limited contribution to central tolerance: thymic fibroblasts and B cells (Perera et al., 2016; Nitta et al., 2011). High avidity interactions between thymic APCs and autoreactive thymocytes lead to thymocyte deletion (negative selection) or generation of regulatory T cells (Treg) (Malhotra et al., 2016).

The main drivers of central tolerance, mTECs and DCs, display considerable phenotypic and functional heterogeneity. Indeed, recent single-cell RNA-seq (scRNA-seq) studies have identified several subpopulations of mTECs: immature mTEC(I) that stimulate thymocyte migration to the medulla via chemokine secretion (Lkhagvasuren et al., 2013), mTEC(II) that express high levels of MHC and are essential to tolerance induction, fully differentiated corneocyte-like mTEC(III) that foster a pro-inflammatory microenvironment (Laan et al., 2021), and finally mimetic mTECs that express peripheral tissue antigens (Michelson et al., 2022). Three different proteins whose loss of function leads to severe autoimmunity, AIRE, FEZF2, and CHD4, have been shown to drive the expression of non-redundant sets of peripheral tissue antigens in mTECs (Ramsey et al., 2002; Takaba et al., 2015; Tomofuji et al., 2020). DCs, on the other hand, are separated into three main populations. Conventional DC 1 and 2 (cDC1 and cDC2) have an unmatched ability to present both endogenous antigens and exogenous antigens acquired via cross-presentation or cross-dressing (Ginhoux et al., 2022). Plasmacytoid DC (pDC) are less effective APCs than cDCs, their primary role being to produce interferon alpha (IFNɑ) (Ginhoux et al., 2022). Notably, thymic pDCs originate from intrathymic IRF8hi precursors, and, in contrast to extrathymic pDCs, they constitutively secrete high amounts of IFNɑ (Colantonio et al., 2011; Lavaert et al., 2020; Le et al., 2020). This constitutive IFNɑ secretion by thymic pDCs regulates the late stages of thymocyte development by promoting the generation of Tregs and innate CD8 T cells (Xing et al., 2016; Hanabuchi et al., 2010; Martín Gayo et al., 2010; Martinet et al., 2015; Epeldegui et al., 2015).

Transposable elements (TEs) are repetitive sequences representing ~45% of the human and mouse genomes (Treangen and Salzberg, 2011; Deniz et al., 2019). Most TEs can be grouped into three categories: the long and short interspersed nuclear elements (LINE and SINE, respectively) and the long terminal repeats (LTRs). These broad categories are subdivided into over 800 subfamilies based on sequence homology (Bourque et al., 2018). TE expression is typically repressed in host cells to prevent deleterious integrations of TE sequences in protein-coding genes (Argueso et al., 2008). Unexpectedly, TEs were recently found to be expressed at higher levels in human mTECs than in any other MHC-expressing tissues and organs (i.e., excluding the testis) (Larouche et al., 2020; Carter et al., 2022), suggesting a role for TEs in thymopoiesis. Since some TEs are translated and generate MHC I-associated peptides (MAP) (Larouche et al., 2020), they might induce TE-specific central tolerance (Kassiotis, 2023). Additionally, TEs provide binding sites to transcription factors (TFs) and stimulate cytokine secretion via the formation of double-stranded RNA (dsRNA) (Chuong et al., 2016; Bogdan et al., 2020; Adoue et al., 2019; Lefkopoulos et al., 2020; Lima-Junior et al., 2021). Hence, TEs could have pleiotropic effects on thymopoiesis. To evaluate the role of TEs in thymopoiesis, we adopted a multipronged strategy beginning with scRNA-seq of human thymi and culminating in MS analyses of the MAP repertoire of mouse mTECs.

Results

LINE, LTR, and SINE expression shows extensive variations during ontogeny of the human thymus

We first profiled TE expression in various thymic cell populations during development. To do so, we quantified the expression of 809 TE subfamilies (classified according to the RepeatMasker annotations) in the scRNA-seq dataset of human thymi created by Park et al., 2020. Cells were clustered in 19 populations representing the main constituents of the thymic hematolymphoid and stromal compartments (Figure 1a, Figure 1—figure supplement 1). The expression of TE subfamilies was quantified at all developmental stages available, ranging from 7 post-conception weeks (pcw) to 40 years of age (Supplementary file 1a). Unsupervised hierarchical clustering revealed three clusters of TE subfamilies based on their pattern of expression during thymic development (Figure 1b, upper panel): (i) maximal expression at early embryonic stages persisting, albeit at lower levels, throughout ontogeny (cluster 1), (ii) an expression specific to a given timepoint (cluster 2), or (iii) a high expression at early embryonic stages that decreases rapidly at later timepoints (cluster 3). LINE and SINE subfamilies were enriched in cluster 1, whereas LTR subfamilies were significantly enriched in clusters 2 and 3 (Figure 1b, lower panel). Expression of individual LINE and SINE subfamilies was highly shared among different cell types (Figure 1d). In contrast, the LTR subfamilies' expression pattern was shared by fewer cell subsets and adopted a quasi-random distribution (Figure 1d). The pattern of expression assigned to TE subfamilies (Figure 1c, innermost track) was not affected by the proportion of cells of different developmental stages (embryonic or postnatal) (Figure 1c, outermost track, and Figure 1—figure supplement 2). This suggests that our observations do not result from a bias in the composition of the dataset. To gain further insights into the expression of TE subfamilies, we studied two biological processes known to regulate TE expression in other contexts: cell proliferation and expression of KRAB zinc-finger proteins (KZFP) (Brocks et al., 2018; Imbeault et al., 2017). Cell cycling scores negatively correlated with TE expression in various thymic cell subsets, particularly for LINE and SINE subfamilies shared among cell types (Figure 1—figure supplement 3 and Supplementary file 1b), whereas analysis of KZFP expression identified ZNF10 as a probable repressor of L1 subfamilies in Th17 and NK cells (Figure 1—figure supplement 4 and Supplementary file 1c). Thus, we conclude that the expression of the three main classes of TEs shows major divergences as a function of age and thymic cell types.

Figure 1 with 4 supplements see all
Download asset Open asset
Long interspersed nuclear elements (LINE), short interspersed nuclear elements (SINE), and long terminal repeats (LTRs) exhibit distinct expression profiles in human thymic cell populations.

(a) UMAP depicting the cell populations present in human thymi (CD4 SP, CD4 single positive thymocytes; CD8 SP, CD8 single positive thymocytes; cTEC, cortical thymic epithelial cells; DC, dendritic cells; DN, double negative thymocytes; DP, double positive thymocytes; Mono/Macro, monocytes and macrophages; mTEC, medullary thymic epithelial cells; NK, natural killer cells; NKT, natural killer T cells; pro/pre-B, pro-B and pre-B cells; Th17, T helper 17 cells; Treg, regulatory T cells; VSMC, vascular smooth muscle cell). Cells were clustered in 19 populations based on the expression of marker genes from Lefkopoulos et al., 2020. (b) Upper panel: heatmap of transposable element (TE) expression during thymic development, with each column representing the expression of one TE subfamily in one cell type. Unsupervised hierarchical clustering was performed, and the dendrogram was manually cut into three clusters (red dashed line). Lower panel: the class of TE subfamilies and significant enrichments in the three clusters (Fisher’s exact tests; ****p≤0.0001). pcw, post-conception week; m, month; y, year. (c) Circos plot showing the expression pattern of TE subfamilies across thymic cells. From outermost to innermost tracks: (i) proportion of cells in embryonic and postnatal samples, (ii) class of TE subfamilies, and (iii) expression pattern of TE subfamilies identified in (b). TE subfamilies are in the same order for all cell types. (d) Histograms showing the number of cell types sharing the same expression pattern for a given TE subfamily. LINE (n = 171), LTR (n = 577), and SINE (n = 60) were compared to a randomly generated distribution (n = 809) (Kolmogorov–Smirnov tests, ****p≤0.0001).

TEs form interactions with transcription factors regulating thymic development and function

TEs provide binding sites to TFs (Chuong et al., 2016; Kunarso et al., 2010; Sundaram et al., 2014), and T-cell development is driven by the coordinated timing of multiple changes in transcriptional regulators (Hosokawa and Rothenberg, 2021). We, therefore, investigated interactions between TE subfamilies and TFs during the development of the human thymus. Two criteria defined an interaction: (i) a significant and positive correlation between the expression of a TF and a TE subfamily in a given cell population, and (ii) the presence of the TF binding motif in the loci of the TE subfamily (Figure 2a). Additionally, we validated the correlations we obtained using a bootstrap procedure to ascertain their reproducibility (see ‘Materials and methods’ for details). This procedure removed weakly correlated TF-TE pairs (Figure 2b). TF-TE interactions were observed in all thymic cell populations (Figure 2c and d, Figure 2—figure supplement 1, and Supplementary file 1d). Numerous TF-TE interactions were conserved between hematolymphoid and stromal cell subsets (Figure 2e). However, the number of interactions and the complexity of the interaction networks were much higher in mTECs than in other cell populations (Figure 2c and d, Figure 2—figure supplements 1 and 2).

Figure 2 with 3 supplements see all
Download asset Open asset
Transposable elements (TEs) shape complex gene regulatory networks in human thymic cells.

(a) The flowchart depicts the decision tree for each TE promoter or enhancer candidate. (b) Density heatmap representing the correlation coefficient and the empirical p-value determined by bootstrap for transcription factor (TF) and TE pairs in each cell type of the dataset. The color code shows density (i.e., the occurrence of TF-TE pairs at a specific point). (c) Connectivity map of interactions between TEs and TFs in medullary thymic epithelial cells (mTECs). For visualization purposes, only TF-TE pairs with high positive correlations (Spearman correlation coefficient ≥ 0.3 and p-value adjusted for multiple comparisons with the Benjamini–Hochberg procedure ≤ 0.05) and TF binding sites in ≥1% of TE loci are shown. (d) Number of TF-TE interactions for each thymic cell population. (e) Sharing of TF-TE pairs between thymic cell types. (f) Number of promoter (top) or enhancer (bottom) TE candidates per TF in hematopoietic cells of the thymus. (g) The proportion of statistically significant peaks overlapping with TE sequences in ETS1 ChIP-seq data from NK cells. (h) Genomic tracks depicting the colocalization of ETS1 occupancy (i.e., read coverage) and TE sequences (in red) in the upstream region of two genes in ETS1 ChIP-seq data from NK cells. Statistically significant ETS1 peaks are indicated by the black rectangles.

Several TFs instrumental in thymus development and thymopoiesis interacted with TE subfamilies (Figure 2—figure supplement 2 and Supplementary file 1d). These TFs include the NFKB1 and REL subunits of the NF-κB complex and PAX1 in mTECs (Baik et al., 2016; Akiyama et al., 2008; Yamazaki et al., 2020) and JUND in thymocytes (Meixner et al., 2004). In DCs, the most notable TF-TE interactions involved interferon regulatory factors (IRF), which regulate the late stages of T-cell maturation, and TCF4, which is essential for pDC development (Xing et al., 2016; Cisse et al., 2008). This observation is consistent with evidence that TEs have shaped the evolution of IFN signaling networks (Chuong et al., 2016). Finally, we found significant interactions between CTCF and TE subfamilies in mTECs and endothelial cells, suggesting that the binding of CTCF to TE sequences affects the tridimensional structure of the chromatin in the thymic stroma (Choudhary et al., 2020). Interestingly, LINE and SINE subfamilies that occupy more genomic space interacted with higher numbers of transcription factors (Figure 2—figure supplement 3).

Using data from the ENCODE consortium for hematopoietic cells (Consortium, 2012; Luo et al., 2020), we looked at the histone marks at the TE loci identified as TF interactors by our analyses (i.e., correlated with TF expression and containing the TF binding motif). The objective was to determine if they could act as promoters or enhancers (Figure 2a and Supplementary file 1e). We found several TE promoter and enhancer candidates in all eight hematopoietic cell types analyzed, with a striking overrepresentation of LINE and SINE compared to LTR sequences (Figure 2f and Supplementary file 1f). Finally, we analyzed publicly available ChIP-seq data of ETS1, an important TF for NK cell development (Taveirne et al., 2020), to confirm its ability to bind TE sequences. Indeed, 19% of ETS1 peaks overlap with TE sequences (Figure 2g). Notably, ETS1 peaks overlapped with TE sequences (Figure 2h, in red) in the promoter regions of PRF1 and KLRD1, two genes critical for NK cells’ effector functions (Kim et al., 2014; Gunturi et al., 2004). Hence, our data suggest that TEs affect thymic development and function by providing binding sites to multiple TFs.

TEs are highly and differentially expressed in human thymic APC subsets

We next sought to determine whether the high expression of TEs reported in mTECs (Bourque et al., 2018; Argueso et al., 2008) was limited to this cell subset or was found in other thymic cell types. Since several thymic stromal cells reach maturity after birth (Bornstein et al., 2018), we selected postnatal samples for the following analyses. We computed two distinct Shannon entropy indices: one for the global diversity of TEs expressed by all cells of a given population and another for the median value of TE diversity expressed by individual cells of a population (Figure 3a). Then, we computed a linear model to represent the diversity of TEs expressed by a cell population based on the diversity of TEs expressed by individual cells (Figure 3a, blue curve). Two salient findings emerged from this analysis. First, the diversity of TEs expressed in the T-cell lineage decreases during differentiation according to the following hierarchy: DN thymocytes > DP thymocytes > SP thymocytes (Figure 3a, Figure 3—figure supplement 1). Second, among the populations of thymic APCs implicated in positive and negative selection (Figure 3a, orange dots), cTECs, mTECs, and DCs expressed broader repertoires of TEs than B cells and fibroblasts. While cTECs and DCs expressed highly diverse TE repertoires at both the population and individual cell levels, the breadth of TE expression in mTECs was found only at the population level (Figure 3a). Accordingly, intercellular heterogeneity (i.e., deviation from the linear model) was higher for mTECs than other cell populations (Figure 3b).

Figure 3 with 3 supplements see all
Download asset Open asset
Human plasmacytoid dendritic cells (pDCs) and mTEC(II) express diverse and distinct repertoires of transposable element (TE) sequences.

(a) Diversity of TEs expressed by thymic populations measured by Shannon entropy. The x and y axes represent the median diversity of TEs expressed by individual cells in a population and the global diversity of TEs expressed by an entire population, respectively. The equation and blue curve represent a linear model summarizing the data. Thymic antigen-presenting cell (APC) subsets are indicated in orange. (b) Difference between the observed diversity of TEs expressed by cell populations and the one expected by the linear model in (A). (c) UMAP showing the subsets of thymic APCs (aDC, activated DC; cDC1, conventional DC1; cDC2, conventional DC2). (d) Bar plot showing the number and class of differentially expressed TE subfamilies between APC subsets. (e) Frequency of expression of TE subfamilies by the different APC subsets. The distributions for pDCs and mTEC(II) are highlighted in bold. mTEC, medullary thymic epithelial cell.

We next focused on thymic APCs expressing the broadest TE repertoires: cTECs, mTECs, and DCs (Figure 3a). To this end, we annotated these APC subpopulations based on previously published lists of marker genes (Figure 3c, Figure 3—figure supplement 2; Park et al., 2020; Bautista et al., 2021). We performed differential expression analyses to determine whether some TE subfamilies were overexpressed in specific APC subsets. pDCs and mTEC(II) overexpressed a broader TE repertoire than other APCs: 32.01% of subfamilies were overexpressed in pDCs and 10.88% in mTEC(II) (Figure 3d and Supplementary file 1g). The nature of the overexpressed TEs differed between pDCs and other thymic APC subsets. Indeed, pDCs overexpressed LTRs, LINEs, and SINEs, including several Alu and L1 subfamilies (Figure 3d and Supplementary file 1g). In contrast, other thymic APCs predominantly overexpressed LTRs.

TE expression showed wildly divergent levels of intercellular heterogeneity in APC subsets. Indeed, whereas most TE subfamilies were expressed by <25% of cells of the mTEC(II) population, an important proportion of TEs were expressed by >75% of pDCs (Figure 3e). To evaluate this question further, we compared TE expression between metacells of thymic APCs; metacells are small clusters of cells with highly similar transcription profiles. This analysis revealed that overexpression of TE subfamilies was shared between pDC metacells but not mTEC(II) metacells, reinforcing the idea that TE expression adopts a mosaic pattern in the mTEC(II) population (Figure 3—figure supplement 3). We conclude that cTECs, mTECs, and DCs express broad TE repertoires. However, two subpopulations of thymic APCs clearly stand out. pDCs express an extremely diversified repertoire of LTRs, SINEs, and LINEs, showing limited intercellular heterogeneity, whereas the mTEC(II) population shows a highly heterogeneous overexpression of LTR subfamilies.

TE expression in human pDCs is associated with dsRNA structures

The high expression of a broad repertoire of TE sequences in thymic pDCs was unexpected (Figure 3d). LINE and SINE subfamilies, in particular, were highly and homogeneously expressed by thymic pDCs (Figure 4a). Constitutive IFNɑ secretion is a feature of thymic pDCs not found in extrathymic pDCs. We, therefore, hypothesized that this constitutive IFNɑ secretion by thymic pDCs might be mechanistically linked to their TE expression profile. We first assessed whether thymic and extrathymic pDCs have similar TE expression profiles by reanalyzing scRNA-seq data from human spleens published by Madissoon et al., 2019; Figure 4—figure supplement 1a and b. This revealed that extrathymic pDCs express TE sequences at similar or lower levels than other splenic cells (Figure 4—figure supplement 1c and d). We then used pseudobulk RNA-seq methods to perform a differential expression analysis of TE subfamilies between thymic and splenic pDCs. This analysis confirmed that TE expression was globally higher in thymic than in extrathymic pDCs (Figure 4b). Since TE overexpression can lead to the formation of dsRNA (Lefkopoulos et al., 2020; Lima-Junior et al., 2021), we investigated if such structures were found in thymic pDCs. pDCs were magnetically enriched from primary human thymi following labeling with anti-CD303 antibody (a marker of pDCs). Then, pDC-enriched thymic cells were stained with an antibody against CD123 (another marker of pDCs) and the J2 antibody that stains dsRNA. The intensity of the J2 signal was more than tenfold higher in CD123+ relative to CD123- cells (Figure 4c and d). We conclude that thymic pDCs contain large amounts of dsRNAs. To evaluate if these dsRNAs arise from TE sequences, we analyzed in thymic APC subsets the proportion of the transcriptome assigned to two groups of genomic sequences known as important sources of dsRNAs: TEs and mitochondrial genes (Sadeq et al., 2021). Strikingly, whereas the percentage of reads from mitochondrial genes was typically lower in pDCs than in other thymic APCs, the proportion of the transcriptome originating from TEs was higher in pDCs (~22%) by several orders of magnitude (Figure 4—figure supplement 2). Finally, we performed gene set enrichment analyses to ascertain if the high expression of TEs by thymic pDCs was associated with specific gene signatures. These analyses highlighted signatures of antigen presentation, immune response, and interferon signaling in thymic pDCs (Figure 4e and Supplementary file 1h). Notably, thymic pDCs harbored moderate yet significant enrichment of gene signatures of RIG-I and MDA5-mediated IFN ɑ/β signaling compared to all other thymic APCs (Figure 4e and Supplementary file 1h). Altogether, these data support a model in which the high and ubiquitous expression of TEs in thymic pDCs would lead to the formation of dsRNAs triggering innate immune sensors, which might explain their constitutive secretion of IFN ɑ/β.

Figure 4 with 2 supplements see all
Download asset Open asset
Transposable element (TE) expression in human plasmacytoid dendritic cells (pDCs) is associated with dsRNA formation and type I IFN signaling.

(a) Frequency of long interspersed nuclear elements (LINE), long terminal repeat (LTR), and short interspersed nuclear elements (SINE) subfamilies expression in thymic pDCs. (b) Differential expression of TE subfamilies between splenic and thymic pDCs. TE subfamilies significantly upregulated or downregulated by thymic pDCs are indicated in red and blue, respectively (Upregulated, log2(Thymus/Spleen) ≥ 1 and adj. p≤0.05; Downregulated, log2(Thymus/Spleen) ≤ –1 and adj. p≤0.05). (c, d) Immunostaining of dsRNAs in human thymic pDCs (CD123+) using the J2 antibody (n = 3). (c) One representative experiment. Three examples of CD123 and J2 colocalization are shown with white arrows. (d) J2 staining intensity in CD123+ and CD123- cells from three human thymi (Wilcoxon rank-sum test, ****p-value≤0.0001). (e) UpSet plot showing gene sets enriched in pDCs compared to the other populations of thymic antigen-presenting cells (APCs). On the lower panel, black dots represent cell populations for which gene signatures are significantly depleted compared to pDCs. All comparisons where gene signatures were significantly enriched in pDCs are shown.

AIRE, CHD4, and FEZF2 regulate distinct sets of TE sequences in murine mTECs

The essential role of mTECs in central tolerance hinges on their ability to ectopically express tissue-restricted genes, whose expression is otherwise limited to specific epithelial lineage (Sansom et al., 2014; Pierre et al., 2017). This promiscuous gene expression is driven by AIRE, CHD4, and FEZF2 (Ramsey et al., 2002; Takaba et al., 2015; Tomofuji et al., 2020). We, therefore, investigated the contribution of these three genes to the expression of TE subfamilies in the mTEC(II) population (Figure 3d). First, we validated that mTEC(II) express AIRE, CHD4, and FEZF2 in the human scRNA-seq dataset (Figure 5a). Next, we analyzed published murine mTEC RNA-seq data to assess the regulation of TE sequences by AIRE, CHD4, and FEZF2. Differential expression analyses between knock-out (KO) and wild-type (WT) mice showed that these three factors regulate TE sequences, but the magnitude and directionality of this regulation differed (Figure 5b and Supplementary file 1i). Indeed, while CHD4 had the biggest impact on TE expression by inducing 433 TE loci and repressing 463, FEZF2’s impact was minimal, with 97 TE loci induced and 60 repressed (Figure 5b). Besides, AIRE mainly acted as a repressor of TE sequences, with 326 loci repressed and 171 induced (Figure 5b). Interestingly, there was minimal overlap between the TE sequences regulated by AIRE, CHD4, and FEZF2, indicating that they have non-redundant roles in TE regulation (Figure 5c). Additionally, AIRE, CHD4, and FEZF2 preferentially targeted LTR and LINE elements, with significant enrichment of specific subfamilies such as MTA_Mm-int and RLTR4_Mm that are induced by Aire and Fezf2, respectively (Figure 5d and Figure 5—figure supplement 1a). While AIRE and CHD4 preferentially targeted evolutionary young TE sequences, the age of the TE sequence did not seem to affect the regulation by FEZF2 (Figure 5—figure supplement 1b). We also noticed that the distance between regulated TE loci was smaller than the distributions of randomly selected TEs (Figure 5—figure supplement 1c). This suggests that AIRE, CHD4, and FEZF2 nonrandomly affect the expression of TE sequences located in specific genomic regions. We observed no significant differences in the genomic localization of TE loci targeted by AIRE, CHD4, and FEZF2 relative to the genomic localization of all TE sequences in the murine genome: most TE loci were located in intronic and intergenic regions (Figure 5—figure supplement 1d). Enrichment for intronic TEs could not be ascribed to induction of global intron retention: the intron retention ratio was similar for TEs regulated or not by AIRE, CHD4, and FEZF2 (Figure 5—figure supplement 1e). ChIP-seq-based analysis of permissive histone marks showed that TE loci induced by AIRE, CHD4, and FEZF2 were all marked by H3K4me3 (Figure 5e). As a proof of concept, we validated that 31.42% of AIRE peaks overlap with TE sequences by reanalyzing ChIP-seq data, confirming AIRE’s potential to bind TE sequences (Figure 5f). Hence, AIRE, CHD4, and FEZF2 regulate the expression of small yet non-redundant repertoires of TE sequences associated with permissive histone marks.

Figure 5 with 1 supplement see all
Download asset Open asset
AIRE, FEZF2, and CHD4 regulate non-redundant sets of transposable elements (TEs) in murine medullary thymic epithelial cells (mTECs).

(a) Expression of AIRE, CHD4, and FEZF2 in human TEC subsets. (b) Differential expression of TE loci between wild-type (WT) and Aire-, Chd4-, or Fezf2-knockout (KO) mice (Induced, log2(WT/KO) ≥ 2 and adj. p≤0.05; Repressed, log2(WT/KO) ≤ –2 and adj. p≤0.05). p-Values were corrected for multiple comparisons with the Benjamini–Hochberg procedure. The numbers of induced (red) and repressed (blue) TE loci are indicated on the volcano plots. (c) Overlap of TE loci repressed or induced by AIRE, FEZF2, and CHD4. (d) Proportion of TE classes and subfamilies in the TE loci regulated by AIRE, FEZF2, or CHD4, as well as all TE loci in the murine genome for comparison (chi-squared tests with Bonferroni correction, **adj. p≤0.01, ***adj. p≤0.001). (e) Plots for the tag density of H3K4me3 and H3K4me2 on the sequence and flanking regions (3000 base pairs) of TE loci induced by AIRE, FEZF2, and CHD4. (f) Proportion of statistically significant peaks overlapping TE sequences in AIRE ChIP-seq data from murine mTECs.

TEs are translated and presented by MHC class I molecules in murine TECs

Several TEs are translated and generate MAPs (Larouche et al., 2020). Hence, the expression of TEs in cTECs and even more in mTECs raises a fundamental question: do these TEs generate MAPs that would shape the T cell repertoire? Mass spectrometry (MS) is the only method that can faithfully identify MAPs (Shapiro and Bassani-Sternberg, 2023; Vizcaíno et al., 2020; Kubiniok et al., 2022). Despite its quintessential role in central tolerance, the MAP repertoire of mTECs has never been studied by MS because of the impossibility of obtaining sufficient mTECs for MS analyses: mTECs represent ≤1% of thymic cells, and they do not proliferate in vitro. To get enough cTECs and mTECs for MS analyses, we used transgenic mice that express cyclin D1 under the control of the keratin 5 promoter (K5D1 mice). These mice develop dramatic thymic hyperplasia, but their thymus is morphologically and functionally normal (Robles et al., 1996; Klug et al., 2000; Ohigashi et al., 2019). Primary cTECs and mTECs (two replicates of 70 × 106 cells from 121 and 90 mice, respectively) were isolated from the thymi of K5D1 mice as described (Dumont-Lagacé et al., 2019). Following cell lysis and MHC I immunoprecipitation, MAPs were analyzed by liquid chromatography MS/MS (Figure 6a). To identify TE-coded MAPs, we generated a TE proteome by in silico translation of TE transcripts expressed by mTECs or cTECs, and this TE proteome was concatenated with the canonical proteome. MS analyses enabled the identification of a total of 1636 and 1714 MAPs in mTECs and cTECs, respectively. From these, we identified four TE-derived MAPs in mTECs and two in cTECs, demonstrating that TEs can be translated and presented by MHC I in the thymic cortex and medulla (Figure 6b and Supplementary file 1j). These MAPs were coded by the three major groups of TE: LINEs (n = 1), LTRs (n = 1), and SINEs (n = 4). Next, we evaluated whether the low number of TE MAPs identified could result from mass spectrometry detection limits (Ghosh et al., 2020; Nanaware et al., 2021). We measured the level and frequency of TE expression in two subsets of cTECs (Figure 6c, left) or mTECs (Figure 6c, right) using scRNA-seq data from Baran-Gale et al., 2020. TE subfamilies generating MAPs in cTECs or mTECs are highlighted in red in their respective plots. Strikingly, TECs highly and ubiquitously expressed the MAP-generating TE subfamilies. These results suggest that the contribution of TEs to the MAP repertoire of cTECs and mTECs might be significantly underestimated by the limits of detection of MS. This is particularly true for mTECs because they express high levels of TEs (Figure 3d), but their TE profile displays considerable intercellular heterogeneity (Figure 3e, Figure 3—figure supplement 2). Nonetheless, our data provide direct evidence that TEs can generate MAPs presented by cTECs and mTECs, which can contribute to thymocyte education.

Figure 6
Download asset Open asset
Murine cortical thymic epithelial cells (cTEC) and medullary thymic epithelial cell (mTEC) present transposable element (TE) MAPs.

(a) mTECs and cTECs were isolated from the thymi of K5D1 mice (n = 2). The peptide-MHC I complexes were immunoprecipitated independently for both populations, and MAPs were sequenced by MS analyses. (b) Number of long interspersed nuclear elements (LINE)-, long terminal repeat (LTR)-, and short interspersed nuclear elements (SINE)-derived MAPs in mTECs and cTECs from K5D1 mice. (c) Distributions of TE subfamilies in murine TECs subsets based on expression level (x-axis) and frequency of expression (y-axis).

Discussion

TEs are germline-integrated parasitic DNA elements that comprise about half of mammalian genomes. Over evolutionary timescales, TE sequences have been co-opted for host regulatory functions. Mechanistically, TEs encode proteins and noncoding RNAs that regulate gene expression at multiple levels (Bourque et al., 2018; Frank and Feschotte, 2017). Regulation of IFN signaling and triggering innate sensors are the best-characterized roles of TEs in the mammalian immune system (Kassiotis, 2023). TEs are immunogenic and can elicit adaptive immune responses implicated in autoimmune diseases (Larouche et al., 2020; Kassiotis, 2023; Gröger and Cynis, 2018; Volkman and Stetson, 2014). Pervasive TE expression in various somatic organs means that co-evolution with their host must depend on establishing immune tolerance, a concept supported by the highly diversified TE repertoire expressed in mTECs (Larouche et al., 2020). This observation provided the impetus to perform multiomic studies of TE expression in the thymus. At the whole organ level, we found that TE expression showed extensive age- and cell lineage-related variations and was negatively correlated with cell proliferation and expression of KZFPs. The negative correlation between TE expression and cell cycle scores in the thymus is coherent with recent data showing that transcriptional activity of L1s is increased in senescent cells (De Cecco et al., 2019). A potential rationale for this could be to prevent deleterious transposition events during DNA replication and cell division. On the other hand, the contribution of KZFPs to TE regulation in the thymus is likely underestimated due to their typically low expression (Huntley et al., 2006) and scRNA-seq detection limit. Additionally, TEs interact with multiple TFs in all thymic cell subsets. This is particularly true for the LINE and SINE subfamilies that occupy larger genomic spaces. Notably, TEs appear to play particularly important roles in two cell types located in the thymic medulla: mTECs and pDCs.

As mTECs are the APC population crucial to central tolerance induction, their high and diverse TE expression is poised to impact the T cell repertoire’s formation profoundly. The extent and complexity of TF-TE interactions were higher in mTECs than in all other thymic cell subsets. These interactions included PAX1 and subunits of the NF-κB complex (e.g., RELB). PAX1 is essential for the development of TEC progenitors (Yamazaki et al., 2020), and RELB is for the development and differentiation of mTECs (Mouri et al., 2014). RelB-deficient mice have reduced thymic cellularity, markedly fewer mTECs, lack Aire expression, and suffer from autoimmunity (Akiyama et al., 2008; O’Sullivan et al., 2018). Under the influence of Aire, Fezf2, and Chd4, mTECs collectively express almost the entire exome (Sansom et al., 2014; Pierre et al., 2017). However, the expression of all genes in each mTEC would cause proteotoxic stress (Pierre et al., 2017). Hence, promiscuous expression of tissue-restricted genes in mTECs adopts a mosaic pattern: individual tissue-restricted genes are expressed in a small fraction of mTECs (Michelson et al., 2022; Klein et al., 2014). The present work shows that mTECs also express an extensive repertoire of TEs in a mosaic pattern (i.e., with considerable intercellular heterogeneity). Aire, Fezf2, and Chd4 regulate non-redundant sets of TEs and preferentially induce TE sequences associated with permissive histone marks. The immunopeptidome of thymic stromal cells is responsible for thymocyte education and represents one of the most fundamental ‘known unknowns’ in immunology. Inferences on the immunopeptidome of thymic stromal cells are based on transcriptomic data. However, (i) TCRs interact with MAPs, not transcripts, and (ii) the MAP repertoire cannot be inferred from the transcriptome (Shapiro and Bassani-Sternberg, 2023; Caron et al., 2011; Admon, 2023). Using K5D1 mice presenting prominent thymic hyperplasia, we conducted MS searches of TE MAPs, identifying four TE MAPs in mTECs and two in cTECs. These results demonstrate that cTECs and mTECs present TE MAPs and suggest they present different TE MAPs. However, the correlation between transcriptomic and immunopeptidomic data suggests that TECs can present many more TE MAPs. Their profiling will require MS analyses of enormous numbers of TECs or the development of more sensitive MS techniques. As TE MAPs have been detected in normal and neoplastic extrathymic cells (Larouche et al., 2020; Laumont et al., 2018; Burbage et al., 2023; Shah et al., 2023), the presentation of TEs by mTECs is likely essential to central tolerance. In line with vibrant plaidoyers for a collaborative Human Immunopeptidome Project (Vizcaíno et al., 2020; Shao et al., 2018), our work suggests that immunopeptidomic studies should not be limited to protein-coding genes (2% of the genome) but also encompass non-coding sequences such as TEs.

The second population of cells exhibiting high TE expression, pDCs, are mainly seen as producers of IFN ɑ/β and potentially as APCs (Ginhoux et al., 2022). Thymic and extrathymic pDCs are ontogenically and functionally different. They develop independently from each other from different precursor cells (Lavaert et al., 2020; Le et al., 2020; Weijer et al., 2002). IFN ɑ/β secretion is inducible in extrathymic pDCs but constitutive in thymic pDCs (Ginhoux et al., 2022; Colantonio et al., 2011). In line with the location of pDCs in the thymic medulla, their constitutive IFN ɑ/β secretion is instrumental in the terminal differentiation of thymocytes and the generation of Tregs and innate CD8 T cells (Xing et al., 2016; Hanabuchi et al., 2010; Martín Gayo et al., 2010; Martinet et al., 2015; Epeldegui et al., 2015). We report here that high TE expression is also a feature of thymic, but not extrathymic, pDCs. Thus, the present study provides a rationale for the constitutive IFN ɑ/β secretion by thymic pDCs: they homogeneously express large numbers of TEs (in particular LINEs and SINEs), leading to the formation of dsRNAs that trigger RIG-I and MDA5 signaling that causes the constitutive secretion of IFN ɑ/β. As such, our data suggest that recognition of TE-derived dsRNAs by innate immune receptors promotes a pro-inflammatory environment favorable to the establishment of central tolerance in the thymic medulla.

At first sight, the pleiotropic effects of TEs on thymic function may look surprising. It should be reminded that the integration of genetic parasites such as TEs is a source of genetic conflicts with the host. Notably, the emergence of adaptive immunity gave rise to higher-order conflicts between TEs and their vertebrate hosts (Kassiotis, 2023; Boehm et al., 2023). The crucial challenge for the immune system is developing immune tolerance towards TEs to prevent autoimmune diseases that affect up to 10% of humans (Harroud and Hafler, 2023) without allowing selfish retrotransposition events that hinder genome integrity. The resolution of these conflicts has been proposed to be a determining factor in shaping the function of the immune system (Boehm et al., 2023). Our data suggest that the thymus is the central battlefield for conflict resolution between TEs and T cells in vertebrates. Consistent with the implication of TEs in autoimmunity, more than 90% of putative causal variants associated with autoimmune diseases are in allegedly noncoding regions of the genome (Harroud and Hafler, 2023). In this context, our study illustrates the complexity of interactions between TEs and the vertebrate immune system and should provide impetus to explore them further in health and disease. We see two limitations to our study. First, as with all multiomic systems immunology studies, our work provides a roadmap for many future mechanistic studies that could not be realized at this stage. Second, our immunopeptidomic analyses of TECs prove that TECs present TE MAPs but certainly underestimate the diversity of TE MAPs presented by cTECs and mTECs.

Materials and methods

Experimental design

Request a detailed protocol

This study aimed to understand better the impacts of TE expression on thymus development and function. Thymic populations are complex and heterogeneous, so we opted for single-cell RNA-seq data to draw a comprehensive profile of TE expression in the thymus. To better understand the impact of AIRE, FEZF2, and CHD4 on TE expression in the mTEC(II) population, RNA-seq data from WT and KO murine mTEC, as well as ChIP-seq for different histone marks in murine mTECs, were reanalyzed to characterize the TE sequences regulated by these three proteins. Unless stated otherwise, studies were done in human cells. For MS analyses, two replicates of 70 million cells from K5D1 mice (Robles et al., 1996) were injected for both cTECs and mTECs. All experiments in mice were in accordance with the Canadian Council on Animal Care guidelines and approved by the Comité de Déontologie de l’Expérimentation sur des Animaux of Université de Montréal. Primary human thymi were obtained from 4-month-old to 12-year-old children undergoing cardiovascular surgery at the CHU Sainte-Justine. This project was approved by the CHU Sainte-Justine Research Ethics Board (protocol and biobank #2126), which also granted permission to publish results. Written informed consent was obtained from parents of all children involved.

Transcriptomic data processing

Request a detailed protocol

Preprocessing of the scRNA-seq data was performed with kallisto (Bray et al., 2016), which uses an expectation-maximization algorithm to reassign multimapping reads based on the frequency of unique mappers at each sequence and bustools workflow. For human thymic data from Park et al., 2020 and splenic data from Madissoon et al., 2019, two different indexes were built for the pseudoalignment of reads with kallisto (version 0.46.0): one containing Ensembl 88 (GRCh38.88) transcripts used for the annotation of cell populations, and a second containing Ensembl 88 transcripts and human TE sequences (LINE, LTR, SINE) from RepeatMasker (Smit et al., 2013) which was used for all subsequent analyses of TE expression. For murine data from Baran-Gale et al., 2020, cell-type annotations from the original publication were used, and an index containing mm10 transcripts and murine TE sequences from RepeatMasker was used to analyze TE expression. The cell barcodes were corrected, and the feature-barcode matrices were generated with the correct count functions of bustools (version 0.39.3) (Melsted et al., 2019). For murine bulk RNA-seq data, an index composed of mm10 (GRCm38) transcripts and murine TE sequences from RepeatMasker was used for quantification with kallisto.

ChIP-seq data reanalysis

Request a detailed protocol

ChIP-seq data for (i) ETS1 in human NK cells, (ii) AIRE in murine mTECs, and (iii) several histone marks of mTECs from WT mice were reanalyzed (see ‘Data availability statement’ for the complete list). ETS1 ChIP-seq reads were aligned to the reference Homo sapiens genome (GRCh38) using bowtie2 (version 2.3.5) (Langmead and Salzberg, 2012) with the --very-sensitive parameter. Multimapping reads were removed using the samtools view function with the -q 10 parameter, and duplicate reads were removed using the samtools markdup function with the -r parameter (Danecek et al., 2021). Peak calling was performed with macs2 with the -m 5 50 parameter (Zhang et al., 2008). Peaks overlapping with the ENCODE blacklist regions (Amemiya et al., 2019) were removed with bedtools intersect (Quinlan and Hall, 2010) with default parameters. Overlap of ETS1 peaks with TE sequences was determined using bedtools intersect with default parameters. BigWig files were generated using the bamCoverage function of deeptools2 (Ramírez et al., 2016), and genomic tracks were visualized in the USCS Genome Browser (Kent et al., 2002). For the murine histone marks and AIRE data, reads were aligned to the reference Mus musculus genome (mm10) using bowtie2 with the --very-sensitive parameter. Multimapping reads were removed using the samtools view function with the -q 10 parameter, and duplicate reads were removed using the samtools markdup function with the -r parameter. For histone marks, read coverage at the sequence body and flanking regions (±3000 base pairs) of TE loci induced by AIRE, FEZF2, and CHD4 was visualized using ngs.plot.r (version 2.63) (Shen et al., 2014). For AIRE, peaks overlapping with the ENCODE blacklist regions were removed with bedtools intersect, and overlap of peaks with TE sequences was determined using bedtools intersect with default parameters.

Cell population annotation

Request a detailed protocol

Feature-barcode matrices were imported in R with SingleCellExperiment (version 1.12.0) (Amezquita et al., 2020). As a quality control, cells with less than 2000 UMI detected, less than 500 genes detected, or more than 5% reads assigned to mitochondrial genes were considered low quality and removed from the dataset with scuttle (version 1.0.4) (McCarthy et al., 2017). Cells with more than 7000 genes detected were considered doublets and removed. Normalization of cell size factors was performed with scran (version 1.18.7) (Lun et al., 2016), and log-normalization of read counts was done with scuttle with default parameters. Variable regions of TCR and IG genes, as well as ribosomal and cell cycle genes (based on Park et al., 2020), were removed, and highly variable features were selected based on a mean-variance trend based on a Poisson distribution of noise with scran. Adjustment of sequencing depths between batches and mutual nearest neighbors (MNN) correction were computed with batchelor (version 1.6.3) (Haghverdi et al., 2018). Cell clustering was performed with scran using the Jaccard index for edge weighting and the Louvain method for community detection. Lists of marker genes for human thymic cell populations and TEC subsets were taken from Park et al., 2020 and Bautista et al., 2021, whereas marker genes of splenic populations were based on Madissoon et al., 2019.

TE expression throughout thymic development

Request a detailed protocol

The expression of TE subfamilies was obtained by summing the read counts of loci based on the RepeatMasker annotations. For each TE subfamily in each cell population, expression levels amongst developmental stages were normalized by dividing them with the maximal expression value. Next, the Euclidean distance between each TE subfamily in each cell population (based on their normalized expression across developmental stages) was computed, followed by unsupervised hierarchical clustering. The tree was then manually cut into three clusters, and enrichment of LINE, LTR, and SINE elements in these three clusters was determined using Fisher’s exact tests. The cluster assigned to each TE subfamily in each cell population was visualized in a circos plot using the circlize package (version 0.4.14) (Gu et al., 2014) in R, and the percentage of each cell population found in embryonic or postnatal samples. Finally, we computed the frequency that each TE family was assigned to the three clusters, and the maximal value was kept. As a control, a random distribution of the expression of 809 TE subfamilies in 18 cell populations was generated. A cluster (cluster 1, 2, or 3) was randomly attributed for each combination of TE subfamily and cell type, and the maximal occurrence of a given cluster across cell types was then computed for each TE subfamily. Finally, the LINE, LTR, and SINE element distributions were compared to the random distribution with Kolmogorov–Smirnov tests.

Regulation of TE expression by cell proliferation and KZFPs

Request a detailed protocol

Proliferation scores were generated for each dataset cell using the CellCycleScoring function of Seurat (version 4.1.0). As per Cowan et al., 2019, we combined previously published lists of G2M and S phase marker genes (Kowalczyk et al., 2015) to compute the proliferation scores. For each thymic cell population, we calculated the Spearman correlation between proliferation scores and the expression of TE subfamilies. The Benjamini–Hochberg method was applied to correct for multiple comparisons. Correlations were considered positive if the correlation coefficient was ≥0.2 and the adjusted p-value≤0.05, and negative if the coefficient was ≤–0.2 and the adjusted p-value≤0.05. We also computed the median of all correlation coefficients for each cell population. We then assigned the class of each TE subfamily correlated with cell proliferation and compared this distribution to the distribution of classes of all TE subfamilies in the human genome. The percentage of overlap of the sets of TE subfamilies significantly correlated with cell proliferation was determined. A list of 401 human KZFPs was downloaded from Imbeault et al., 2017. Spearman correlations between KZFP and TE expression were independently computed in each cell population with the same methodology as the cell proliferation analysis, and Benjamini–Hochberg correction for multiple comparisons was applied. The information on the enrichment of KZFPs within TE subfamilies was downloaded from Imbeault et al., 2017. Sharing of KZFP-TE pairs between cell populations was represented using the circlize package.

Estimation of TE sequences’ age

Request a detailed protocol

The sequence divergence (defined as the number of mismatches per thousand) was given by the milliDiv value in RepeatMasker. The milliDiv values of each TE locus were divided by the substitution rate of its host’s genome (2.2 × 10–9 mutation/year for Homo sapiens and 4.5 × 10–9 mutation/year for Mus musculus; Lander et al., 2001; Waterston et al., 2002). Finally, the age of each TE subfamily was determined by averaging the age of all loci of the subfamily.

Interactions between TE subfamilies and transcription factors

Request a detailed protocol

We downloaded a list of 1638 transcription factors (TFs) manually curated by Lambert et al., 2018. For each cell population of the thymus, Spearman correlations were computed for each possible pair of TF and TE subfamily, and the Benjamini–Hochberg method was applied to correct the p-values for multiple comparisons. Correlations were considered significant if (i) the correlation coefficient was ≥0.2, (ii) the adjusted p-value was ≤0.05, and (iii) the TF was expressed by ≥10% of the cells of the population. The correlations were validated using a bootstrap procedure (1000 iterations) to ensure their reproducibility. Briefly, we randomly selected n cells out of the n cells of a given population (while allowing cells to be selected multiple times). The empirical p-value was determined by dividing the number of iterations with a correlation coefficient <0.2 by the total number of iterations (1000). In parallel, the curated binding motifs of 945 TFs were downloaded from the JASPAR database. We then used the Find Individual Motif Occurrences (FIMO) software (Grant et al., 2011) to identify the 100,000 genomic positions with the most significant matches for the TF binding motif. These lists of binding motif positions were then intersected with the positions of TE loci with the intersect function of BEDTools (version 2.29.2) (Quinlan and Hall, 2010), and the percentage of TE loci of each subfamily harboring TF binding motifs was determined. Thus, in a specific cell population of the thymus, a TF was considered as interacting with a TE subfamily if it satisfied two criteria: (i) its expression was correlated with the one of the TE family (Spearman coefficient ≥ 0.2, adjusted p-value≤0.05, and expression of TF in ≥10% of cells), and (ii) at least one locus of the TE subfamily contained a binding motif of the TF. For each cell population, networks of interactions between TF and TE subfamilies were generated with the network package (version 1.17.1) (Butts, 2008) in R and represented with the ggnetwork package. For the sake of clarity, only the most significant interactions were illustrated for each cell type (i.e., correlation coefficient ≥ 0.3, TF binding sites in ≥1% of the loci of the TE subfamily, and TF expression in ≥10% of cells of the population). Sharing of TF-TE interactions between cell populations was represented with a chord diagram using the circlize package. For each TE subfamily, the number of interactions with TFs and the number of loci of the TE subfamily in the human genome were determined. Wilcoxon–Mann–Whitney tests were used to compare the number of interactions with TF of LTR, LINE, and SINE elements, whereas Kendall tau correlation was calculated between the number of interactions with TF and the number of loci of TE subfamilies.

Identification of TE promoter and enhancer candidates

Request a detailed protocol

From the previously identified list of TF-TE interactions, we isolated the specific loci containing TF binding sites from the subfamilies whose expression was positively correlated with the TF. To determine if these TE loci could act as promoters or enhancers, we used histone ChIP-seq data from the ENCODE consortium for H3K27ac, H3K4me1, and H3K4me3. BED files from the ENCODE consortium were downloaded for eight immune cell populations: B cells, CD4 single positive T cells (CD4 SP), CD8 single positive T cells (CD8 SP), DCs, monocytes and macrophages (mono/macro), NK cells, Th17, and Treg. TE loci colocalizing with peaks in histone ChIP-seq data were identified using the intersect function of BEDTools (version 2.29.2). To be considered enhancer candidates, TE loci had to colocalize with H3K27ac and H3K4me1 but not H3K4me3. To be considered as promoter candidates, TE loci had to colocalize with H3K27ac and H3K4me3, but not H3K4me1, and be located at ≤1000 nucleotides from a transcription start site (TSS) annotated in the refTSS database (Abugessaisa et al., 2019).

Diversity of TE expression

Request a detailed protocol

The human thymic scRNA-seq dataset was subsampled to retain only postnatal cells, as it was shown by Taveirne et al., 2020 that thymic APCs are mainly found in postnatal samples. The diversity of TE sequences expressed by thymic populations was assessed using Shannon entropy. Using the vegan package (version 2.5–7) (Tikhonov et al., 2020) in R, two distinct Shannon entropy metrics were computed for each cell population. First, the Shannon entropy was computed based on the expression level (i.e., log(read count)) of TE subfamilies for each cell individually. The median entropy was calculated for each cell population. In parallel, the diversity of TE sequences expressed by an entire population was also assessed. For this purpose, a binary code was generated to represent the expression status of TE subfamilies in each cell (where 1 is expressed and 0 is not expressed). For each population separately, the binary codes of individual cells were summed to obtain the frequency of expression of each TE subfamily in the population, which was used to compute the Shannon entropy of TE sequences expressed by the population. A linear model was generated with the lm function of the stats package in R to summarize the data distribution. The deviation (Δy) from the observed population’s TE diversity and the one expected by the linear model was computed for each cell population.

TE expression in thymic APC

Request a detailed protocol

A differential expression analysis of TE subfamilies between the subsets of thymic APC was performed with the FindAllMarkers function with default parameters of Seurat (Satija et al., 2015) with the MAST model. Finally, the heterogeneity of TE expression inside thymic APC subsets was evaluated with the MetaCell package (version 0.3.5) (Baran et al., 2019). The composition of the metacells was validated based on manual annotation (see the ‘Single-cell RNA-seq preprocessing’ section), and only metacells with >50% of cells belonging to the same subset of thymic APCs were kept. Differential expression of TE subfamilies between metacells was performed as described above, and the percentage of overlap between the sets of TEs overexpressed by the different metacells was computed.

Isolation of human thymic pDCs and immunostaining of dsRNAs

Request a detailed protocol

Primary human thymi were obtained from 4-month-old to 12-year-old children undergoing cardiovascular surgeries at the CHU Sainte-Justine. This project was approved by the CHU Sainte-Justine Research Ethics Board (protocol and biobank #2126). Thymi from 4-month-old to 12-year-old individuals were cryopreserved in liquid nitrogen in the following solution: 95% (PBS-5% Dextran 40 (Sigma-Aldrich)) – 5% DMSO (Fisher Scientific). Protocol for thymic pDCs isolation was based on Stoeckle et al., 2013. Briefly, thymic samples were cut in ~2 mm pieces, followed by three rounds of digestion (40 min, 180 RPM at 37℃) in RPMI 1640 (Gibco) supplemented with 2 mg/ml of Collagenase A (Roche) and 0.1 mg/ml of DNase I (Sigma-Aldrich). APCs were then enriched using Percoll (Sigma-Aldrich) density centrifugation (3500 × g, 35 min at 4℃), followed by an FBS cushion density gradient (5 ml of RPMI 1640 containing enriched APCs layered on 5 ml of heat-inactivated FBS [Invitrogen, 12483020], 1000 RPM for 10 min at 4℃) to remove cell debris. Finally, thymic pDCs were magnetically enriched using the QuadroMACS Separator (Miltenyi). Cells were stained with a CD303 (BDCA-2) MicroBead Kit (Miltenyi), and labeled cells were loaded on LS columns (Miltenyi) for magnetic-activated cell sorting.

Purified thymic pDCs were pipetted on poly-l-lysine (Sigma-Aldrich, 1:10 in dH2O) coated 15μ-Slide 8 well (ibid) and incubated for 2 hr at 37℃ in RPMI 1640 supplemented with 10% BSA (Sigma-Aldrich). Cells were fixed using 1% (w/v) paraformaldehyde (PFA, Sigma-Aldrich) in PBS 1× (Sigma-Aldrich) for 30 min at room temperature. Cells were permeabilized for 30 min at room temperature with 0.1% (v/v) Triton X-100 (Sigma-Aldrich) in PBS 1×, followed by blocking using 5% (w/v) BSA (Sigma-Aldrich) in PBS 1× for 30 min at room temperature. Immunostaining was performed in four steps to avoid unspecific binding of the secondary antibodies: (i) incubation overnight at 4℃ with the mouse monoclonal IgG2a J2 antibody anti-dsRNA (Jena Bioscience, Cat# RNT-SCI-10010500, dilution 1:200), (ii) incubation with the donkey anti-mouse IgG (H+L) antibody coupled to Alexa Fluor 555 (Invitrogen, Cat# A-31570, dilution 1:500) for 30 min at room temperature, (iii) incubation with the mouse monoclonal IgG1 clone 6H6 anti-CD123 (eBioscience, Cat# 14-1239-82, 1:100) for 1 hr at room temperature, and (iv) incubation with the goat anti-mouse IgG1 polyclonal Alexa Fluor 488 antibody (Invitrogen, Cat# A-21121, 1:1000) for 30 min at room temperature. Finally, cells were stained with DAPI (Invitrogen, Cat# D3571, 1:1000) for 5 min at room temperature. All antibodies and DAPI were diluted in a blocking solution. Image acquisition was made with an LSM 700 laser scanning confocal microscope (Zeiss) using a ×40 oil objective (Zeiss, Plan-Neofluar N.A. 1.4) and the ZEN software. Using the whiteTopHat function of the EBImage package and the sigmoNormalize function of the MorphoR package in R, the background of the DAPI signal was removed. The nuclei were segmented on the resulting images as circular shapes based on the DAPI signal. The mean intensity of CD123 and J2 staining was determined for each cytoplasm, defined as 19 nm rings around nuclei. Based on the distribution of the CD123 signal across cells, a threshold between CD123- and CD123+ cells was set up for each replicate independently. J2 signal intensity was compared between CD123- and CD123+ cells using the Wilcoxon rank-sum test in R.

Gene set enrichment analysis

Request a detailed protocol

Gene set enrichment analyses were performed to determine which biological processes are enriched in mTEC(II) and pDCs. Differential gene expression analyses were performed between each possible pair of thymic APCs subsets using MAST with the FindMarkers function of Seurat. The gene set enrichment analysis was performed using the iDEA package (version 1.0.1) (Ma et al., 2020) in R. As per Ma et al., 2020, the fold change and standard error of gene expression were used as input for iDEA, in addition to predefined lists of gene sets compiled in the iDEA package. Gene sets associated with antigen presentation, interferon signaling, and immune response were manually annotated. iDEA was launched with default parameters, except for the 500 iterations of the Markov chain Monte Carlo algorithm, and p-values were corrected with the Louis method. We also visualized the expression of AIRE, FEZF2, and CHD4 in the TEC lineage to validate their expression in mTEC(II).

TE loci regulated by AIRE, FEZF2, and CHD4

Request a detailed protocol

A differential expression analysis of TE subfamilies between WT and Aire-, Fezf2-, or Chd4-KO mice was performed with the voom method of the limma package (version 3.46.0) (Law et al., 2014; Ritchie et al., 2015). Stringent criteria (i.e., an expression below 2 transcripts per million [TPM] in all samples) were applied to remove lowly expressed TEs. TE subfamilies with (i) a fold change ≥2 and an adjusted p-value≤0.05 or (ii) a fold change ≤–2 and an adjusted p-value≤0.05 were considered as induced and repressed, respectively. The percentage of overlap between the sets of TE loci induced or repressed by AIRE, FEZF2, and CHD4 was computed. The class and subfamily were assigned to each regulated TE locus, and the distributions of classes and subfamilies across all TE sequences of the murine genome were used as controls. Significant enrichment of classes or subfamilies was determined with chi-squared tests, and a Bonferroni correction for multiple comparisons was performed to enrich subfamilies in induced or repressed TEs. The distance between TE loci induced or repressed by AIRE, FEZF2, or CHD4 was defined as the minimal distance between the middle position of TE loci on the same chromosome. As a control, distributions of randomly selected TE loci whose expression is independent of AIRE, FEZF2, and CHD4 and equal size to the sets of regulated TEs were generated (e.g., if 433 TE loci are induced by CHD4, 433 independent TE loci were randomly selected). Wilcoxon rank-sum tests were used to compare random and regulated distributions. Genomic positions of exons, introns 3′ and 5′ untranslated transcribed region (UTR) were downloaded from the UCSC Table Browser. The genomic localization of regulated TEs was determined using the intersect mode of the BEDTools suite version 2.29.2. TE loci not located in exons, introns, 3′UTR, or 5′UTR were considered intergenic. The percentage of regulated TE loci in each type of genomic region was determined and compared to the genomic localization of all TE loci in the murine genome with chi-square tests. Finally, we estimated the frequency of intron retention events for introns containing TE loci regulated by AIRE, FEZF2, or CHD4 with S-IRFindeR (Broseus and Ritchie, 2020). Sequencing reads were aligned to the reference Mus musculus genome (mm10) using STAR version 2.7.1a (Dobin et al., 2013) with default parameters. Each intron’s Stable Intron Retention ratio (SIRratio) was computed with the computeSIRratio function of S-IRFindeR. Introns containing TE loci induced by AIRE, FEZF2, or CHD4 were filtered using BEDTools intersect. Random distributions of equivalent sizes of introns containing TE sequences independent of AIRE, FEZF2, and CHD4 were generated as control. A SIRratio of 0.1 was used as a threshold of significant intron retention events.

Enzymatic digestion and isolation of murine TECs

Request a detailed protocol

Thymic stromal cell enrichment was performed as previously described (Dumont-Lagacé et al., 2019; Kim et al., 2015). Briefly, thymi from 16- to 22-week-old K5D1 mice were mechanically disrupted and enzymatically digested with papain (Worthington Biochemical Corporation), DNase I (Sigma-Aldrich), and collagenase IV (Sigma-Aldrich) at 37°C. Next, the single-cell suspension obtained after enzymatic digestion was maintained at 4°C in FACS buffer (PBS, 0.5% [w/v] BSA, 2 mM EDTA) and enriched in thymic epithelial cells using anti-EpCAM (CD326) or anti-CD45 microbeads (mouse, Miltenyi) and LS columns (Miltenyi). Then, the enriched epithelial cell suspension was stained for flow cytometry cell sorting with the following antibodies and dyes: anti-EpCAM-APC-Cy7 clone G8.8 (BioLegend, Cat# 118218), anti-CD45-APC clone 30-F11 (BD Biosciences, Cat# 559864), anti-UEA1-biotinylated (Vector Laboratories, Cat# B-1065), anti-I-A/I-E-Alexa Fluor 700 clone M5/114.15.2 (BioLegend, Cat# 107622), anti-Ly51-FITC clone 6C3 (BioLegend, Cat# 553160), anti-streptavidin-PE-Cy7 (BD Biosciences, Cat# 557598), and 7-AAD (BD Biosciences, Cat# 559925). Cell sorting was performed using a BD FACSAria (BD Biosciences), and data were analyzed using the FACSDiva. TECs were defined as EpCAM+CD45−, while the cTEC and mTEC subsets were defined as UEA1–Ly51+ and UEA1+Ly51- TEC, respectively.

RNA-sequencing

Request a detailed protocol

Total RNA from 80,000 mTECs or cTECs was isolated using TRIzol and purified with an RNeasy micro kit (QIAGEN). Total RNA was quantified using Qubit (Thermo Scientific), and RNA quality was assessed with the Agilent 2100 Bioanalyzer (Agilent Technologies). Transcriptome libraries were generated using a KAPA RNA HyperPrep kit (Roche) using a poly(A) selection (Thermo Scientific). Sequencing was performed on the Illumina NextSeq 500, obtaining ∼200 million paired-end reads per sample.

Preparation of CNBR-activated Sepharose beads for MHC I immunoprecipitation

Request a detailed protocol

CNBR-activated Sepharose 4B beads (Sigma-Aldrich, Cat# 17-0430-01) were incubated with 1 mM HCl at a ratio of 40 mg of beads per 13.5 ml of 1 mM HCl for 30 min with tumbling at room temperature. Beads were spun at 215 × g for 1 min at 4°C, and supernatants were discarded. 40 mg of beads were resuspended with 4 ml of coupling buffer (0.1 M NaHC03/0.5 M NaCl pH 8.3), spun at 215 × g for 1 min at 4°C, and the supernatants were discarded. Mouse antibodies Pan-H2 (clone M1/42), H2-Kb (clone Y-3), and H2-Db (clone 28-14-8S) were coupled to beads at a ratio of 1 mg of antibody to 40 mg of beads in coupling buffer for 120 min with tumbling at room temperature. Beads were spun at 215 × g for 1 min at 4°C, and supernatants were discarded. 40 mg of beads were resuspended with 1 ml of blocking buffer (0.2 M glycine), incubated for 30 min with tumbling at room temperature, and the supernatants were discarded. Beads were washed by centrifugation twice with PBS pH 7.2, resuspended at a concentration of 1 mg of antibody per ml of PBS pH 7.2, and stored at 4°C.

Immuno-isolation of MAPs

Request a detailed protocol

Frozen pellets of mTECs (90 mice, 191 million cells total) and cTECs (121 mice, 164 million cells total) were thawed, pooled, and resuspended with PBS pH 7.2 up to 4 ml and then solubilized by adding 4 ml of detergent buffer containing PBS pH 7.2, 1% (w/v) CHAPS (Sigma, Cat# C9426-5G) supplemented with Protease inhibitor cocktail (Sigma, Cat# P8340-5mL). Solubilized cells were incubated for 60 min with tumbling at 4°C and then spun at 16,600 × g for 20 min at 4°C. Supernatants were transferred into new tubes containing 1.5 mg of Pan-H2, 0.5 mg of H2-Kb, and 0.5 mg of H2-Db antibodies covalently-cross-linked CNBR-Sepharose beads per sample and incubated with tumbling for 180 min at 4°C. Samples were transferred into Bio-Rad Poly prep chromatography columns and eluted by gravity. Beads were first washed with 11.5 ml PBS, then with 11.5 ml of 0.1× PBS, and finally with 11.5 ml of water. MHC I complexes were eluted from the beads by acidic treatment using 1% trifluoroacetic acid (TFA). Acidic filtrates containing peptides were separated from MHC I subunits (HLA molecules and β–2 microglobulin) using home-made stage tips packed with two 1 mm diameter octadecyl (C-18) solid phase extraction disks (EMPORE). Stage tips were pre-washed with methanol, then with 80% acetonitrile (ACN) in 0.1% TFA, and finally with 1% TFA. Samples were loaded onto the stage tips and washed with 1% TFA and 0.1% TFA. Peptides were eluted with 30% ACN in 0.1% TFA, dried using vacuum centrifugation, and then stored at –20°C until MS analysis.

MS analyses

Request a detailed protocol

Peptides were loaded and separated on a home-made reversed-phase column (150 μm i.d. by 200 mm) with a 106 min gradient from 10 to 38% B (A: formic acid 0.1%; B: 80% CAN 0.1% formic acid) and a 600 nl/min flow rate on an Easy nLC-1200 connected to an Orbitrap Exploris 480 (Thermo Fisher Scientific). Each full MS spectrum acquired at a resolution of 240,000 was followed by tandem-MS (MS-MS) spectra acquisition on the most abundant multiply charged precursor ions for a maximum of 3 s. Tandem-MS experiments were performed using higher energy collision-induced dissociation (HCD) at a collision energy of 34%. The generation of the personalized proteome containing TE sequences, as well as the identification of TE-derived MAPs, was performed as per Larouche et al., 2020 with the following modifications: the mm10 murine reference genome was downloaded from the UCSC Genome Browser, the annotations for murine genes and TE sequences were downloaded from the UCSC Table Browser, and the UniProt mouse database (16,977 entries) was used for the canonical proteome. MAPs were identified using PEAKS X Pro (Bioinformatics Solutions, Waterloo, ON). The level and frequency of expression of TE subfamilies generating MAPs or not were determined in thymic epithelial cells were determined by averaging the expression values across cells of a TEC subset and dividing the number of cells with a positive (i.e., >0) expression of the TEs by the total number of cells of the TEC subset, respectively.

Data availability

scRNA-seq data of human thymi and spleen were downloaded from ArrayExpress (accession number E-MTAB-8581) and the NCBI BIOPROJECT (accession code PRJEB31843), respectively. scRNA-seq data of murine thymi were downloaded from ArrayExpress (accession number E-MTAB-8560). RNA-seq data from WT, Aire-KO, Fezf2-KO, and Chd4-KO murine mTECs were downloaded from the Gene Expression Omnibus (GEO) under the accession code GSE144880. ChIP-seq data of ETS1 in human NK cells and AIRE in murine mTECs were downloaded from GEO (accession codes GSE124104 and GSE92654, respectively). ChIP-seq data for different histone marks in murine mTECs were also downloaded from GEO: H3K4me3 for mTECs (GSE53111); H3K4me1 and H3K27ac from MHCIIhi mTECs (GSE92597); H3K4me2 in mTEC-II (GSE103969); and H3K4ac and H3K9ac in mTECs (GSE114713). Transcriptomic and immunopeptidomic data of K5D1 mice mTECs and cTECs generated in this study are available on the Gene Expression Omnibus (GEO) under the accession GSE232011 and on the Proteomics Identification Database (PRIDE) under the accession PXD042241, respectively.

The following data sets were generated
    1. Courcelles M
    2. Thibault P
    (2024) PRIDE
    ID PXD042241. LC-MSMS of the MHC-I immunopeptidome of murine thymic epithelial cells.
    https://www.ebi.ac.uk/pride/archive/projects/PXD042241
    1. Larouche JD
    2. Perreault C
    (2024) NCBI Gene Expression Omnibus
    ID GSE232011. Transposable elements regulate T cell maturation and education in the thymus.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE232011
The following previously published data sets were used
    1. Park JE
    (2020) EMBL-EBI ArrayExpress
    ID E-MTAB-8581. A cell atlas of human thymic development defines T cell repertoire formation.
    https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-8581?query=E-MTAB-8581
    1. Madissoon E
    2. Meyer KB
    (2019) NCBI BioProject
    ID PRJEB31843. Ischaemic sensitivity of human tissue by single cell RNA seq.
    https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB31843
    1. Baran-Gale J
    2. Morgan M
    (2020) EMBL-EBI ArrayExpress
    ID E-MTAB-8560. Single-cell RNA-sequencing of mouse thymic epithelial cells across the first year of life.
    https://www.ebi.ac.uk/biostudies/arrayexpress/studies/E-MTAB-8560?query=E-MTAB-8560
    1. Bansal K
    2. Yoshida H
    3. Benoist C
    4. Mathis D
    (2017) NCBI Gene Expression Omnibus
    ID GSE92654. Aire, guardian of immunological tolerance, binds to and activates super-enhancers.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92654
    1. Sansom SN
    2. Shikama-Dorn N
    3. Zhanybekova S
    4. Nusspaumer G
    5. Heger A
    6. Ponting CP
    7. Hollander GA
    (2014) NCBI Gene Expression Omnibus
    ID GSE53111. Aire secures the transcription of polycomb marked genes to complete a comprehensive program of promiscuous gene expression in the thymic epithelial cell lineage.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE53111
    1. Bansal K
    2. Mathis D
    3. Benoist C
    (2017) NCBI Gene Expression Omnibus
    ID GSE92597. Aire, guardian of immunological tolerance, binds to and activates super-enhancers [ChIP-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE92597
    1. Handel AE
    2. Holländer GA
    (2018) NCBI Gene Expression Omnibus
    ID GSE114713. Comprehensively profiling the chromatin architecture of tissue restricted antigen expression in thymic epithelial cells.
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE114713
    1. Amit I
    2. Abramson J
    3. Bornstein C
    4. Nevo S
    5. Giladi A
    6. Kadouri N
    (2018) NCBI Gene Expression Omnibus
    ID GSE103969. Large-scale single cell mapping of the thymic stroma identifies a new thymic epithelial cell lineage [ChIP-seq].
    https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103969

References

    1. Lander ES
    2. Linton LM
    3. Birren B
    4. Nusbaum C
    5. Zody MC
    6. Baldwin J
    7. Devon K
    8. Dewar K
    9. Doyle M
    10. FitzHugh W
    11. Funke R
    12. Gage D
    13. Harris K
    14. Heaford A
    15. Howland J
    16. Kann L
    17. Lehoczky J
    18. LeVine R
    19. McEwan P
    20. McKernan K
    21. Meldrim J
    22. Mesirov JP
    23. Miranda C
    24. Morris W
    25. Naylor J
    26. Raymond C
    27. Rosetti M
    28. Santos R
    29. Sheridan A
    30. Sougnez C
    31. Stange-Thomann Y
    32. Stojanovic N
    33. Subramanian A
    34. Wyman D
    35. Rogers J
    36. Sulston J
    37. Ainscough R
    38. Beck S
    39. Bentley D
    40. Burton J
    41. Clee C
    42. Carter N
    43. Coulson A
    44. Deadman R
    45. Deloukas P
    46. Dunham A
    47. Dunham I
    48. Durbin R
    49. French L
    50. Grafham D
    51. Gregory S
    52. Hubbard T
    53. Humphray S
    54. Hunt A
    55. Jones M
    56. Lloyd C
    57. McMurray A
    58. Matthews L
    59. Mercer S
    60. Milne S
    61. Mullikin JC
    62. Mungall A
    63. Plumb R
    64. Ross M
    65. Shownkeen R
    66. Sims S
    67. Waterston RH
    68. Wilson RK
    69. Hillier LW
    70. McPherson JD
    71. Marra MA
    72. Mardis ER
    73. Fulton LA
    74. Chinwalla AT
    75. Pepin KH
    76. Gish WR
    77. Chissoe SL
    78. Wendl MC
    79. Delehaunty KD
    80. Miner TL
    81. Delehaunty A
    82. Kramer JB
    83. Cook LL
    84. Fulton RS
    85. Johnson DL
    86. Minx PJ
    87. Clifton SW
    88. Hawkins T
    89. Branscomb E
    90. Predki P
    91. Richardson P
    92. Wenning S
    93. Slezak T
    94. Doggett N
    95. Cheng JF
    96. Olsen A
    97. Lucas S
    98. Elkin C
    99. Uberbacher E
    100. Frazier M
    101. Gibbs RA
    102. Muzny DM
    103. Scherer SE
    104. Bouck JB
    105. Sodergren EJ
    106. Worley KC
    107. Rives CM
    108. Gorrell JH
    109. Metzker ML
    110. Naylor SL
    111. Kucherlapati RS
    112. Nelson DL
    113. Weinstock GM
    114. Sakaki Y
    115. Fujiyama A
    116. Hattori M
    117. Yada T
    118. Toyoda A
    119. Itoh T
    120. Kawagoe C
    121. Watanabe H
    122. Totoki Y
    123. Taylor T
    124. Weissenbach J
    125. Heilig R
    126. Saurin W
    127. Artiguenave F
    128. Brottier P
    129. Bruls T
    130. Pelletier E
    131. Robert C
    132. Wincker P
    133. Smith DR
    134. Doucette-Stamm L
    135. Rubenfield M
    136. Weinstock K
    137. Lee HM
    138. Dubois J
    139. Rosenthal A
    140. Platzer M
    141. Nyakatura G
    142. Taudien S
    143. Rump A
    144. Yang H
    145. Yu J
    146. Wang J
    147. Huang G
    148. Gu J
    149. Hood L
    150. Rowen L
    151. Madan A
    152. Qin S
    153. Davis RW
    154. Federspiel NA
    155. Abola AP
    156. Proctor MJ
    157. Myers RM
    158. Schmutz J
    159. Dickson M
    160. Grimwood J
    161. Cox DR
    162. Olson MV
    163. Kaul R
    164. Raymond C
    165. Shimizu N
    166. Kawasaki K
    167. Minoshima S
    168. Evans GA
    169. Athanasiou M
    170. Schultz R
    171. Roe BA
    172. Chen F
    173. Pan H
    174. Ramser J
    175. Lehrach H
    176. Reinhardt R
    177. McCombie WR
    178. de la Bastide M
    179. Dedhia N
    180. Blöcker H
    181. Hornischer K
    182. Nordsiek G
    183. Agarwala R
    184. Aravind L
    185. Bailey JA
    186. Bateman A
    187. Batzoglou S
    188. Birney E
    189. Bork P
    190. Brown DG
    191. Burge CB
    192. Cerutti L
    193. Chen HC
    194. Church D
    195. Clamp M
    196. Copley RR
    197. Doerks T
    198. Eddy SR
    199. Eichler EE
    200. Furey TS
    201. Galagan J
    202. Gilbert JG
    203. Harmon C
    204. Hayashizaki Y
    205. Haussler D
    206. Hermjakob H
    207. Hokamp K
    208. Jang W
    209. Johnson LS
    210. Jones TA
    211. Kasif S
    212. Kaspryzk A
    213. Kennedy S
    214. Kent WJ
    215. Kitts P
    216. Koonin EV
    217. Korf I
    218. Kulp D
    219. Lancet D
    220. Lowe TM
    221. McLysaght A
    222. Mikkelsen T
    223. Moran JV
    224. Mulder N
    225. Pollara VJ
    226. Ponting CP
    227. Schuler G
    228. Schultz J
    229. Slater G
    230. Smit AF
    231. Stupka E
    232. Szustakowki J
    233. Thierry-Mieg D
    234. Thierry-Mieg J
    235. Wagner L
    236. Wallis J
    237. Wheeler R
    238. Williams A
    239. Wolf YI
    240. Wolfe KH
    241. Yang SP
    242. Yeh RF
    243. Collins F
    244. Guyer MS
    245. Peterson J
    246. Felsenfeld A
    247. Wetterstrand KA
    248. Patrinos A
    249. Morgan MJ
    250. de Jong P
    251. Catanese JJ
    252. Osoegawa K
    253. Shizuya H
    254. Choi S
    255. Chen YJ
    256. Szustakowki J
    (2001) Initial sequencing and analysis of the human genome
    Nature 409:860–921.
    https://doi.org/10.1038/35057062
    1. Waterston RH
    2. Lindblad-Toh K
    3. Birney E
    4. Rogers J
    5. Abril JF
    6. Agarwal P
    7. Agarwala R
    8. Ainscough R
    9. Alexandersson M
    10. An P
    11. Antonarakis SE
    12. Attwood J
    13. Baertsch R
    14. Bailey J
    15. Barlow K
    16. Beck S
    17. Berry E
    18. Birren B
    19. Bloom T
    20. Bork P
    21. Botcherby M
    22. Bray N
    23. Brent MR
    24. Brown DG
    25. Brown SD
    26. Bult C
    27. Burton J
    28. Butler J
    29. Campbell RD
    30. Carninci P
    31. Cawley S
    32. Chiaromonte F
    33. Chinwalla AT
    34. Church DM
    35. Clamp M
    36. Clee C
    37. Collins FS
    38. Cook LL
    39. Copley RR
    40. Coulson A
    41. Couronne O
    42. Cuff J
    43. Curwen V
    44. Cutts T
    45. Daly M
    46. David R
    47. Davies J
    48. Delehaunty KD
    49. Deri J
    50. Dermitzakis ET
    51. Dewey C
    52. Dickens NJ
    53. Diekhans M
    54. Dodge S
    55. Dubchak I
    56. Dunn DM
    57. Eddy SR
    58. Elnitski L
    59. Emes RD
    60. Eswara P
    61. Eyras E
    62. Felsenfeld A
    63. Fewell GA
    64. Flicek P
    65. Foley K
    66. Frankel WN
    67. Fulton LA
    68. Fulton RS
    69. Furey TS
    70. Gage D
    71. Gibbs RA
    72. Glusman G
    73. Gnerre S
    74. Goldman N
    75. Goodstadt L
    76. Grafham D
    77. Graves TA
    78. Green ED
    79. Gregory S
    80. Guigó R
    81. Guyer M
    82. Hardison RC
    83. Haussler D
    84. Hayashizaki Y
    85. Hillier LW
    86. Hinrichs A
    87. Hlavina W
    88. Holzer T
    89. Hsu F
    90. Hua A
    91. Hubbard T
    92. Hunt A
    93. Jackson I
    94. Jaffe DB
    95. Johnson LS
    96. Jones M
    97. Jones TA
    98. Joy A
    99. Kamal M
    100. Karlsson EK
    101. Karolchik D
    102. Kasprzyk A
    103. Kawai J
    104. Keibler E
    105. Kells C
    106. Kent WJ
    107. Kirby A
    108. Kolbe DL
    109. Korf I
    110. Kucherlapati RS
    111. Kulbokas EJ
    112. Kulp D
    113. Landers T
    114. Leger JP
    115. Leonard S
    116. Letunic I
    117. Levine R
    118. Li J
    119. Li M
    120. Lloyd C
    121. Lucas S
    122. Ma B
    123. Maglott DR
    124. Mardis ER
    125. Matthews L
    126. Mauceli E
    127. Mayer JH
    128. McCarthy M
    129. McCombie WR
    130. McLaren S
    131. McLay K
    132. McPherson JD
    133. Meldrim J
    134. Meredith B
    135. Mesirov JP
    136. Miller W
    137. Miner TL
    138. Mongin E
    139. Montgomery KT
    140. Morgan M
    141. Mott R
    142. Mullikin JC
    143. Muzny DM
    144. Nash WE
    145. Nelson JO
    146. Nhan MN
    147. Nicol R
    148. Ning Z
    149. Nusbaum C
    150. O’Connor MJ
    151. Okazaki Y
    152. Oliver K
    153. Overton-Larty E
    154. Pachter L
    155. Parra G
    156. Pepin KH
    157. Peterson J
    158. Pevzner P
    159. Plumb R
    160. Pohl CS
    161. Poliakov A
    162. Ponce TC
    163. Ponting CP
    164. Potter S
    165. Quail M
    166. Reymond A
    167. Roe BA
    168. Roskin KM
    169. Rubin EM
    170. Rust AG
    171. Santos R
    172. Sapojnikov V
    173. Schultz B
    174. Schultz J
    175. Schwartz MS
    176. Schwartz S
    177. Scott C
    178. Seaman S
    179. Searle S
    180. Sharpe T
    181. Sheridan A
    182. Shownkeen R
    183. Sims S
    184. Singer JB
    185. Slater G
    186. Smit A
    187. Smith DR
    188. Spencer B
    189. Stabenau A
    190. Stange-Thomann N
    191. Sugnet C
    192. Suyama M
    193. Tesler G
    194. Thompson J
    195. Torrents D
    196. Trevaskis E
    197. Tromp J
    198. Ucla C
    199. Ureta-Vidal A
    200. Vinson JP
    201. Von Niederhausern AC
    202. Wade CM
    203. Wall M
    204. Weber RJ
    205. Weiss RB
    206. Wendl MC
    207. West AP
    208. Wetterstrand K
    209. Wheeler R
    210. Whelan S
    211. Wierzbowski J
    212. Willey D
    213. Williams S
    214. Wilson RK
    215. Winter E
    216. Worley KC
    217. Wyman D
    218. Yang S
    219. Yang SP
    220. Zdobnov EM
    221. Zody MC
    222. Lander ES
    (2002) Initial sequencing and comparative analysis of the mouse genome
    Nature 420:520–562.
    https://doi.org/10.1038/nature01262

Article and author information

Author details

  1. Jean-David Larouche

    1. Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    2. Department of Medicine, Université de Montréal, Montréal, Canada
    Contribution
    Conceptualization, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-3056-0714

  2. Céline M Laumont

    1. Deeley Research Centre, BC Cancer, Victoria, Canada
    2. Department of Medical Genetics, University of British Columbia, Vancouver, Canada
    Contribution
    Conceptualization, Visualization, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared

  3. Assya Trofimov

    1. Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    2. Department of Computer Science and Operations Research, Université de Montréal, Montréal, Canada
    3. Fred Hutchinson Cancer Center, Seattle, United States
    4. Department of Physics, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Visualization, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8748-3735

  4. Krystel Vincent

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Conceptualization, Methodology, Writing - original draft, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared

  5. Leslie Hesnard

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  6. Sylvie Brochu

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  7. Caroline Côté

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation
    Competing interests
    No competing interests declared

  8. Juliette F Humeau

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  9. Éric Bonneil

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  10. Joel Lanoix

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  11. Chantal Durette

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  12. Patrick Gendron

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared

  13. Jean-Philippe Laverdure

    Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    Contribution
    Investigation, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9245-6413

  14. Ellen R Richie

    Department of Epigenetics and Molecular Carcinogenesis, University of Texas M.D. Anderson Cancer Center, Houston, United States
    Contribution
    Resources, Writing – review and editing
    Competing interests
    No competing interests declared

  15. Sébastien Lemieux

    1. Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    2. Department of Biochemistry and Molecular Medicine, Université de Montréal, Montreal, Canada
    Contribution
    Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  16. Pierre Thibault

    1. Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    2. Department of Chemistry, Université de Montréal, Montréal, Canada
    Contribution
    Supervision, Writing – review and editing
    Competing interests
    No competing interests declared

  17. Claude Perreault

    1. Institute for Research in Immunology and Cancer, Université de Montréal, Montreal, Canada
    2. Department of Medicine, Université de Montréal, Montréal, Canada
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing – review and editing
    For correspondence
    claude.perreault@umontreal.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9453-7383

Funding

Canadian Institutes of Health Research (FDN 148400)

  • Claude Perreault

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Christian Charbonneau and Raphaëlle Lambert from IRIC’s bio-imaging and genomics platforms, respectively. We also thank Allan Sauvat for the help with the microscopy quantification. We thank Mathilde Soulez, Bernhard Lehnertz, Biljana Culjkovic, Brian Wilhelm, and Michaël Imbeault for insightful discussions. We are indebted to Kathie Béland and Elie Haddad from the CHU Sainte-Justine Research Center for providing the primary thymic samples.

Ethics

This project was approved by the CHU Sainte-Justine Research Ethics Board (protocol and biobank #2126). Written informed consent was obtained from parents of all 391 children involved.

This project was approved by the Comité de Déontologie de l’Expérimentation sur des Animaux of Université de Montréal (certificate 19-107).

Version history

  1. Preprint posted: July 14, 2023 (view preprint)
  2. Sent for peer review: August 20, 2023
  3. Preprint posted: November 3, 2023 (view preprint)
  4. Preprint posted: April 5, 2024 (view preprint)
  5. Version of Record published: April 18, 2024 (version 1)

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.91037. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2023, Larouche et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 760
    views
  • 87
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jean-David Larouche
  2. Céline M Laumont
  3. Assya Trofimov
  4. Krystel Vincent
  5. Leslie Hesnard
  6. Sylvie Brochu
  7. Caroline Côté
  8. Juliette F Humeau
  9. Éric Bonneil
  10. Joel Lanoix
  11. Chantal Durette
  12. Patrick Gendron
  13. Jean-Philippe Laverdure
  14. Ellen R Richie
  15. Sébastien Lemieux
  16. Pierre Thibault
  17. Claude Perreault
(2024)
Transposable elements regulate thymus development and function
eLife 12:RP91037.
https://doi.org/10.7554/eLife.91037.3

Share this article

https://doi.org/10.7554/eLife.91037

Categories and tags

Research organisms

Further reading

    1. Genetics and Genomics

    Mono-methylated histones control PARP-1 in chromatin and transcription

    Gbolahan Bamgbose, Guillaume Bordet ... Alexei Tulin
    Research Article

    PARP-1 is central to transcriptional regulation under both normal and stress conditions, with the governing mechanisms yet to be fully understood. Our biochemical and ChIP-seq-based analyses showed that PARP-1 binds specifically to active histone marks, particularly H4K20me1. We found that H4K20me1 plays a critical role in facilitating PARP-1 binding and the regulation of PARP-1-dependent loci during both development and heat shock stress. Here, we report that the sole H4K20 mono-methylase, pr-set7, and parp-1 Drosophila mutants undergo developmental arrest. RNA-seq analysis showed an absolute correlation between PR-SET7- and PARP-1-dependent loci expression, confirming co-regulation during developmental phases. PARP-1 and PR-SET7 are both essential for activating hsp70 and other heat shock genes during heat stress, with a notable increase of H4K20me1 at their gene body. Mutating pr-set7 disrupts monomethylation of H4K20 along heat shock loci and abolish PARP-1 binding there. These data strongly suggest that H4 monomethylation is a key triggering point in PARP-1 dependent processes in chromatin.

    1. Cancer Biology
    2. Genetics and Genomics

    Germline cis variant determines epigenetic regulation of the anti-cancer drug metabolism gene dihydropyrimidine dehydrogenase (DPYD)

    Ting Zhang, Alisa Ambrodji ... Steven M Offer
    Research Article

    Enhancers are critical for regulating tissue-specific gene expression, and genetic variants within enhancer regions have been suggested to contribute to various cancer-related processes, including therapeutic resistance. However, the precise mechanisms remain elusive. Using a well-defined drug-gene pair, we identified an enhancer region for dihydropyrimidine dehydrogenase (DPD, DPYD gene) expression that is relevant to the metabolism of the anti-cancer drug 5-fluorouracil (5-FU). Using reporter systems, CRISPR genome-edited cell models, and human liver specimens, we demonstrated in vitro and vivo that genotype status for the common germline variant (rs4294451; 27% global minor allele frequency) located within this novel enhancer controls DPYD transcription and alters resistance to 5-FU. The variant genotype increases recruitment of the transcription factor CEBPB to the enhancer and alters the level of direct interactions between the enhancer and DPYD promoter. Our data provide insight into the regulatory mechanisms controlling sensitivity and resistance to 5-FU.

    1. Computational and Systems Biology
    2. Genetics and Genomics

    Imputation of 3D genome structure by genetic–epigenetic interaction modeling in mice

    Lauren Kuffler, Daniel A Skelly ... Gregory W Carter
    Research Article

    Gene expression is known to be affected by interactions between local genetic variation and DNA accessibility, with the latter organized into three-dimensional chromatin structures. Analyses of these interactions have previously been limited, obscuring their regulatory context, and the extent to which they occur throughout the genome. Here, we undertake a genome-scale analysis of these interactions in a genetically diverse population to systematically identify global genetic–epigenetic interaction, and reveal constraints imposed by chromatin structure. We establish the extent and structure of genotype-by-epigenotype interaction using embryonic stem cells derived from Diversity Outbred mice. This mouse population segregates millions of variants from eight inbred founders, enabling precision genetic mapping with extensive genotypic and phenotypic diversity. With 176 samples profiled for genotype, gene expression, and open chromatin, we used regression modeling to infer genetic–epigenetic interactions on a genome-wide scale. Our results demonstrate that statistical interactions between genetic variants and chromatin accessibility are common throughout the genome. We found that these interactions occur within the local area of the affected gene, and that this locality corresponds to topologically associated domains (TADs). The likelihood of interaction was most strongly defined by the three-dimensional (3D) domain structure rather than linear DNA sequence. We show that stable 3D genome structure is an effective tool to guide searches for regulatory elements and, conversely, that regulatory elements in genetically diverse populations provide a means to infer 3D genome structure. We confirmed this finding with CTCF ChIP-seq that revealed strain-specific binding in the inbred founder mice. In stem cells, open chromatin participating in the most significant regression models demonstrated an enrichment for developmental genes and the TAD-forming CTCF-binding complex, providing an opportunity for statistical inference of shifting TAD boundaries operating during early development. These findings provide evidence that genetic and epigenetic factors operate within the context of 3D chromatin structure.