Abstract
15
Transposable elements (TEs) are genomic elements that are found in multiple copies in 16
mammalian genomes. TEs were previously thought to have little functional relevance but recent 17
studies have reported TE roles in multiple biological processes, particularly in embryonic 18
development. To investigate the expression dynamics of TEs during human early development, we 19
used long-read sequence data generated from in vitro differentiation of human pluripotent stem 20
cells (hPSCs) to endoderm, mesoderm , and ectoderm lineages to construct lineage -specific 21
transcriptome assemblies and accurately place TE sequences in their transcript context. Our 22
analysis revealed that specific TE types, such as LINEs and LTRs, exhibit distinct expression 23
patterns across different lineages. Notably, an expression outburst was observed in the ectoderm 24
lineage, with multiple TE types showing dynamic expression trajectories. Additionally, certain 25
LTRs, including HERVH and LTR7Y, were highly expressed in hPSCs and endodermal cells, but 26
these HERVH and LTR7 Y sequences originated from completely different transcripts. 27
Interestingly, TE-containing transcripts exhibit distinct levels of transcript stability and subcellular 28
localization across different lineages . Moreover, we showed a consistent trend of increased 29
chromatin association of TE-containing transcripts in germ lineage cells compared to hPSCs. This 30
study suggests that TEs contribute to human embryonic development through dynamic chromatin 31
interaction. 32
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
2
33
Key findings: 34
• Different loci of the same TEs are independently regulated in different cell states 35
• Ectoderm has the highest frequency of TE-containing transcripts 36
• The presence of TEs dynamically drives transcripts to different sub-cellular compartments 37
in different cell states 38
• hPSCs have the least stable TE transcripts with the weakest TE chromatin association , 39
highlighting loose hPSC chromatin and potential roles in cell differentiation 40
41
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
3
Main 42
TEs are genomic elements with multiple copies resulting from autonomous and non -autonomous 43
duplication in genomes 1,2. About half of the human genome consists of TEs of different types and 44
properties1,3,4. Because of their repetitive nature, TEs were previously considered to have little 45
functional importance5,6. However, recent studies have shown that TEs are functionally important 46
in diverse biological processes 2,4, including transcription 7–9, post -transcription expression 47
regulation10,11, transcript processing and stability 10,12,13, chromatin regulation 14,15, development16,17 48
and disease progression18,19. While TEs have been co-opted in normal developmental processes20,21, 49
TEs still pose a risk to genomic integrity 4,22. Therefore, TEs have both positive and negative 50
aspects6. Although TEs have been studied in pluripotent stem cells 21,23 and some terminally 51
differentiated somatic tissues 18,24, their expression patterns and their roles in post -implantation 52
human development and gastrulation have not been well explored. 53
The presence of multiple copies of TEs in the genome makes the investigation of their 54
functions difficult25,26. This is especially the case in TE expression quantification based on short -55
read RNA-seq in which a TE fragment often cannot be uniquely mapped to a genomic location27. 56
Thus, assembling transcripts from only short-reads is fraught with difficulty even in well-annotated 57
species such as human 28,29, and attempting to include assembly of TE -containing transcripts is 58
doubly difficult27. Previous efforts have considered global TE expression without consideration for 59
the specific TE loci 18,30,31, nor with their transcript context. Efforts are now being made to 60
investigate TE expression at specific loci 23,25, with an emphasis on understanding how TE 61
sequences are spliced into transcripts. Locus -level expression quantification has been reported to 62
be improved by transcript assembly 32. Although both short -read and long -read can be used for 63
transcript assembly 23,26, the quality of the transcript assembly based on long -read data is 64
superior23,33. Therefore, transcript assembly based on long -read data would substantially benefit 65
TE expression studies. 66
We have previously shown that TE sequences are richly expressed in hPSCs, including 67
inside coding and noncoding transcripts, and their presence is associated with changes in transcript 68
biophysical processes10. It is generally assumed that TEs are most active in the pluripotent stage 69
of development and decline in somatic tissues. However, this has not been explored, particularly 70
in the context of TE sequences inside transcripts. In this study, we investigated human early 71
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
4
development transcripts (hEDTs) assembled exclusively from long -read data to investigate TE 72
expression dynamics in human early development (hED). We firstly differentiated hPSCs in vitro 73
to the three germ layers; endoderm, mesoderm , and ectoderm, to mimic human gastrulation and 74
conducted long-read RNA sequencing for each cell type. By integrating deepCAGE and polyA 75
data, we showed that our transcripts are virtually complete from the 5’ to the 3’ ends of the 76
transcripts. Using the assemblies, w e showed that TE -containing transcripts were dynamically 77
regulated across different loci and cell states. We next investigated the expression dynamics, 78
subcellular localization, and transcript stability of the assembled hEDTs with a focus on the TE -79
containing transcripts and found that the stability and chromatin interaction of TE -containing 80
transcripts increases during the differentiation process. Overall, this study presents the expression 81
dynamics of TEs in early development and implicates TEs in cell state changes. 82
83
hPSC transcriptome assembled from long-read RNA-Seq data 84
We have previously described an hPSC-specific transcriptome based on a combined long-read and 85
short-read pipeline23. However, here we relied exclusively on deeply sequenced long-read data to 86
assemble transcripts. This strategy has the advantage that TEs will be placed into their exact 87
transcript context based only on long-reads, which will improve the accuracy of TE annotation in 88
transcripts. To this end, w e first generated PacBio long -read RNA-sequence data for hPSCs. 89
Consensus reads were generated from the PacBio subreads, and then processed to produce 90
noiseless reads using a standard PacBio read preparation pipeline (see Supplementary Methods) 91
(https://github.com/nf-core/isoseq). The noiseless reads were then mapped to the human genome 92
using Minimap234. Next, StringTie235,36 was used to assemble transcripts. After quality assessment, 93
33,375 hPSC transcripts, corresponding to 16,814 genes, were assembled. 94
Reference-based transcript assembly with long-read data are not completely error-free28,37,38. 95
Errors might arise because of the long -read region -biased error 39, reference -induced error 40, or 96
random RNA shearing during extraction. We, therefore, checked if our assembled transcripts were 97
full-length transcripts supported by long -read data, as opposed to annotation -aided assembly. 98
Indeed, the assembled transcripts were covered entirely by multiple consensus reads ( Fig. 1a), 99
suggesting that the transcripts were expressed in hPSCs. One advantage of transcript assembly 100
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
5
based on long-read data only is the ability to sequence full -length mostly intact transcripts26,27. In 101
support of this, most of the assembled transcripts were completely covered by at least one 102
consensus read ( Fig. 1b), suggesting that full -length transcripts were assembled. Interestingly, 103
almost all the consensus reads used for the assembly were used without trimming (Extended Data 104
Fig. 1 a). As a negative control, we generated a ‘random GTF set’ with the same number and 105
structure as the assembled transcripts (see Supplementary Methods). This random GTF was not 106
supported by our long-read data (Extended Data Fig. 1b). The pileup of consensus reads shows 107
that the reads map from the transcription start sites (TSSs) to the transcription end sites (TESs) of 108
the hPSC assembly (Fig. 1c). The pileup of short-read data equally validates end-to-end transcript 109
coverage for the hPSC assembly, but not for random GTF coordinates ( Extended Data Fig. 1c). 110
These data highlight the power of long-reads in complete transcript assembly. 111
To further assess the completeness of the assembled transcripts, we utilized published 112
hPSC deepCAGE41 and polyadenylated (polyA) data that marks the TSS and TES, respectively. 113
The long-read assembled TSSs were enriched for deepCAGE signal tags, while the TESs were 114
enriched for the polyA signal (Fig. 1d). These enrichments were not found in the random GTF set 115
(Extended Data Fig. 1 d). Similarly, ChIP-Seq data of hPSC POL2R, H3K27ac and H3K4me3 116
showed the expected promoter enrichment in the hPSC assembly, but not in the random GTF (Fig. 117
1e, Extended Data Figs. 1e, f). Additionally, H3k36me3 signal is enriched on the hPSC transcript 118
bodies (Extended Data Fig. 1g). The specific RNA-seq and histone modification data extensively 119
validate the completeness of the hPSC assembly. 120
We next checked if the splice junctions of the assembled transcripts were supported. The 121
evolutionary conservation scores across the junctions were higher relative to intron s. The highest 122
conservation was in the dinucleotide at the splice junctions for both the 5’ and 3’ ends ( Fig. 1f). 123
Indeed, the nucleotide frequencies showed the typical GT dinucleotide at the donor site, and AG 124
at the acceptor site. Analysis of short-read RNA-seq data with an anchor length of 10 bp, showed 125
that more than 99% of the splice junctions have at least one short-read junction support, and 85% 126
have 10 or more reads ( Fig. 1g, Extended Data Fig. 1h). Taken together, these results, from 127
multiple data, demonstrate the fidelity of the terminals and the splice junctions of the assembled 128
transcripts. 129
130
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
6
In vitro differentiation drives consistent and specific transcriptome trajectory in germ layers 131
To investigate transcript expression dynamics during human embryonic development, w e 132
performed in vitro differentiation of the hPSC to three germ layers using previously published 133
protocols42 (Fig. 2a). Lineage-specific marker genes were specifically up-regulated in each lineage 134
(Extended Data Fig. 2). Long-read sequencing was done for two technical replicates for each cell 135
state, and each sample had 20 -45 million PacBio raw subreads ( Extended Data Fig. 3a). The 136
consensus reads generated from each replicate ranged from 509 thousand to 1.2 million. Individual 137
assembly produced 31,245 endoderm, 37,161 mesoderm , and 56,063 ectoderm transcripts ( Fig. 138
2b). These corresponded to 14,989 endoderm, 17,520 mesoderm , and 24,766 ectoderm genes 139
(Extended Data Fig. 3b). Since the number of assembled transcripts can potentially be affected 140
by the sequencing depth 27, we asked if the higher transcript count in ectoderm cells is a genuine 141
phenomenon or a consequence of greater sequencing depth in these samples. Interestingly, 142
transcript assembly with individual ectoderm samples produced more transcripts even though the 143
sequencing depths were lower than the combined sequencing depths in other cell states (Extended 144
Data Figs. 3a-b), suggesting that transcript complexity is higher in ectoderm cells . To further 145
confirm this, we sampled 300,000 consensus reads from the replicates of each cell state. Transcript 146
assembly from the uniform read counts revealed that more transcripts were assembled in ectoderm 147
(Extended Data Fig s. 3c). As the impact of sequencing depth tends to be more pronounced in 148
lowly expressed transcripts27, we checked if the higher transcript count in ectoderm would persist 149
at a higher expression threshold. Indeed, that was the case up to 20 counts per million (CPM) 150
expression threshold (Extended Data Figs. 3d), supporting the increased complexity of ectoderm 151
transcripts. 152
Transcript quantification using the merged transcript reference showed that samples from 153
the same cell state tended to cluster together and separated from other samples (Fig. 2c, Extended 154
Data Fig. 3e). Importantly, the differentiated samples were all separated from hPSC samples, and 155
did not cluster with samples from other cell states, suggesting that different cell states had a specific 156
transcriptomic landscape. Indeed, previously generated long -read data of H1 and c11 cell lines 23 157
consistently clustered with hPSC data, demonstrating the robustness of our pipelines to technical 158
replicates and different hPSC cell lines. These results demonstrated the reliability of our in vitro 159
differentiation experiments and transcript assembly and quantification procedures. 160
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
7
We then proceeded to investigate the features of hEDT across different cell types. Ectoderm 161
transcripts tended to have longer spliced transcripts with fewer exons per transcript (Mann-162
Whitney U p value < 1 × 10-5) (Figs. 2d, e ). However other parameters, including un -spliced 163
transcript lengths, exon lengths and numbers of transcripts per gene were not substantially different 164
across different cell states ( Extended Data Fig. 3f). We next predicted coding and noncoding 165
potential using FEELnc 43, and transcript distributions based on protein coding potential showed 166
that ectoderm tended to have more noncoding transcripts and expressed more noncoding RNA 167
reads (Fig. 2f). The comparisons of the assemblies of each cell states to the GENCODE (v42) 168
transcript set 44 revealed that ectoderm indeed had a lower proportion of GENCODE -known 169
(matching) transcripts ( Extended Data Fig. 3g). Noncoding transcripts tended to be novel 170
transcripts in all cell types, with a higher proportion of novel transcripts in ectoderm ( Fig. 2g). 171
Interestingly, base-level conservation levels, nucleotide frequencies of the splice junctions and the 172
distributions of alternative splicing events were similar across each cell state (Extended Data Figs. 173
3h-j). 174
To investigate the heterogeneity of transcript expression at the single cell level, previously 175
published scRNA -seq data 45 containing hPSCs, and differentiated cell types was reanalyzed 176
against our hEDT assembly. UMAP clustering revealed that cells of similar states tended to cluster 177
together (Extended Data Fig. 4a). Meanwhile, hEDTs specifically expressed in each cell state and 178
the top 20 uniformly expressed transcripts were extracted from the long -read data (Fig. 2h). We 179
then checked the expression of the selected hEDTs and found that the expression patterns were 180
consistently reflected in scRNA -seq cell states ( Fig. 2i, Extended Data Fig s. 3b-f). Both long-181
read and scRNA -seq data largely consistently showed that the in vitro differentiation and 182
bioinformatics procedures reliably captured the expressed hEDT set. 183
184
Different TE types are dynamically spliced into coding and noncoding transcripts during in 185
vitro differentiation 186
Having established the reliability of the in vitro differentiation experiments and the assembled 187
transcripts, we then investigated TE splicing patterns across the cell states. We used nhmmer to 188
annotate TEs inside transcript sequences based on a previously described approach 2,46. The 189
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
8
proportions of TE-containing transcripts and RNA read in the assembled transcripts revealed that 190
ectoderm expressed more TE -containing transcripts ( Extended Data Fig. 5a). As previously 191
reported23, most noncoding transcripts tended to contain more TEs than the coding transcripts (Figs. 192
3a, b). This observation was true when we considered the number of TE-containing transcripts and 193
total TE -containing RNA count, highlighting the coding potential disruption ability of TEs. 194
Notably, in both coding and noncoding transcripts, ectoderm had more TEs. Intriguingly, most 195
TE-containing transcripts contained less than half TE -derived sequence, and only a small 196
proportion of the transcripts contained >50% TE-derived sequences (Figs. 3c). Again, there was a 197
stark difference between coding and noncoding transcripts. TE -containing noncoding transcripts 198
contained more TE -derived sequences, although most noncoding transcripts still contained less 199
than 50% TE-derived sequences. The distributions of TE coverage were more similar in coding 200
transcripts than in noncoding transcripts. Specifically, ectoderm TE -containing noncoding 201
transcripts tended to contain more intermediate values. However, hPSC and endoderm contained 202
the highest proportion of transcripts consisting entirely of TE -derived sequence ( Fig. 3c) even 203
though they had fewer overall TE-containing transcripts (Fig. 3a). 204
To investigate the specific TE splicing pattern further , we checked which TE types were 205
contained in hEDTs. We found that the most frequently spliced TE types included SINE, 206
retroposon, LINE, LTR, and DNA (Fig. 3d). It is important to note that these TE frequencies did 207
not reflect the genomic TE distribution (Fig. 3e), suggesting the TEs are not just randomly spliced 208
from genomic sequences . Across the major TE types, frequent TE splicing in ectoderm was 209
observed with almost half of the transcripts containing at least one SINE, compared to about a 210
third in hPSC transcripts (Fig 3d). Indeed, the higher frequency of SINE and retroposon splicing 211
is observed across cell states in both coding and noncoding transcripts ( Extended Data Figs. 5b, 212
c). On the contrary, less than 10% of the coding transcripts contained LTR-containing transcripts. 213
The coding bias of the TEs was evident as the proportion of coding transcripts varied across cell 214
states and TE types, ranging from 33% for LTR in ectoderm to 63% for SINE in endoderm 215
(Extended Data Fig. 5d). While most of the TE-free transcripts were either known GENCODE 216
transcripts (matching) or new isoforms of known GENCODE transcripts (variants), TE-containing 217
transcripts, especially the noncoding transcripts, contained more novel transcripts (that do not 218
overlap any known GENCODE transcript) with the highest novel proportion (~60%) for DNA TE-219
containing transcripts in ectoderm noncoding transcripts (Extended Data Fig. 5e). As TEs could 220
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
9
potentially create a minor transcript as an alternative isoform of a major functional transcript , we 221
investigated if TE-containing transcripts tended to contain fewer major transcripts (transcripts with 222
the highest expression among the isoforms of a single gene). For coding transcripts, the TE -free 223
sequences tended to have a higher proportion of major transcripts, compared to the TE-containing 224
transcripts (Extended Data Fig. 5f). Conversely, TE-containing noncoding transcripts tended to 225
contain a higher proportion of major transcripts, suggesting that many TE -containing coding 226
transcripts were the products of alternative splicing. The data demonstrated that our long -read 227
assembly uncovered multiple previously unannotated transcripts and the bias of the insertion of 228
various TE types into transcripts. 229
TEs might contribute to cell states or serve as functional modules for cell state 230
transitions47,48. As previously reported 49, noncoding transcripts had higher expression variability, 231
as quantified by the coefficient of variation, across cell types ( Fig. 3f). Further, TE -containing 232
coding and noncoding transcripts tended to have higher expression variability than TE -free 233
transcripts. Interestingly, the Tau coefficient which measures expression specificity showed that 234
both TE-containing coding and noncoding transcripts had higher cell-specific expression than TE-235
free transcripts, suggesting that TE expression is regulated in a cell type-specific manner. We next 236
checked the expression variability within cell states using scRNA -seq data and found that hPSC 237
transcripts had lower expression variability than differentiated cells ( Fig. 3g). Except for 238
mesoderm LTR-containing transcripts, noncoding transcripts tended to have higher expression 239
variability than the coding transcripts. Notably, TE presence led to higher expression variability in 240
coding transcripts but lower variability in noncoding transcripts. The differences in the expression 241
variability might be influenced by the expression detectability, especially for lowly expressed 242
transcripts. We checked the expression detectability and found that TE -containing transcripts 243
tended to have higher detectability (Fig. 3h). Interestingly, while the transcript detectability in TE-244
free noncoding transcripts was lower than that of the TE-free coding transcripts, the detectability 245
of the TE-containing coding and noncoding transcripts were more similar. These data showed that 246
the low expression variability in TE-containing noncoding transcripts was not just because of the 247
inability to detect due to the low transcript expression. Taken together, these results suggest that 248
the trajectory of TE expression is coordinated during differentiation, but cell -level coding and 249
noncoding transcript expression is dynamic. 250
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
10
251
Biased TE splicing patterns and frequencies of TE sequences in coding and noncoding 252
transcripts 253
Having established the dynamic expression of TE -containing transcript in different lineages, we 254
investigated the splicing patterns and insertion frequencies of different TEs in the four cell states. 255
Because of the differences in coding and noncoding transcripts, we investigated the splicing 256
patterns separately. Analysis of TE-splicing patterns in coding transcripts showed that TE 257
sequences are rare in coding sequence (CDS) regions ( Fig. 4a, Extended Data Fig. 6a), 258
presumably reflecting the evolutionary cost of a TE inserting into and disrupting a coding sequence. 259
Across TE types and cell states, the overrepresentation of TEs in the 5’ and 3’ untranslated regions 260
(UTRs) was high, compared to the CDS, with higher frequencies in the 3’ UTRs. However, 261
investigation of specific TE s showed differences between different regions and cell states. For 262
example, while HERVH and MSTB are both LTRs, their TE splicing patterns were different, with 263
HERVH sequences being rare in the 3’UTR, whilst MSTBs were enriched ( Fig. 4a). In addition, 264
there was a higher frequency of HERVH TE -derived sequences in endoderm cells. Interestingly, 265
the TE pattern of X24_LINE was very similar to that of MSTB, and also showed higher splicing 266
in the UTRs of ectoderm, revealing the heterogeneity of the splicing patterns of specific TEs. The 267
comparison of the MER3 and LTR6A frequencies further highlighted differences between these 268
TEs and different cell states (Extended Data Fig. 6a). 269
We next looked at splicing at the TSS and TES boundaries of transcripts, and interestingly, 270
while the TE splicing patterns across the TSS and TES were different, there were no substantial 271
differences across cell states (Fig. 4b). Across the TE types, TE proportion was low at the TSS. In 272
contrast, TEs were also reduced at the TES, although this reduction extended 5’ into the transcript 273
body (Extended Data Fig. 6b). Similarly, TE frequencies at both the donor and acceptor splice 274
junctions were low (Fig. 4c, Extended Data Fig. 6c). Although there are subtle differences across 275
cell lines and TE types, the general patterns were similar (Extended Data Fig. 6c). Surprisingly, 276
TE frequencies of the exons at the splice sites suggest that TE insertion suppression at the exon 277
starts is stronger than the suppression at the exon ends. This supports the idea that TEs are impaired 278
from splicing into TSSs, TESs, or splice junctions directly to prevent the disruption of important 279
motifs for transcription initiation and termination and transcript splicing. 280
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
11
We next studied the properties of TEs in noncoding transcripts and found that TE splicing 281
patterns in noncoding transcripts were dynamic between TE types, cell states , and even positions 282
within the transcript (Fig. 4d, Extended Data Fig. 6d). Specifically, while SINE and retroposon 283
TEs are enriched towards the 5’ and 3’ ends of the transcripts, LINE, LTR , and DNA TEs are 284
uniformly distributed. Additionally, the overall frequencies were different across cell states ( Fig. 285
4d, Extended Data Fig. 6d). Ectoderm noncoding transcripts were enriched for LINE and DNA 286
TEs, whilst LTRs were enriched in hPSCs, and were also high in ectoderm and endoderm ( Fig. 287
4d). The dynamics of TE expression was more obvious for transcripts containing LTRs. For 288
example, while HERVFH21 was concentrated in the center of noncoding transcripts, endoderm 289
LTR6A had higher enrichment around 5’ and 3’ ends, suggesting these transcripts are noncoding 290
remnants of intact ERVs. HERVH, conversely, was uniformly distributed across the noncoding 291
transcripts in hPSC and endoderm cells ( Extended Data Fig. 6d), highlighting the differences 292
among TEs. 293
TE frequencies around noncoding transcript TSSs were also dynamic. SINE, LTRs , and 294
retroposon were rare at the TSS ( Fig. 4e, Extended Data Fig. 6e). Conversely at the transcript 295
ends, SINEs and LTRs were enriched up to the TES, and then their enrichment dropped 296
substantially, suggesting that they marked transcript ends. TE sequences were absent at the TSSs 297
of noncoding transcripts, suggesting TE sequences rarely act as core promoters of TSSs. TE 298
frequencies at noncoding transcript splice junctions were also low ( Fig. 4e, Extended Data Fig. 299
6f), but not as low as in coding transcripts ( Extended Data Fig. 6f). Across the noncoding 300
transcript bodies, TSS, TES, and splice junctions ( Figs. 4d-f, Extended Data Figs. 6d-f), splice 301
pattern differences based on the TE type and cell states were obvious. 302
To figure out the cell type -specific TE activities, we explored their frequencies in coding 303
and noncoding transcripts. T he comparisons of the frequencies of TE -containing coding 304
transcripts identified 218 TEs that were significantly different in their frequencies upon 305
differentiation into endoderm, mesoderm, or ectoderm (Fig. 4g). The majority of the TEs that were 306
significantly different were enriched in ectoderm cells and were LINEs and LTRs and a few DNA 307
TEs. In contrast , some LTRs like LTR7, HERVH, HERVH48 , and MTSB2 were enriched in 308
endoderm coding transcripts. On the other hand, t he frequency comparison for noncoding 309
transcripts identified 331 specific TEs that were differentially enriched in at least one of the four 310
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
12
cell states (Fig. 4h). Similar to the coding transcripts, the majority of the differentially spliced TEs 311
in noncoding transcripts were LINEs and LTRs, and they were more frequent in the ectoderm. 312
Some LTRs such as LTR7B and HERVH were more frequent in hPSCs as previously reported50,51, 313
and here we observed many of them in endoderm cells, suggesting these are not wholly specific 314
features of hPSCs. Indeed, LTR7B/Y/H and HERVH were enriched in endoderm cells ( Fig. 4h), 315
supporting the idea that HERVH expression extends to cell types beyond hPSCs. HERVH 316
modulates 3-dimensional genome structure in hPSCs 52 and may also be performing this function 317
in somatic cells. Interestingly, there was a substantial overlap in the significantly enriched TE types 318
in any cell type for both coding and noncoding transcripts ( Fig. 4i). Breaking this down by TE 319
type, the overlap is high for LINEs, with just 21-specific for noncoding transcripts and 134 LINE 320
types were enriched in both coding and noncoding. For example, although the frequencies were 321
different, the overall abundance patterns for LINE L1MD_orf2 were similar for coding and 322
noncoding transcripts ( Fig. 4j). In contrast to LINEs, LTRs were different between coding and 323
noncoding transcripts, and only 49 LTR-types were common. This was exemplified by the LTRs 324
MSTB2 and LTR7Y which were specific to cell types (Fig. 4j). Overall, these data indicate that 325
the TE content in different lineages is highly dynamic, whilst the TE sequences in coding and 326
noncoding transcripts are more similar, except for LTRs, which are cell type and transcript -type 327
specific. 328
329
TE-containing coding and noncoding transcripts expressed from different genomic loci are 330
independently regulated 331
Above, we mainly consider TEs at the family or type level. However, a significant advantage of 332
long-read sequence data is that we can place TEs into their specific transcript context. Although 333
there is a substantial overlap in the types of TEs that are found in coding and noncoding transcripts, 334
the exact genomic loci or the TE-containing transcripts are likely to be different. Potentially there 335
are two scenarios, the first, is that the same TE type is contained in the same transcript, or that the 336
same TE type is expressed, but from different transcripts. To explore this, we first needed to 337
understand the relationship between TE type and expression level at the individual transcript level. 338
We used bulk short-read RNA-seq data for each lineage with the long-read based hEDT reference, 339
as the dynamic range of expression is higher for short -read data than long-read data27. In all four 340
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
13
cell states, the expression levels of TE-free coding transcripts were significantly higher (T-test p-341
value < 0.05) than those of TE -containing coding transcripts ( Extended Data Fig. 7a). TEs 342
inserted into the 5’ UTRs resulted in transcripts with the lowest expression levels while 3’ UTR 343
insertion led to comparatively higher expression levels, but the levels were not as high as those of 344
the TE-free transcripts. 345
We next investigated the aggregate expression patterns and found that the overall 346
expression patterns for SINE, LTR, LINE, retroposon , and DNA coding transcripts were not 347
substantially different across the four states (Fig. 5a). However, differential transcript expression 348
between hPSC and the three differentiated cell state s revealed that hPSC -to-endoderm 349
differentiation induced the least number of transcript changes while mesoderm differentiation 350
drove the largest coding transcript expression changes ( Extended Data Fig. 7b) . Differentiation 351
in all three lineages led to both upregulation and downregulation of TE -containing coding 352
transcripts, suggesting that transcript regulation is locus -specific and that TE expression is not 353
restricted to hPSCs. We checked the overlaps of the differentiation-induced differential transcripts 354
and found that many transcripts were regulated in a lineage-specific manner (Extended Data Fig. 355
7c). The aggregated expression of specific TE coding transcripts identified seven LTR TEs 356
(LTR7Y , LTR7C, LTR7B, LTR7, HERVH, HERVH48 and HERVFH21) that were significantly 357
downregulated in both ectoderm and mesoderm cells, and two LINE TEs (L1P4b_5end and 358
L1HS_5end) that were substantially upregulated in ectoderm cells (Fig. 5a, Extended Data Figs. 359
7d, e). This indicates that although the TE types are similar in all three cell types, they originate 360
from different transcripts. 361
We then examined the expression levels of the individual TE-containing coding transcripts. 362
Surprisingly, the majority of the HERVH-containing differentially expressed coding transcripts 363
were more highly expressed in endoderm cells but not in hPSCs (Fig. 5b). Different coding 364
transcripts containing L2 and L1HS_5end were activated in different cell states, with activation of 365
many L1HS_5end-containing coding transcripts in ectoderm cells. These data support the second 366
of our scenarios, at least for coding transcripts, that although similar TE types are expressed, they 367
are derived from different loci and are expressed in cell type-specific transcripts. 368
As noncoding transcripts are more likely to be expressed in a cell type -specific 369
manner27,53,54, we expected that noncoding TE-containing transcripts would also be expressed from 370
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
14
different transcript loci. However, in comparison to the coding transcripts ( Extended Data Fig. 371
7b), fewer differentially expressed noncoding transcripts were induced by the differentiation 372
process (Extended Data Fig. 7g), reflecting the overall reduced number of noncoding transcripts 373
in each cell type ( Fig. 2f ). Among the TE types, LTRs changed the most, with more LTR-374
containing noncoding transcripts downregulated in mesoderm cells. Also, many of the 375
differentially expressed transcripts were regulated in a cell type-specific manner, and only a few 376
differentially expressed transcripts were shared among cell states (Extended Data Fig. 7h). As for 377
coding transcripts, in all cell states, TE-containing transcripts had lower overall expressions levels 378
than TE-free transcripts (T-test p-value < 0.05) (Extended Data Fig. 7i). The comparison of the 379
aggregate expression for noncoding transcripts containing specific TEs showed that 126 TEs were 380
differentially expressed between the hPSC and at least one of the differentiated cell states ( Fig. 381
5c). These TEs included multiple LTRs that are activated in hPSC s and endoderm cells. 382
Interestingly, several SINE TEs were also differentially expressed. LTR6A and a few other LTR 383
TEs were also found to be enriched for endoderm cells. The analyses of the differentially expressed 384
LTR7-containing transcripts showed that many of the transcripts are downregulated in ectoderm 385
and mesoderm cells, whilst LTR6A-containing noncoding transcripts were up in endoderm ( Fig. 386
5d). Also, specific LINE-containing transcripts were specifically expressed in different lineages. 387
Interestingly, half of the TEs with significant expression differences also had significantly different 388
splicing frequencies (Extended Data Fig. 7j). These data support the model that similar TE types 389
are expressed from different transcripts in different cell types. 390
We next investigated the expression dynamics of TE -containing coding transcripts at the 391
single-cell level using scRNA -seq data from in vitro differentiated germ lineage cells 45. The 392
scRNA-seq data recaptured the higher expression of multiple TEs in ectoderm cells and 393
downregulation of some LTRs such as LTR6A, LTR7, LTR7C , and HERVH in ectoderm and 394
mesoderm cells (Extended Data Fig. 8a). The scRNA-seq data also revealed that the expression 395
differences were not obvious at the level of TE types but became distinct at individual transcripts 396
(Fig. 5e, Extended Data Fig. 8b). Also, single-cell analyses of the noncoding transcripts were 397
consistent with the observations from bulk data, and revealed upregulation of multiple TEs in 398
ectoderm cells (Extended Data Figs. 8a, b). Several LTR TEs were downregulated in ectoderm 399
and mesoderm, while LTR6A was enriched in endoderm cells (Fig. 5e, Extended Data Figs. 8c-400
d). As with bulk data , lineage -specific expression became more obvious when TE s were 401
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
15
considered at the transcript level (Fig. 5e, Extended Data Fig. 8d). These results demonstrated 402
that the observations in bulk RNA-seq were recaptured in single cells. 403
To further elaborate on the expression dynamics of TE-containing transcripts we extended 404
our analysis to scRNA-seq from in vivo gastrulation data55 (Fig. 5f, Extended Data Fig. 9a). As 405
expected, ectoderm cells were marked by more TE-containing coding and noncoding transcripts 406
(Extended Data Fig s. 9a, b) . Indeed, compared to other cell types, more ectoderm-expressed 407
transcripts were found in almost all of the TE types in both coding and noncoding transcripts 408
(Extended Data Fig s. 9c). Also, many transcripts that were identified in the in vitro data were 409
recaptured in the in vivo data ( Fig. 5d, Extended Data Fig. 9d), indicating that in vitro 410
observations reflected in vivo phenomena. Taken together, these data revealed that whilst similar 411
TE types are expressed in all cell types, they originate from multiple lineage -specific transcripts, 412
suggesting that TEs of the same type from multiple loci are independently regulated. 413
414
TE-type switching of transcripts during differentiation into somatic cells 415
As the same TE types are expressed from different transcripts, we next took advantage of 416
differentiation time course data to investigate how the expressions of TE -containing transcripts 417
vary during the endoderm differentiation process. As many LTR6A and HERVH -containing 418
transcripts are specifically expressed in the endoderm or hPSCs (Figs. 5b, d), we focused our 419
initial analyses on the transcripts containing these two TEs. The analyses of the time -course bulk 420
RNA-seq data showed that the activation of LTR6A noncoding transcripts happened at 72 hours 421
in the differentiation process (Figs. 6a, b). HERVH-containing noncoding transcripts with higher 422
endoderm expression were upregulated as early as 24 hours but did not reach maximum expression 423
until 72 hours. Conversely, HERVH -containing noncoding transcripts with high expression in 424
hPSCs were downregulated in endoderm at as early as 12 hours and were essentially undetectable 425
at 72 hours of differentiation to endoderm. The situation was similar for HERVH -containing 426
coding transcripts (Figs. 6a, Extended Data Fig. 10a). These results demonstrate how HERVH 427
and LTR6A can both be expressed in hPSCs and endoderm, but these TE sequences are contained 428
in lineage-specific transcripts. 429
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
16
Next, we analyzed scRNA-seq data of time-course hPSC-to-endoderm differentiation. Cell 430
clustering by identity showed that the cells formed three distinct clusters representing early (0h to 431
24h) middle (36h) and late (72h and 96h) stages of differentiation (Fig. 6c). In agreement with the 432
bulk RNA -seq analysis, endoderm -upregulated LTR6A , and HERVH -containing noncoding 433
transcripts had high expression only at 72 - and 96-hour of differentiation. In contrast, endoderm-434
downregulated noncoding transcripts were shut down at 12 hours of differentiation (Fig. 6d ). 435
These data showed that the timing for the activation and downregulation of TE transcripts varied 436
during differentiation. 437
To further investigate if the expression of the endoderm -induced upregulated transcripts 438
were maintained in terminally differentiated somatic cells, we first checked their expression in an 439
array of hepatocyte -related samples. Surprisingly, the expression of noncoding transcripts with 440
endoderm-induced upregulation was not sustained in hepatocyte-related terminally differentiated 441
somatic cells (Fig. 6e, Extended Data Fig. 10b, c). Similar observations were found for HERVH-442
containing endoderm-upregulated coding transcripts (Extended Data Fig. 10d). Importantly, the 443
expression of endoderm-specific HERVH-containing noncoding transcripts was not maintained in 444
endoderm-derived terminally differentiated somatic cells (Extended Data Fig. 11). Conversely to 445
HERVH, the expression of noncoding transcripts containing L1P1_orf2 and LTR that were 446
specifically expressed in the ectoderm were maintained in multiple ectoderm-derived somatic cells, 447
such as neurons (Extended Data Fig. 11). Similarly, the expression of many ectoderm-activated 448
coding transcripts containing L2, L1P1_orf2 and LTR remained active in ectoderm -derived 449
terminally differentiated somatic tissues. In summary, these results demonstrate how TE types can 450
remain active across cell types, but they are often expressed from different transcripts in distinct 451
tissues. 452
453
TE presence influences transcript subcellular localization 454
Previous studies have shown that TE-containing transcripts tend to localize in the nucleus 23,56. To 455
this end, we generated RNA -seq data for subcellular fractions of hPSCs, mesoderm, endoderm, 456
and ectoderm differentiated cells for the nucleus and cytoplasm. Western blots for proteins known 457
to localize to the nucleus and cytoplasm confirmed that the intended cellular fractions were purified 458
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
17
(Extended Data Fig. 12a). We confirmed our previous observation that both coding and 459
noncoding TE-containing transcripts preferentially localize to the nucleus in hPSCs 23 (Fig. 7a). 460
This pattern was the same for mesoderm and endoderm. Surprisingly, this was not the case for 461
ectoderm, where the opposite pattern was seen: TE-containing transcripts were instead more likely 462
to localize in the cytoplasm. Indeed, analyses of TE types revealed that they all tended to 463
preferentially localize to the cytoplasm in ectoderm cells, and to the nucleus in all other cell states 464
(Extended Data Fig. 12b). Greater details were revealed at the transcript level. For example, LTR-465
containing noncoding transcripts tended to have different cytoplasm/nucleus enrichment when 466
compared to other TE types in PSC and endoderm (Extended Data Fig. 12c). 467
Noncoding RNAs that localize to the nucleus are often recruited to chromatin 57,58. Hence, 468
we performed RNA-seq on the nucleoplasm and chromatin fractions to understand the sub-nuclear 469
localization of TE -containing transcripts ( Extended Data Fig. 12a). TE-containing hPSC 470
transcripts were reduced in the chromatin compartment, versus the nucleoplasm (Fig. 7b). On the 471
contrary, the TE containing transcripts of endoderm, mesoderm and ectoderm tended to be more 472
enriched in the chromatin fraction ( Fig. 7b, Extended Data Fig. 12d). Transcript-level analyses 473
of nucleoplasm/chromatin enrichment revealed heterogeneity across coding potentials, cell states 474
and TE types (Extended Data Fig. 12e). These data reveal that the sub -cellular and sub-nuclear 475
localization of TE -containing transcripts is cell type -specific, and retention of TE -containing 476
transcripts in the nucleus is not a general feature. 477
To further investigate the subcellular transcript localization, we compared the expression 478
levels of different transcripts to obtain transcripts enriched in specific sub -cellular fractions (see 479
Supplementary methods). Pairwise comparisons of nucleus -cytoplasm and chromatin -480
nucleoplasm expression patterns identified transcripts that were enriched in specific sub -cellular 481
fractions (Extended Data Fig. 13). The number of transcripts localized to the sub -cellular 482
compartments varied across cell types (Fig. 7c). Importantly, compared to the differentiated cells, 483
hPSCs had fewer transcripts differentially localized to the nucleoplasm and chromatin. The 484
distribution of coding and noncoding transcripts revealed both cell state and location -specific 485
differences. Specifically, while nucleoplasm -localized transcripts had the highest noncoding 486
proportion in hPSC, chromatin -localized transcripts had the highest noncoding proportion in 487
differentiated cells ( Extended Data Fig. 14a). Similarly, the distribution of TE -containing 488
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
18
transcripts varied across cell types, subcellular location, and protein -coding ability ( Extended 489
Data Fig. 14b). This effect was also transcript-specific, as few transcripts overlapped in the distinct 490
sub-cellular compartments between the different cell types (Extended Data Fig. 15). Indeed, using 491
a Jaccard similarity measure, transcripts that localized to the nucleus or cytoplasm were more 492
consistent in hPSC, endoderm and mesoderm but more divergent in ectoderm (Figs. 7d). 493
Interestingly, we found that the distribution of transcripts based on the number of cells in which 494
they were localized varied across subcellular structures ( Fig. 7e ). The most unique chromatin 495
localization was found in the mesoderm (Figs. 7d). Overall, nucleoplasm localization was the least 496
consistent across the cell states. To further explore the relationship between transcript subcellular 497
localization and cell state conversion, we checked the proportion of differentiation-induced DETs 498
among the transcripts localized to various subcellular structures across cell states and found that 499
subcellular localization significantly influenced the proportion of DETs ( Fig. 7 f), as the 500
distribution of DETs were significantly different in localized transcripts compared to overall 501
transcripts (Chi -square p-value < 0.05). For hPSCs, nucleoplasm -enriched transcripts had the 502
lowest proportion of DETs. For differentiated cells, however, the lowest proportion of DETs was 503
found in chromatin-localized transcripts. These data indicate that the localization of TE-containing 504
transcripts is cell type -specific. Overall, compared to the TE -free transcripts, we found that TE -505
containing transcripts were more cell type-specific (Fig. 7e), with hPSCs having fewer chromatin-506
enriched TE transcripts. 507
508
TE presence influences transcript stability 509
As subcellular localization can influence degradation, we wondered if this would impact transcript 510
half-life in a TE, sub-cellular localization, and cell state-dependent manner. We performed RNA-511
seq in differentiated cells treated with Actinomycin D to block transcription for 1 and 8 hours and 512
measured the transcript levels relative to the 0-hour time-point. Surprisingly, there were divergent 513
patterns of transcript half -life across TE type s, coding potentials, and cell states ( Fig. 7g). The 514
lowest stability was found in hPSCs and endoderm. Interestingly, the impact of TE presence was 515
also dynamic. While the TE-free transcripts were more stable in hPSC, ectoderm , and endoderm, 516
this was not the case in mesoderm (Fig. 7g, Extended Data Fig. 16a, b). TE-containing ectoderm 517
transcripts also contrasted with hPSCs and endoderm, as whilst TE -containing coding transcripts 518
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
19
were less stable, TE -containing noncoding transcripts were more stable than TE -free transcripts 519
(Fig. 7g, Extended Data Fig. 16a, b). These revealed that transcript stability varied across cell 520
states and transcript coding ability. 521
Finally, we checked if TE presence had any substantial influence on the stability of 522
transcripts with differential subcellular localization. The stability at 8 hours after transcription stop 523
showed differences across transcripts localized to different subcellular structures (Fig. 7h). While 524
nuclear-localized transcripts were more stable than cytoplasm-localized transcripts in the ectoderm, 525
the opposite was true in the other cell states. Conversely, chromatin -localized transcripts were 526
more stable than nucleoplasm -localized transcripts in hPSCs, but the nucleoplasm-localized 527
transcripts were more stable than chromatin -localized transcripts in differentiated cells. 528
Interestingly, some differences in stability can be observed as early as 1 hour after transcription 529
stop (Extended Data Fig. 16c). Moreover, the comparison of stabilities of TE-free transcripts to 530
those containing TEs identified multiple TEs associated with significant differences in stability 531
after 8 hours (Extended Data Fig. 16d). In hPSCs, TE presence mostly led to lower stability. In 532
differentiated samples, mosaic patterns were found. For example, TE -containing transcripts 533
localized to cytoplasm and nucleoplasm in mesoderm were more stable , along with cytoplasm-534
localized transcripts in ectoderm. These data suggest that TE presence might have different effects 535
on transcript localization and stability across different cell states. Overall, the data suggest that 536
stability and subcellular localization of TE-containing transcripts are dynamic across different cell 537
state conversions 538
539
References
600
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
22
1. Hutchins, A. P . & Pei, D. Transposable elements at the center of the crossroads between 601
embryogenesis, embryonic stem cells, reprogramming, and long non-coding RNAs. Sci Bull 602
(Beijing) 60, 1722–1733 (2015). 603
2. Ma, G., Babarinde, I. A., Zhou, X. & Hutchins, A. P . Transposable Elements in Pluripotent Stem 604
Cells and Human Disease. Front Genet (2022). 605
3. Biémont, C. A brief history of the status of transposable elements: from junk DNA to major 606
players in evoluXon. Gene3cs 186, 1085–93 (2010). 607
4. Bourque, G. et al. Ten things you should know about transposable elements. Genome Biol 19, 1–608
12 (2018). 609
5. Kim, Y .-J., Lee, J. & Han, K. Transposable Elements: No More ‘Junk DNA’. Genomics Inform 10, 226–610
33 (2012). 611
6. Naville, M. et al. Not so bad ader all: retroviruses and long terminal repeat retrotransposons as a 612
source of new genes in vertebrates. Clinical Microbiology and Infec3on 22, 312–323 (2016). 613
7. Wang, J., Huang, J. & Shi, G. Retrotransposons in pluripotent stem cells. Cell Regenera3on 9, 614
(2020). 615
8. Wang, J. et al. Primate-specific endogenous retrovirus-driven transcripXon defines naive-like stem 616
cells. Nature 516, (2014). 617
9. Elbarbary, R. A., Lucas, B. A. & Maquat, L. E. Retrotransposons as regulators of gene expression. 618
Science (1979) 351, (2016). 619
10. Babarinde, I. A. et al. Transposable element sequence fragments incorporated into coding and 620
noncoding transcripts modulate the transcriptome of human pluripotent stem cells. Nucleic Acids 621
Res 49, 9132–9153 (2021). 622
11. Fueyo, R., Judd, J., Feschoie, C. & Wysocka, J. Roles of transposable elements in the regulaXon of 623
mammalian transcripXon. Nat Rev Mol Cell Biol 23, 481–497 (2022). 624
12. Farré, D., Engel, P . & Angulo, A. Novel Role of 3’UTR-Embedded Alu Elements as Facilitators of 625
Processed Pseudogene Genesis and Host Gene Capture by Viral Genomes. PLoS One 11, 626
e0169196 (2016). 627
13. Gong, C. & Maquat, L. E. lncRNAs transacXvate STAU1-mediated mRNA decay by duplexing with 3’ 628
UTRs via Alu elements. Nature 470, 284–290 (2011). 629
14. Raviram, R. et al. Analysis of 3D genomic interacXons idenXfies candidate host genes that 630
transposable elements potenXally regulate. Genome Biol 19, (2018). 631
15. Sun, L., Fu, X., Ma, G. & Hutchins, A. P . ChromaXn and EpigeneXc Rearrangements in Embryonic 632
Stem Cell Fate TransiXons. Front Cell Dev Biol 9, 174 (2021). 633
16. Garcia-Perez, J. L., Widmann, T. J. & Adams, I. R. The impact of transposable elements on 634
mammalian development. Development 143, 4101–4114 (2016). 635
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
23
17. Hackei, J. A., Kobayashi, T., Dietmann, S. & Surani, M. A. AcXvaXon of Lineage Regulators and 636
Transposable Elements across a Pluripotent Spectrum. Stem Cell Reports 8, (2017). 637
18. Criscione, S. W., Zhang, Y ., Thompson, W., Sedivy, J. M. & Nerem, N. TranscripXonal landscape of 638
repeXXve elements in normal and cancer human cells. BMC Genomics 15, (2014). 639
19. LaRocca, T. J., Cavalier, A. N. & Wahl, D. RepeXXve elements as a transcriptomic marker of aging: 640
Evidence in mulXple datasets and models. Aging Cell 19, (2020). 641
20. Kruse, K. et al. Transposable elements drive reorganisaXon of 3D chromaXn during early 642
embryogenesis. bioRxiv 523712 (2019) doi:10.1101/523712. 643
21. Nellåker, C. et al. The genomic landscape shaped by selecXon on transposable elements across 18 644
mouse strains. Genome Biol 13, R45 (2012). 645
22. Gray, Y . H. M. It takes two transposons to tango:transposable-element-mediated chromosomal 646
rearrangements. Trends in Gene3cs 16, 461–468 (2000). 647
23. Babarinde, I. A. et al. Transposable element sequence fragments incorporated into coding and 648
noncoding transcripts modulate the transcriptome of human pluripotent stem cells. Nucleic Acids 649
Res 49, 9132–9153 (2021). 650
24. Chishima, T. et al. IdenXficaXon of Transposable Elements ContribuXng to Tissue-Specific 651
Expression of Long Non-Coding RNAs. Genes (Basel) 9, 23 (2018). 652
25. Berrens, R. V. et al. Locus-specific expression of transposable elements in single cells with CELLO-653
seq. Nat Biotechnol 40, (2022). 654
26. Babarinde, I. A., Li, Y . & Hutchins, A. P . ComputaXonal Methods for Mapping, Assembly and 655
QuanXficaXon for Coding and Non-coding Transcripts. Comput Struct Biotechnol J 17, (2019). 656
27. Babarinde, I. A. & Hutchins, A. P . The effects of sequencing depth on the assembly of coding and 657
noncoding transcripts in the human genome. BMC Genomics 23, 487 (2022). 658
28. Shumate, A., Wong, B., Pertea, G. & Pertea, M. Improved transcriptome assembly using a hybrid 659
of long and short reads with StringTie. PLoS Comput Biol 18, (2022). 660
29. Lagarde, J. et al. High-throughput annotaXon of full-length long noncoding RNAs with Capture 661
Long-Read Sequencing. Nat Genet 49, 1731 (2017). 662
30. He, J. et al. IdenXfying transposable element expression dynamics and heterogeneity during 663
development at the single-cell level with a processing pipeline scTE. Nat Commun 12, (2021). 664
31. Rodríguez-Quiroz, R. & Valdebenito-Maturana, B. SoloTE for improved analysis of transposable 665
elements in single-cell RNA-Seq data using locus-specific expression. Commun Biol 5, (2022). 666
32. Shao, W. & Wang, T. Transcript assembly improves expression quanXficaXon of transposable 667
elements in single-cell RNA-seq data. Genome Res 31, 88–100 (2021). 668
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
24
33. Bayega, A., Fahiminiya, S., Oikonomopoulos, S. & Ragoussis, J. Current and future methods for 669
mRNA analysis: A drive toward single molecule sequencing. in Methods in Molecular Biology vol. 670
1783 209–241 (Humana Press, New York, NY , 2018). 671
34. Li, H. Minimap2: pairwise alignment for nucleoXde sequences. Bioinforma3cs 34, 3094–3100 672
(2018). 673
35. Pertea, M. et al. StringTie enables improved reconstrucXon of a transcriptome from RNA-seq 674
reads. Nat Biotechnol 33, 290–295 (2015). 675
36. Pertea, M., Kim, D., Pertea, G. M., Leek, J. T. & Salzberg, S. L. Transcript-level expression analysis 676
of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc 11, 1650–67 (2016). 677
37. Chen, Y ., Zhang, Y ., Wang, A. Y ., Gao, M. & Chong, Z. Accurate long-read de novo assembly 678
evaluaXon with Inspector. Genome Biol 22, (2021). 679
38. Watson, M. & Warr, A. Errors in long-read assemblies can criXcally affect protein predicXon. 680
Nature Biotechnology 2019 37:2 37, 124–126 (2019). 681
39. Ono, Y ., Hamada, M. & Asai, K. PBSIM3: a simulator for all types of PacBio and ONT long reads. 682
NAR Genom Bioinform 4, (2022). 683
40. Nip, K. M. et al. Reference-free assembly of long-read transcriptome sequencing data with RNA-684
Bloom2. Nat Commun 14, (2023). 685
41. Abugessaisa, I. et al. FANTOM5 CAGE profiles of human and mouse reprocessed for GRCh38 and 686
GRCm38 genome assemblies. Sci Data 4, 170107 (2017). 687
42. Tsankov, A. M. et al. TranscripXon factor binding dynamics during human ES cell differenXaXon. 688
Nature 518, 344–9 (2015). 689
43. Wucher, V. et al. FEELnc: a tool for long non-coding RNA annotaXon and its applicaXon to the dog 690
transcriptome. Nucleic Acids Res 45, e57 (2017). 691
44. Frankish, A. et al. GENCODE: reference annotaXon for the human and mouse genomes in 2023. 692
Nucleic Acids Res 51, D942–D949 (2023). 693
45. Chu, L. F. et al. Single-cell RNA-seq reveals novel regulators of human embryonic stem cell 694
differenXaXon to definiXve endoderm. Genome Biol 17, (2016). 695
46. Wheeler, T. J. & Eddy, S. R. Nhmmer: DNA homology search with profile HMMs. Bioinforma3cs 29, 696
2487–2489 (2013). 697
47. Sakashita, A. et al. TranscripXon of MERVL retrotransposons is required for preimplantaXon 698
embryo development. Nat Genet 55, (2023). 699
48. Modzelewski, A. J. et al. A mouse-specific retrotransposon drives a conserved Cdk2ap1 isoform 700
essenXal for development. Cell 184, (2021). 701
49. Cabili, M. et al. IntegraXve annotaXon of human large intergenic noncoding RNAs reveals global 702
properXes and specific subclasses. Genes Dev 25, 1915–1927 (2011). 703
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
25
50. Santoni, F. A., Guerra, J. & Luban, J. HERV-H RNA is abundant in human embryonic stem cells and 704
a precise marker for pluripotency. Retrovirology 9, (2012). 705
51. Lu, X. et al. The retrovirus HERVH is a long noncoding RNA required for human embryonic stem 706
cell idenXty. Nat Struct Mol Biol 21, (2014). 707
52. Zhang, Y . et al. TranscripXonally acXve HERV-H retrotransposons demarcate topologically 708
associaXng domains in human pluripotent stem cells. Nat Genet 51, (2019). 709
53. Mamck, J. S. et al. Long non-coding RNAs: definiXons, funcXons, challenges and 710
recommendaXons. Nat Rev Mol Cell Biol 24, (2023). 711
54. Harrow, J. et al. GENCODE: The reference human genome annotaXon for the ENCODE project. 712
Genome Res 22, 1760–1774 (2012). 713
55. Tyser, R. C. V. et al. Single-cell transcriptomic characterizaXon of a gastrulaXng human embryo. 714
Nature 600, (2021). 715
56. Mas-Ponte, D. et al. LncATLAS database for subcellular localizaXon of long noncoding RNAs. Rna 716
23, 1080–1087 (2017). 717
57. Yin, Y . et al. U1 snRNP regulates chromaXn retenXon of noncoding RNAs. Nature 580, (2020). 718
58. Mishra, K. & Kanduri, C. Understanding Long Noncoding RNA and ChromaXn InteracXons: What 719
We Know So Far. Non-Coding RNA 5, 54 (2019). 720
59. Mustafin, R. N. The role of transposable elements in the differenXaXon of stem cells. Molecular 721
Gene3cs Microbiology and Virology (Russian version) 37, (2019). 722
60. Castro-Diaz, N. et al. EvoluXonally dynamic L1 regulaXon in embryonic stem cells. Genes Dev 28, 723
(2014). 724
61. Pereira Fernandes, D., Bitar, M., Jacobs, F. & Barry, G. Long Non-Coding RNAs in Neuronal Aging. 725
Noncoding RNA 4, 12 (2018). 726
62. Reilly, M. T., Faulkner, G. J., Dubnau, J., Ponomarev, I. & Gage, F. H. The role of transposable 727
elements in health and diseases of the central nervous system. J Neurosci 33, 17577–86 (2013). 728
63. Lapp, H. E. & Hunter, R. G. Early life exposures, neurodevelopmental disorders, and transposable 729
elements. Neurobiol Stress 11, 100174 (2019). 730
64. Kelley, D. & Rinn, J. Transposable elements reveal a stem cell-specific class of long noncoding 731
RNAs. Genome Biol 13, R107 (2012). 732
65. Kelley, D. R., Hendrickson, D. G., Tenen, D. & Rinn, J. L. Transposable elements modulate human 733
RNA abundance and splicing via specific RNA-protein interacXons. Genome Biol 15, 537 (2014). 734
735
736
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
-1 . 0 TSS TES 1.0Kb
0
5
10
15
deepCAGE
polyA
deepCAGE
-1. 0 TSS TES 1.0Kb
polyA
0.0
0.5
1.0
1.5
2.0
2.5
-1.0 TSS TES 1.0Kb
0.2
0.4
0.6
POLR2A
pPOLR2A
POLR2A
-1.0 TSS TES 1.0Kb
pPOLR2A
0.0
0.2
0.4
0.6
0.8
1.0
01 0 2 0 3 0
Ranked transcripts (,000)
0
50
100Coverage (%)
rep1
rep2
Merged
0 10 20 30
Ranked transcripts (,000)
0
50
100Coverage (%)
rep1
rep2
Merged
Completeness Transcription
Average CPM
Average CPM
Average CPM
-1.0 TSS TES 1.0Kb
20
40
rep1
rep2
Replicate 1
-1. 0 TSS TES 1.0Kb
Replicate 2
0
25
50
75
100
125
150
175
200
Expression
CPM CPMCPM
1 2 3 4 5 6 7 8 9 10
Minimum read coverage
Junction coverage (%)
99.06
98.04
96.76
95.33
93.76
92.1
90.42
88.64
86.87
85.11
99.15
98.24
97.16
95.89
94.5
93.01
91.47
89.87
88.26
86.61
99.64
99.27
98.85
98.38
97.85
97.26
96.6
95.9
95.22
94.46 rep1
rep2
Merged
0
50
100
Exon Intron
ExonIntron
a
b
c d e
f g
Transcript coverage
Longest read coverage
Fig. 1
Distance to junction (bp)
Acceptor splice junctions
Distance to junction (bp)
Donor splice junctions
G
A
C
C
G
T
A
C
T
A
GGT
C
T
GA
C
T
G
A
C
T
A
G
C
A
G
T
A
G
C
T
A
G
C
T
A
G
C
T
G
A
C
T
G
A
C
T
G
A
C
T
A
TCAG
T
C
A
G
T
0.0
1.0
2.0
Bits
10 0 10
0.0
2.5
5.0
phyloP score
10 0 10
2.5
5.0
0.0
1.0
2.0
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
26
Fig. 1: Quality assessment of the assembled hPSC long-read transcriptome. a, Consensus read 737
coverage of the assembled PSC transcripts using all consensus reads . b, The longest consensus 738
read coverage for the assembled PSC transcripts. c, Pileups for the consensus long reads for the 739
assembled hPSC transcripts across the exons from the TSS to the TES. Each transcript is scaled to 740
a uniform size. d, DeepCAGE (5’ end) and polyA tail (3’ end) pileups of the assembled transcripts, 741
from the TSS to the TES scaled to a uniform size and not including intronic regions . PolyA data 742
were from GSE138759 and GSE111134 while deepCAGE data were from GSE34448 and 743
GSE61264. e, Density pileup of POLR2A and phosphorylated POLR2A (pPOLR2A) for the 744
assembled hPSC transcripts. The data were from GSE242645. f, Position weight matrix nucleotide 745
frequencies and phyloP conservation scores of the 10-bp positions around the donor and acceptor 746
splice junctions of the assembled multi -exon hPSC transcripts g, Proportion of hPSC junctions 747
covered by indicated minimum number of reads. 748
749
750
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
ACTIVIN, WNT3A
TGF-
, WNT, PMB
inhibition
BMP4, VEGF, FGF2
PSC
Mesoderm
Ectoderm
Endoderm
0
50000
100000
33375
31245
37161
56063
98350
Transcript countPSC
Endoderm
Ectoderm
Mesoderm
Merged
PSC
Endoderm
Mesoderm
Ectoderm
Merged
26672
25392
29509
39266
67777
6701
5853
7651
16795
30568
Transcript counts
Matching Variant Novel
PSC
Endoderm
Mesoderm
Ectoderm
Merged
15928
14264
16711
20081
27485
9957
10296
11923
17648
36630
786
832
875
1537
3666
Coding
0 50 100
Transcripts (%)
PSC
Endoderm
Mesoderm
Ectoderm
Merged
1493
1191
1412
2282
3807
2916
2742
3299
6035
12425
2292
1920
2940
8478
14375
Noncoding
0 50 100
Transcripts (%)
PSC
Endoderm
Mesoderm
Ectoderm
993599
1107007
960744
1935175
80357
58053
64734
204251
RNA counts
NoncodingCoding
Cell state
Ectoderm
Endoderm
PSC
Mesoderm
ect_rep1
ect_rep2
mes_rep1
mes_rep2
mes_rep1
mes_rep2
end_rep1
end_rep2
ect_rep1
ect_rep2
H1_1
H1_2
H1_rep1
H1_rep2
c11_1
c11_2
Cell type
Cell marked
Cell type
PSC
Endoderm
Mesoderm
Ectoderm
Cell marked
PSC
Endoderm
Mesoderm
Ectoderm
All
0
0.2
0.4
0.6
0.8
1
Expression
Uniformly expressed PSC-specific Endoderm-specific Ectoderm-specific
0 10 20 300 2500 5000 7500 10000
Spliced transcript length Exons per transcript
PSC
Endoderm
Mesoderm
Ectoderm
Merged
PSC
Endoderm
Mesoderm
Ectoderm
Merged
a b c
d e
f g
h
i
Clustering by identity
Ectoderm
Endoderm
Fibroblast
PSC
Mesoderm
Trophoblast
20
10
0
10
10 0 10
UMAP_1
UMAP_2
2
3
4
hEDT_00083255
0.0
0.5
1.0
1.5
2.0
hEDT_00038452
0.0
0.5
1.0
1.5
hEDT_00065605
0.00.10.20.30.40.5
hEDT_00024447
012345
hEDT_00092327
0
1
2
3
4
hEDT_00003308
0.00
0.25
0.50
0.75
1.00
hEDT_00027673
0.00
0.25
0.50
0.75
hEDT_00083803
Fig. 2
end_rep1
end_rep2
H1_1
H1_2
PSC_rep1
PSC_rep2
c11_1
c11_2
ect_rep1
ect_rep2
mes_rep1
mes_rep2
end_rep1
end_rep2
H1_1
H1_2
PSC_rep1
PSC_rep2
c11_1
c11_2
Cell state
Cell state
0
1
0.2
0.4
0.6
0.8
R2
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
27
Fig. 2: Expression dynamics of the assembled human early development transcripts . a, 751
Experimental design for the in vitro differentiation of human hPSC to each of the three embryonic 752
germ layers. b, The numbers of transcripts assembled from each cell state and the merged transcript 753
assembly. c, The expression correlation heatmap of the samples made from long -read 754
quantification. H1_1, H1_2, c11_1 and c11_2 samples were from previously reported data 23. d, 755
Boxplots showing the spliced transcript lengths of the assembled transcripts. Boxplots show the 756
mean (red central bar), and second and third quartiles, and the whiskers show 1.5 times the 757
interquartile ranges, for this and all subsequent boxplots. Kruskal-Wallis p value = 2.37 × 10-28. 758
e, Boxplots showing the number of exons per assembled transcript. Kruskal-Wallis p value = 1.86 759
× 10-173. f, Stacked bar chart showing the proportions and numbers of coding and noncoding 760
transcripts in the hEDT assemblies. The upper plot shows the distribution of the assembled 761
transcript counts while the lower plot shows the distribution of the total RNA counts from the long-762
read quantifications. g, Stacked bar chart showing t he proportions and numbers of assembled 763
transcripts based on the similarity to GENCODE assembly. Matching transcript completely 764
matches to a known transcript (including all exons and splice junctions), while a variant transcript 765
partially overlaps a known GENCODE transcript. Novel transcript does not overlap any known 766
GENCODE transcript exon. h, Heatmap showing the uniformly expressed transcripts and 767
transcripts specific to each cell state based on the long-read quantification. For each transcript, the 768
expression was normalized by the highest expression level. i, Single-cell RNA-seq expression 769
showing uniformly expressed transcripts and marker transcripts specific to the indicated cell states. 770
771
772
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
TranscriptPSC
Endoderm
Mesoderm
Ectoderm
Merged
9050
8481
11418
17374
29875
17622
16911
18091
21892
37902
Coding
0 50 100
%
PSC
Endoderm
Mesoderm
Ectoderm
234816
251578
261646
629472
758783
855429
699098
1305703
RNA
TranscriptPSC
Endoderm
Mesoderm
Ectoderm
Merged
5612
5853
6606
15104
27019
1089
1082
1045
1691
3549
Noncoding
0 50 100
%
PSC
Endoderm
Mesoderm
Ectoderm
63783
42843
53813
177551
16574
15210
10921
26700
RNA
Coding
TE transcript TE contents
TE coverage (%)
Frequency of TE-containing transcripts
Noncoding
PSC
Endoderm
Mesoderm
Ectoderm
Merged
05 0 1 0 0
0
1
2
4
0
a b c
with TE no TE
d e
g
Coding
Noncoding
Proportion of TE transcript (%)
Fig. 3
with TE no TE
scRNA-Seq expression variation scRNA-Seq expression detectability
Coefficient of variation (CV) Cells with detectable expression (%)
Coding
Noncoding
Coding
Noncoding
Comparisons across different cell states
h
LINE
SINE
LTR
DNA
Retroposon
snRNA
0
10
20Coverage (%)
21.82
13.57
9.16
3.46
0.14
0.01
Genome coverage
LINE
SINE
DNA
LTR
Retroposon
noTE
Variation
0 1
CV
LINE
SINE
DNA
LTR
Retroposon
noTE
CodingNoncoding
Specificity
0 1
Tau
0
10
20
30
40
50 hPSC Endoderm
SINE
Retroposon
LINE
LTR
DNA
snRNA
Mesoderm
SINE
Retroposon
LINE
LTR
DNA
snRNA
0
10
20
30
40
50
Ectoderm
SINE
Retroposon
LINE
LTR
DNA
snRNA
Merged assembly
Ectoderm
0 3 6
Ectoderm
0 50 100
Mesoderm
0 3 6
Mesoderm
0 50 100
Endoderm
0 3 6
Endoderm
0 50 100
LINE
SINE
DNA
LTR
Retroposon
noTE
PSC
0 3 6
LINE
SINE
DNA
LTR
Retroposon
noTE
LINE
SINE
DNA
LTR
Retroposon
noTE
PSC
0 50 100
LINE
SINE
DNA
LTR
Retroposon
noTE
Coding
Noncoding
f
CodingNoncoding
CodingNoncoding
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
28
Fig. 3: Divergent TE expression dynamics across coding and noncoding transcripts . a-b, 773
Stacked bar charts showing the proportion and number of TE-containing transcripts in the coding 774
(a) and noncoding ( b) transcripts in the hEDT assemblies. The upper plots show the transcript 775
counts while the lower plots show the RNA counts from long-read based transcript quantification. 776
c, Histograms showing the distributions of TE coverages (Percent of sequence that consists of TE) 777
for both coding and noncoding transcripts. d, Stacked bar charts showing the proportion of TE -778
containing transcripts for coding and noncoding transcripts in the assembled hEDTs. e, Genomic 779
distributions of the top expressed TE types. f, Violin plots based on long -read quantification 780
showing the expression variability of TE-containing transcripts across different cell states. The left 781
charts show the expression variability measured by the coefficient of variation, while the right 782
charts show the expression specificity as measured by Tau. g, Single-cell RNA-seq expression 783
variability of different groups of transcripts within the cell states. h, Single -cell RNA -seq 784
expression detectability of different groups of transcripts within the cell states. T he red and blue 785
dotted lines represent the median values for the coding and noncoding transcripts without TEs. 786
787
788
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
0
5
10
Proportion of TE transcript (%)
All TEs
0
5
SINE
0.0
0.1
0.2
LTR:ERV1:HERVH
0.0
0.2
LTR:ERVL-MaLR:MSTB
TSS
0.0
0.5
1.0
LINE:L2:X24_LINE
TES
hPSC Mesoderm Endoderm Ectoderm
5' UTR CDS 3' UTR 5' UTR CDS 3' UTR
a b
TSS
0
1
2
LTR
TES Gene
Gene
-2 0 0 0 TSS 2000
0
25
All TEs
-2 0 0 0 TES 2000
25
50
All TEs
TSS TES
10
15
LINE
TSS TES
20
40
All TEs
TSS TES
15
20
25
Retroposon
TSS TES
20
30 SINE
TSS TES
7.5
10.0
12.5
LTR
TSS TES
5.0
7.5
DNA
hPSC Mesoderm
Endoderm Ectoderm
ed
-2000 TSS 2000
25
50
All TEs
-2000 TSS 2000
10
20
30
SINE
-2000 TSS 2000
20
40
LTR
-2000 TES 2000
25
50
-2000 TES 2000
10
20
30
-2000 TES 2000
15
20
25
Gene Gene
Gene
g h
17 130201
All TEs
02 1134
LINE
12 75 49
LTR
53 21 8
DNA
0
2
0
20
L1MD_orf2
0
1
0
10
MSTB2
0.0
0.2
0
5
LTR7Y
TE frequency (%)
Coding Noncoding
EctodermEndoderm
PSC Mesoderm
Coding Noncoding
c
Gene Gene
Gene
Exon
Intron
Exon Intron
f
-200 0 200
2.5
5.0
7.5
All TEs
-200 0 200
2.5
5.0
7.5
All TEs
Exon Intron
ExonIntron
-200 0 200
10
20
All TEs
-200 0 200
10
20
30
DonorAcceptor DonorAcceptor
Coding transcriptsNoncoding transcripts
Exon Intron
i
j
Significantly different TE types
Fig. 4
Ectoderm
Mesoderm
hPSC
Endoderm
L1P1_orf2
L1M8_5end
HAL1MEHERVH48
MSTB2
MLT1H2
MLT1A1
LTR33L2L1PB4_3end
L1PA4_3end
L1PA12_3end
L1P3_5end
L1ME3_3endL1ME3C_3end
L1MD_orf2
L1MB5_3end
L1M4_orf2
L1HS_5end
Tigger7Tigger15a
Ricksha
Charlie3
Cell state
TE type
Lineage count
Lineage
Cell state
PSC
Endoderm
Mesoderm
Ectoderm
Lineage
Ectoderm
Endoderm
Mesoderm
Mesoderm, Ectoderm
Lineage count
one
two
TE type
SINE
Retroposon
LINE
LTR
DNA
1.5
-1.5
0.0
Z score
Ectoderm
Mesoderm
hPSC
Endoderm
LTR7Y
LTR7BHERVH48HERVH
L1P1_orf2
L1M8_5end
HAL1ME
MIRb
LTR6A
MSTB2
MLT1H2
MLT1A1
LTR28
L3
L2L1PB4_3endL1PA4_3end
L1PA12_3end
L1P3_5endL1ME3_3end
L1MD_orf2
L1MB5_3end
L1M4_orf2
L1HS_5end
Tigger15a
Ricksha
Charlie3
Cell state
TE type
Lineage count
Lineage
1.5
0.0
-1.5
Z score
Proportion of TE transcript (%)
Proportion of TE transcript (%)
Proportion of TE transcript (%) Proportion of TE transcript (%)
Proportion of TE transcript (%)
Coding transcript TE frequencies Noncoding transcript TE frequencies
hPSC Mesoderm Endoderm Ectoderm
hPSC Mesoderm
Endoderm Ectoderm
hPSC Mesoderm
Endoderm Ectoderm
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
29
Fig. 4: TE splicing patterns and frequencies of TE -containing coding and noncoding 789
transcripts a , The splicing patterns of different groups of TEs in coding transcripts. Coding 790
transcripts were divided into 5’ UTR, CDS and 3’ UTR, and scaled to a uniform size. b, Proportion 791
of TEs 2 kbp upstream and downstream of coding transcript TSS and TES. c, Proportion of TEs 792
around 200 bp upstream and downstream of the donor and acceptor splice sites of coding 793
transcripts. d, The splicing patterns of different TE types in noncoding transcripts. Each transcript 794
is scaled to a uniform size. e, TE contents around the TSS and TES of noncoding transcripts. f, TE 795
contents around the donor and acceptor sites in the noncoding transcripts. For panels a-f, each 796
transcript was divided into 20 bins. The TE content for each bin was then computed. g, Heatmap 797
of TEs with significant difference in coding transcript frequencies between hPSC and any of the 798
three germ layers. h, Heatmap of TEs with significant frequency differences in the conceding 799
transcripts of the hPSC and any of the three germ layers. For panels g and h, the statistical 800
significance was defined as an adjusted p-value 2. For the heatmap, the z-801
score of the proportion of TE -containing coding and noncoding transcripts was computed. For 802
each TE, the lineage(s) with significant difference, and the number of lineage (s) with significant 803
differences are shown. i, Venn diagrams of the overlaps of the TE types with significant frequency 804
differences in coding and noncoding transcripts. j, TE frequencies of selected TEs in the coding 805
and noncoding transcripts. 806
807
808
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
Coding transcripts
LINE:L2
0 2
PSC_rp1
PSC_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
LTR:ERV1:HERVH
(n=139)
0 1
Expression (log2 TPM)
LTR:ERV1:LTR6A
(n=16)
0 1
LINE:L1:L1HS_5end
(n=301)
0 2
LINE (n=9137)
0 2
LTR (n=5809)
0 2
PSC_rp1
PSC_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
SINE (n=22391) LTR:ERV1:HERVH
Cell type
PSC Endoderm Mesoderm Ectoderm
LINE:L1:L1HS_5end
a b
TE transcript average expression
Cell type
PSC
Endoderm
Mesoderm
Ectoderm
hPSC
Endoderm
Mesoderm
Ectoderm
3
4
Bulk and scRNA-seq in vitro differentiation
Bulk and scRNA-seq in vivo
0.4
0.8
1.2
1.6
hPSC
Endoderm
Mesoderm
Ectoderm
1.5
2.0
2.5
hPSC
Endoderm
Mesoderm
Ectoderm
5
0
5
1
2
3
c d
f
Fig. 5
Coding transcripts
Noncoding transcripts Noncoding transcripts
PSC_rp1
PSC_rp2
ect_rp1
ect_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
PSC_rp1
PSC_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
PSC_rp1
PSC_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
Cell type
2
1
0
1
2
end_rp1
end_rp2
mes_rp1
mes_rp2
PSC_rp1
PSC_rp2
ect_rp1
ect_rp2
PSC_rp1
PSC_rp2
end_rp1
end_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
mes_rp1
mes_rp2
ect_rp1
ect_rp2
end_rp1
end_rp2
PSC_rp1
PSC_rp2
Cell type
PSC_rp1
PSC_rp2
end_rp1
end_rp2
ect_rp1
ect_rp2
mes_rp1
mes_rp2
Cell type
2
1
0
1
2
Z score
hESC.rp1
hESC.rp2
ect_rp1
ect_rp2
end.rp1
end.rp2
mes_rp1
mes_rp2
PRIMA41intLTR7YLTR46intHERVI
HERVHHERVFH21
HERV17
MSTB2
LTR6A
FAM
AluY
AluJr4
LTR83L1P3_5endL1HS_5end
Cell state
TE type
Lineage number
Lineage combination
Cell state
hPSC
Endoderm
Mesoderm
Ectoderm
TE_combination
Ectoderm
Endoderm
Mesoderm
Mesoderm, ectoderm
Endoderm, mesoderm
Endoderm, ectoderm
Lineage number
one
two
TE type
SINE
Retroposon
LINE
LTR
DNA
-1 0 1
Z score
Z score
2
1
0
1
2
Z score
LTR:ERV1:LTR7 LTR:ERV1:LTR6A
LINE:L2 LINE:L1:L1P3_5end
e
Ectoderm
Endoderm
Mesoderm
5
0
5
10 50 5
UMAP_1
UMAP_2
0.25
0.50
0.75
PSC
EctMes
EndTro
Fib
0.2
0.4
0.6
0.8PSC
EctMes
EndTro
Fib
0.4
0.8
1.2
1.6PSC
EctMes
EndTro
Fib
PSC: hPSC
End: Endoderm
Mes: Mesoderm
Ect: Ectoderm
Fib: Fibroblast
Tro: Trophoblast
5
0
hPSC
Endoderm
Mesoderm
Ectoderm
5
0
hEDT_00066295 (coding)
Bulk RNA-seq expression (log2 TPM)
5
0
Bulk scRNA-seq
hEDT_00024356 (noncoding)
hEDT_00094264 (noncoding)
0.5
1.0
1.5
2.0
hEDT_00002680 hEDT_00029253 hEDT_00062263
log2 (TPM)
log2 (TPM)
log2 (TPM)
Bulk scRNA-seq Bulk scRNA-seq Bulk scRNA-seq
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
30
Fig. 5: Dynamic expression patterns of TE -containing coding and noncoding transcripts 809
across different cell states. a, Boxplots of the aggregate expression levels of coding transcripts 810
containing specific TEs and TE types. b, Heatmap showing the expression changes of 811
differentiation-induced differentially expressed coding transcripts containing the stated TEs. Each 812
row of the heatmap represents a transcript containing the indicated type of TE (top label). c, 813
Heatmap of aggregate changes in expression of noncoding transcripts containing different TEs. 814
Each row represents the mean of all transcripts containing a TE-type that is significantly different 815
(adjusted p-value =2) between hPSC and any of the differentiated cell 816
states. For each TE, the lineage (s) with significant difference, and the number of lineage (s) with 817
significant differences are shown. d, Heatmap showing the change in expression levels of 818
differentiation-induced differentially expressed noncoding transcripts containing the indicated TEs. 819
Each row represents a TE-containing transcript that is differentially expressed between hPSC and 820
any of the three differentiated states. Bulk RNA -seq (left) and scRNA -seq (right) expression 821
patterns of selected TE-containing transcripts in in vitro (e) and in vivo (f) differentiation data. 822
823
824
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
a
e
c
EPS_cells
EPS1
hESC
end
EPS_HPLCs
hiHeps
EPS1_Heps
EPS2_Heps
F_PH1
F_PHHs
hFLCs
5
0
5 n=47
EPS_cells
EPS1
hESC
end
EPS_HPLCs
hiHeps
EPS1_Heps
EPS2_Heps
F_PH1
F_PHHs
hFLCs
5
0
5
Expression (log2)
n=251
LTR6A (Down in PSC)HERVH (Down in PSC)
Noncoding transcript expression (log2)
Differentiation time (hour)
5
0
5
n=256
0h
12h
24h
36h
72h
96h
HERVH (Up in PSC)
5
0
5 n=245
HERVH (Down in PSC)
5
0
5
n=44
LTR6A (Down in PSC)
PSC
Endoderm
Differentiated cells
Fig. 6
d
0h
12h
24h
36h 72h96h
6
3
0
3
6
50 5
UMAP_1
UMAP_2
0.001
0.002
0.003
0.004
Noncoding LTR6A
(down in PSC)
0.001
0.002
0.003
Noncoding HERVH
(down in PSC)
0.001
0.002
0.003
0.004
Noncoding HERVH
(up in PSC)
PSC to endoderm PSC to endoderm
0h
12h
24h
36h
72h
96h
LTR6A transcriptsHERVH transcripts
Time
Direction
Time
T0h
T12h
T24h
T36h
T72h
T96h
Direction
Up in hPSC
Down in hPSC
2
1
0
1
2
Z score
b PSC to endoderm
noncoding transcripts
5
0
5
hEDT_00015949
5
0
5
hEDT_00030238
0h
12h
24h
36h
72h
96h
5
0
5
hEDT_00066815
0.2
0.4
0.6
0.25
0.50
0.75
1.00
1.25
1
2
3
Bulk RNA-seq scRNA-seq
Bulk RNA-seq expression level (log2 TPM)
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
31
Fig. 6: Expression changes of TE-containing transcripts in endoderm lineage. a, Violin plots 825
of bulk RNA -seq expression levels of HERVH and LTR6A -containing noncoding transcripts 826
during differentiation to endoderm . The data were from GSE75748. b, Heatmap showing the 827
expression dynamics of HERVH and LTR6A noncoding transcripts differentially expressed in 828
endoderm lineage. c, UMAP showing the single-cell expression data for a hPSC to endoderm time 829
course. Hours of differentiation are shown. The three other UMAPs show example aggregate 830
scores for transcripts containing the indicated TE types that are down -regulated in hPSCs or up -831
regulated in hPSCs. d, Expression dynamics in bulk RNA -seq (left) and scRNA -seq (right) of 832
selected transcripts. e, Expression dynamics of endoderm -upregulated HERVH and LTR6A -833
containing noncoding transcripts in hepatocyte -related cells, using transcripts containing the 834
indicated TE types, that are defined as down-regulated in hPSCs. (EPS: Extended pluripotent stem 835
cells; EPS1: Stage 1 EPS; EPS2: Stage 2 EPS; hESC: Human ESC (H1); HPLC: Hepatic 836
progenitor -like cells; EPS_HPLCs: EPS -derived HPLCs; hiHEPs: iPSC -derived hepatocytes; 837
EPS1_Heps: EPS1 -derived hepatocytes; EPS2_Heps: EPS2 -derived hepatocytes; F_PH1: Fetal 838
primary hepatocytes; F_PHHs: Fresh primary human hepatocytes; hFLC: Human fetal liver cells). 839
840
841
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
RCIcytoplasm/nucleus RCInucleoplasm/chromatinCoding
Noncoding
a b
2
0
2
0
Noncoding
2
0
Coding
2
0
PSCEndodermMesodermEctoderm
SINE
Retroposon
LINE
LTR
DNA
noTEc
Transcript stability
Time after transcription stop (hour)
18 1 8
g
0
5000
3483
4475
5755
2130
Cytoplasm
0
5000
5116
5621
7192
1929
Nucleus
0
2500
1667
4156
2618
4315
Chromatin
0
2000
651
3452
1591
2707
Nucleoplasm
PSC
Endoderm
Mesoderm
Ectoderm
Transcript count
d
0 20 40 60 80 100
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
with TE
0 20 40 60 80 100
Transcript (%)
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
no TE
one
two
three
four
Fig. 7
Nucleus Cytoplasm
Chromatin Nucleoplasm
1.00
0.32
0.27
0.16
0.32
1.00
0.27
0.13
0.27
0.27
1.00
0.09
0.16
0.13
0.09
1.00
1.00
0.21
0.04
0.15
0.21
1.00
0.07
0.24
0.04
0.07
1.00
0.06
0.15
0.24
0.06
1.00
1.00
0.33
0.35
0.10
0.33
1.00
0.33
0.11
0.35
0.33
1.00
0.06
0.10
0.11
0.06
1.00
1.00
0.04
0.01
0.06
0.04
1.00
0.07
0.23
0.01
0.07
1.00
0.05
0.06
0.23
0.05
1.00
Cell state
PSC
Endoderm
Mesoderm
Ectoderm
0.0
0.5
0.1
0.2
0.3
0.4
5
0
5 PSC
Cytoplasm
Nucleus
Chromatin
Nucleoplasm
5
0
Ectoderm
2.5
0.0
2.5
Endoderm
Cytoplasm
Nucleus
Chromatin
Nucleoplasm
5
0
Mesoderm
Stability (8 hour relative to 0 hour)
Jaccard PSC
Coding Noncoding
EndodermMesodermEctoderm
LINE
SINE
DNA
LTR
Retroposon
noTE
LINE
SINE
DNA
LTR
Retroposon
noTE
LINE
SINE
DNA
LTR
Retroposon
noTE
5 0 5
LINE
SINE
DNA
LTR
Retroposon
noTE
5 0 5
Coding Noncoding
PSCEndodermMesodermEctoderm
LINE
SINE
DNA
LTR
Retroposon
noTE
LINE
SINE
DNA
LTR
Retroposon
noTE
LINE
SINE
DNA
LTR
Retroposon
noTE
5 0 5
LINE
SINE
DNA
LTR
Retroposon
noTE
5 0 5
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
Ectoderm
Localization in hPSC Localization in differentiated cell states
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
0
50
DET (%)
Endoderm
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
Mesoderm
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
0
50
Endoderm
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
Mesoderm
Nucleoplasm
Chromatin
Nucleus
Cytoplasm
All
Ectoderm
e
f
h
Downregulated
Upregulated
* *
** **
**
****
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint
32
Fig. 7: Subcellular localization and stability of TE -containing transcripts . a, 842
Cytoplasm/nucleus relative concentration index (RCI). RCI is the log2 -transformed relative 843
expression in two tissues. b, Nucleoplasm/chromatin RCI. c, Bar charts showing the number of 844
transcripts localized to different subcellular structures in different cell states using individual cell 845
state assemblies. d, Heatmap showing s ubcellular localization consistency of hEDT transcripts 846
localized to different structure across cell states . e, Stacked bars showing the distribution of the 847
number of cell states in which transcripts are localized to a subcellular structure. (* represents Chi-848
square significance at p value < 0.05). f, The proportion of localized transcripts that were 849
differentially expressed during differentiation process. The left charts show the localization in 850
hPSC while the right charts show the localization in differentiated cell states. The distributions for 851
all localizations were significantly different from the distribution of all hEDT set (Chi-square p 852
value < 0.05). g, Line plots showing the transcript stability of TE-containing transcripts, using time 853
from actinomycin D treated cells. The expression levels at time 1 and 8 hours, relative to the 854
expression at 0 hour are shown for each TE type across different cell states. h, Boxplots showing 855
the stability of transcripts that are localized to different subcellular structures across cell states (* 856
represents Mann-Whitney test significance at p value < 0.05). 857
.CC-BY-ND 4.0 International licenseavailable under a
was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint (whichthis version posted July 6, 2024. ; https://doi.org/10.1101/2024.07.03.602001doi: bioRxiv preprint