反转录转座子在大多数动物中非常普遍，占人类基因组的35%以上。然而，反转录转座子的流行、生物发生机制和功能在很大程度上仍然未知。在这里，我们开发了retroSeeker，这是一种新的计算软件，可以从基因组的成对比对中识别新的逆转录转座子，并解码它们的生物发生、表达、进化和潜在功能。我们发现大多数新的逆转录转座子表现出特定的L1内切酶切割基序，其中一些基序精确地位于插入位点上游的10个核苷酸。我们发现大量的候选功能基因可能通过反转录转位机制产生。重要的是，我们发现了与组蛋白基因、线粒体基因和穹窿rna相关的先前未被表征的反转录转座子类别。此外，我们阐明了逆转录转座子的组织特异性表达，并证明了它们在各种癌症类型中的普遍表达。我们还揭示了反转录转座子的复杂进化模式，并确定了许多物种特异性的反转录转座子事件。总之，我们的研究结果为发现新的反转录转座子类别和阐明它们在任何物种中的新特征建立了一个范例。
反转录转座子是大多数动物基因组中流行率最高的可移动元件，约占人类DNA的38% (Kazazian and Moran 2017;Wells and Feschotte 2020;Fueyo et al. 2022)。反转录转座子包括自主的长穿插元件(LINEs，由LINE-1表示)和非自主的短穿插元件(sin，由Alu和MIR表示)，它们需要line衍生的蛋白质进行反转录转座子(Moran et al. 1996;Cordaux and Batzer 2009)。逆转录转位的主要机制被称为靶引逆转录(TPRT) (Dewannieux et al. 2003)，通常被称为“复制和粘贴”机制(Gilbert et al. 2002)。简而言之，线衍生蛋白ORF1在目标位点切割第一条DNA链(通常由TTAAAA基序标记)，同时捕获转录RNA的聚(A)尾部(Gilbert等人，2002;dewanieux et al. 2003)。同时，另一种蛋白质ORF2利用断裂的基因组DNA作为引物(“复制”过程)，从反转录子转录物中产生互补DNA (cDNA) (Gilbert et al. 2002;dewanieux et al. 2003)。随后，目标DNA的第二条链被切割，反转录转座子的第二条链被合成(“粘贴”过程)(Gilbert et al. 2002;dewanieux et al. 2003)。这一系列事件导致反转录转座子侧翼的靶位点复制(TSDs) (Gilbert等人，2002;dewanieux et al. 2003)。因此，反转录转座子可以通过存在两个特征来区别于其他基因组重排:TSDs和反转录转座子侧面的3 ' -poly(A)信号。然而，由于逆转录转座子主要是随机整合到宿主基因组中(Fujiwara 2015)，对基因组插入的偏好及其对基因组塑造的影响在很大程度上仍然难以捉摸。
线衍生蛋白还促进蛋白质编码RNA (Maestre et al. 1995)和某些类型的非编码RNA (ncrna)的逆转录，如小核核RNA (snoRNA) (Weber 2006)、小核RNA (snRNA) (Doucet et al. 2015)和Y型RNA (Perreault et al. 2005)。逆转录转位产生新的功能性RNA拷贝，当亲本基因因突变失活时，这些RNA拷贝可以恢复亲本基因的作用，位于DKC1基因内含子内的小鼠ACA36 snoRNA就是一个例子(Weber 2006)。此外，与经典的LINE和SINE元件相比，各种反转录转位活性ncrna所表现出的序列多样性增加，使它们成为研究脊椎动物基因组进化的有价值的标记(Weber 2006)。然而，先前的研究主要依赖于同源性搜索和对一种反转录转位活性ncRNA的人工描述，因此缺乏全面和系统的探索。重要的是，其他种类的ncrna是否作为逆转录转座子的来源并通过LINE-1逆转录转座子机制进行逆转录尚不确定。
虽然大多数逆转录事件是有害的或中性的，但它们可以为创造新的功能基因提供原料，并已成为进化新颖性的主要来源。它们作为基因启动子，影响基因的选择性剪接，诱导现有基因的序列改变，甚至有助于基因的新生形成(Kaessmann et al. 2009)。与同样通过反转录产生的加工假基因相比，新功能基因对生理、形态、行为和生殖表型性状的进化具有深远的影响(Carelli et al. 2016)。它们被认为在推动哺乳动物大脑进化方面发挥着重要作用(Ferrari et al. 2021)。然而，由于缺乏适当的基因组尺度分析方法，导致产生新基因的具体逆转录转座事件以及由逆转录转座产生的新基因的实际数量尚不清楚。
在本研究中，我们利用来自不同物种的基因组两两比对来追踪逆转录转座子隐藏的“复制和粘贴”历史，并开发了一个名为retroSeeker的计算工具，促进了对逆转录转座子的全面探索。我们的研究揭示了在反转录转座子附近反复出现的特定基序。此外，我们的研究结果表明，反转录转座子通过两种主要机制促进新基因的出现:基因内区域整合和基因间区域整合。值得注意的是，我们发现了与组蛋白基因、线粒体基因组和穹窿rna相关的新型反转录转座子。有趣的是，许多反转录转座子在某些组织中表现出特定的表达模式，但在各种类型的癌症中表现出普遍的表达。此外，我们成功地确定了每个物种特有的许多反转录转座子，展示了复杂的进化轨迹。总之，我们的研究结果揭示了一大批新型反转录转座子的潜在生物发生机制、表达和进化。
1.1.1 一种新的计算方法设计用于发现反转录转座子的Nal方法
为了全面搜索反转录转座子，我们开发了一个名为retroSeeker的计算工具，该工具使用net文件，由lastz aligner和UCSC实用程序(见方法)生成的两个不同基因组的成对比对，来发现反转录转座子的复制和粘贴事件。在网络文件中，亲本基因产生了一个“填充”区域，复制粘贴逆转录在新的基因组位点内产生了一个“间隙”区域(Weber 2006)(图1A)。retroSeeker工具通过在两个不同物种之间进行成对基因组比对(图1A)来启动该过程，以生成相应的净文件。随后，它扫描隐藏的间隙和填充区域，这些区域源于在这些网络文件中识别的反转置事件。一旦确定了潜在的逆转录区域，我们采用动态规划算法对侧翼的TSD和poly(a)序列进行搜索和评分(方法)。
图1
发现逆转录转座子的新计算方法。retroSeeker算法和工作流程的示意图。左图:网络文件格式的特征。net格式文件来源于两个全基因组的成对比对(例如，鼠标查询到人类参考)。比对将产生两个结果，gap和fill，其中“gap”表示该区域存在于参考基因组中而不存在于查询基因组中，“fill”表示该区域同时存在于参考基因组和查询基因组中。人鼠比对网络格式文件显示了人和鼠基因组的合成部分比对。1级和2级是两个合成水平:1级对应人类反转录转座子宿主基因的小鼠同源物，2级对应人类反转录转座子的小鼠同源物。右图:2级填充区的特征和反转录转座子的测定。通过对人类序列的检查，可以识别聚A尾A(n)(红色区域)和TSD(绿色区域)，并精确定位两个TSD之间的反转录转座子插入点。B retroSeeker计算的逆转录转座子得分的接受者工作特征曲线;曲线下面积(AUC)值显示。C维恩图描述了retroSeeker和UCSC RetroGenes识别的独特和共享项目的数量。逆转录转座子的基因组浏览器可视化。顶部的反转录转座子的整个序列以不同的颜色显示;绿色表示两端的TSDs，蓝色表示基因体，红色表示poly(A)，其中poly(A)区域的核苷酸数量用括号表示，即A(n)。第一个轨道显示已知基因注释的信息。第二轨显示retroSeeker的识别结果。第三条轨道部分展示了retroSeeker的输入网络文件。E retroSeeker对不同的输入网络文件所花费的时间。显示了网络文件的基因组版本
为了评估retroSeeker的敏感性和特异性，我们进行了一系列的模拟研究。我们将一组line1随机插入模拟基因组多次，以产生包含这些新的反转录转座子的基因组(方法)。我们发现，推测的反转录转座子被召回具有高特异性(~ 1.00)和灵敏度(0.93，图1B)。值得注意的是，通过retroSeeker分析，在人类基因组中共发现了380,849个潜在的逆转录基因。这包含了UCSC基因组浏览器识别的约87%的逆转录基因，并将逆转录转座子目录扩展了14倍(图1C)。例如，retroSeeker成功地识别出插入PPP1R12B内含子内的FKBP8反转录转座子的存在，这与UCSC Genome Browser的发现一致(图S1A)。此外，我们还发现了几个众所周知的源自重复元件的反转录转座子品种，如LINE1(图1D)、Alu(图S1B)和SVA(图S1C)反转录转座子。此外，retroSeeker还展示了其揭示源自各种ncrna类型的逆转录基因的能力。值得注意的是，7SL RNA(图S1D)、snoRNA(图S1E)、snRNA(图S1F)和Y RNA(图S1G)都被成功鉴定出来，每个RNA都表现出完美的TSDs和典型多聚(A)区。
为了评估retroSeeker的运行速度，我们使用了来自不同物种的成对比对数据。值得注意的是，我们的大部分分析都在一分钟内成功完成，最快的完成时间为惊人的1.4 s(图1E、图S1H和图I)。这些结果表明，我们的retroSeeker方法不仅在发现各种反转录转座子方面具有高灵敏度、特异性和效率，而且还可以发现大量新的逆转录基因。
1.1.2 反转录转座子在不同进化过程中的独特特征连演化支
为了系统地研究逆转录转座子的特性，我们将retroSeeker应用于人类、小鼠和苍蝇(表S1)，并严格要求TSD长度≥7,poly(a)长度≥5(方法)。我们分别从人类、小鼠和果蝇中鉴定出139,181、121,074和7018个反转录转座子。值得注意的是，反转录转座子主要位于整个人类基因组的重复元件(51%)、内含子(20%)和基因间区域(10%)(图2A和表S2)。大多数与重复元件相关的反转录转座子被映射到灵长类特异性的Alu元件、LINE 1和内源性逆转录病毒(erv，图S2A)。在小鼠中，观察到类似的基因组区域分布模式(图2B和表S3)，并且发现大多数与重复元件相关的反转录转座子主要映射到啮齿动物特异性的B元件，包括B1 (15%)， B2(10%)和B4(3%)(图S2B)。有趣的是，在果蝇中，在蛋白质编码序列中发现了更高比例的反转录转座子(18%)(图2C和表S4)。此外，这些反转录转座子的主要部分与重复元件相关，特别是Gypsy(27%)、Jockey(7%)和Pao(4%)(图S2C)。
图2
反转录转座子的特性。已鉴定的逆转录转座子在人类(A)、小鼠(B)和果蝇(C)中按注释基因类型的分布。D至G箱形图显示了所鉴定的反转录转座子的TSD长度(D)、polyA长度(E)、反转录转座子长度(F)和得分(G)的值。每个框显示第一个四分位数、中位数和第三个四分位数。H, Circos图显示了人类中已鉴定的逆转录转座子的TSD长度、poly(A)长度和物种数量。从外圆到内圆的图例显示。I在人类反转录转座子5′-起始位点的上下游20nt序列中获得的序列基序。J饼图显示L1内切酶切割基序(TTAAAA)在人类反转录转座子中鉴定的百分比。K和L在人类(K)和小鼠(L)中逆转录转座子5′-起始位点的上下行20 nt序列中获得的序列基序。M靶位点引物逆转录(TPRT)事件中第二链切割特异性的示意图
为了揭示逆转录转座子的结构特征，我们检查了它们的单个成分的长度和分数。我们的分析显示，TSD长度从7到18个核苷酸(nt，图2D)， poly(a)长度从6到13 nt(图2E)，整个逆转录转座子长度从100到1500 nt(图2F)，总分从20到30(图2G)。然而，我们的研究没有发现TSD与每个反转录转座子的poly(A)长度之间的任何相关性(图2H，图S2D和E)。
为了研究反转录转座子中是否存在特定的序列元件，我们在5'起始位点的上游或下游的20个核苷酸中重新鉴定了基序。引人注目的是，我们发现大约20.2%的反转录转座子具有L1内切酶切割基序(图2I)。此外，44.4%的反转录转座子具有与L1内切酶切割基序相似的基序(图2J)。为了研究L1内切酶活性的位点特异性，我们计算了基序与插入的反转录转座子之间的距离。有趣的是，我们发现了一个包含尖锐的Alu元件(图2K)或B元件(图2L)的反转录转座子子集，位于L1内切酶切割基序下游10个核苷酸处，这表明在TPRT事件中，第二链切割发生在精确的位置(图2M)。值得注意的是，毗邻5'-TSD的区域往往含有富u元素(图S2F)或富a元素(图S2G)。
1.1.3 反转录转位产生新的基因并促进功能最终进化
据报道，在物种进化过程中，祖先基因的复制或转座因子的驯化可以产生新的基因(Carelli et al. 2016)。为了排除潜在的未表达的假基因，我们要求新基因在来自DNA元件百科全书(ENCODE)联盟的公开RNA-seq数据集中至少两个样本中表现出20多个独特reads的表达。我们发现，人类共有60228个反转录转座子显示出促进新表达基因产生的潜力(图3A和B)。此外，我们在小鼠中鉴定了25,063个反转录转座子(图S3A和C)，在果蝇中鉴定了3,385个(图S3B和D)，它们可能具有类似的意义。例如，一个反转录转座子事件导致了人类新基因RPEL1的出现(图3C)，而另一个反转录转座子事件导致了小鼠新基因Frg2f1和Uchl4的产生，这些基因都包括完美的TSDs和长聚(a)(图S3E和F)。此外，反转录转座子通过产生新的外显子间接地促进了新基因的产生。例如，一个反转录转座子被发现有助于人类ZNF729(图3D)和POU2F1(图S3G)外显子的产生，从而允许产生新的蛋白质翻译。
图3
通过反转录产生了新的基因。Circos图显示了人类中已识别的反转录转座子的最大log2(RPM)和样本数量。图中显示了逆转录转座子表达量最高的1000个。从外圆到内圆的图例显示。Log2(RPM):每百万次读取的Log2。B新发现的人类基因在标注基因类型中的分布。与RPEL1 (C)和ZNF729 (D)相关的反转录转座子的Genome Browser可视化。顶部的反转录转座子的整个序列以各种颜色显示;绿色表示两端的TSDs，蓝色表示基因体，红色表示poly(A)，其中poly(A)区域的核苷酸数量用括号表示，即A(n)。第一个轨道显示已知基因注释的信息。第二轨显示retroSeeker的识别结果。第三条轨道部分展示了retroSeeker的输入网络文件。E Circos图说明了与人类miRNA基因相关的反转录转座子的亲代关系。人类(F)和小鼠(G)中与miRNAs相关的反转录转座子的Genome Browser可视化。顶部的反转录转座子的整个序列以各种颜色显示;绿色表示两端的TSDs，蓝色表示基因体，红色表示poly(A)，其中poly(A)区域的核苷酸数量用括号表示，即A(n)。第一个轨道显示已知基因注释的信息。第二轨显示retroSeeker的识别结果。第三条轨道部分展示了retroSeeker的输入网络文件。“推定亲本基因> >逆转录基因> > miRNA”表示来自另一个位置的推定亲本基因(如KRT19)逆转录到新的位点，产生逆转录基因(如hsa-retrogene-5378)，最终形成miRNA(如MIR492)。人类亲本基因(H)和人类(I)、小鼠(J)和果蝇(K)中与反转录转座子相关的新基因富集的十大基因本体(GO)生物过程
有趣的是，我们检测到一组与源自蛋白质编码基因的mirna相关的反转录转座子(图3E)。例如，来自另一个位置的推定亲本基因KRT19 (chr17)逆转录到新的位点(chr12)，产生hsa-retrogene-5378，最终形成miRNA MIR492 (KRT19 > > hsa-retrogene-5378 > > MIR492，图3F)。为了探索这种miRNA产生机制在不同物种中的保守性，我们在小鼠中进行了相同的分析，我们检测了一组与来自小鼠蛋白质编码基因的miRNA相关的反转录转座子(图S3H)。具体来说，蛋白质编码基因Ftl1的反转录转座子产生了一个新的miRNA MIR692-2 (Ftl1 > > mmu-retrogene-249-7 > > MIR692-2，图3G)。
为了对亲本基因向新基因的功能转变进行初步研究，我们通过BLAST搜索和随后的基因本体(GO)富集分析，追踪了每个新基因的亲本-后代关系。我们的分析显示亲本基因主要参与核糖核蛋白复合物的生物发生和翻译(图3H)，而新基因主要与人类(图3I)、小鼠(图3J)和果蝇(图3K)的形态发生和神经元相关的生物过程相关。总的来说，这些发现强烈表明，逆转录在促进新基因的诞生中起着关键作用，从而加速了功能多样化的进化之旅。
1.1.4 来自多种基因的新型反转录转座子
为了进一步发现反转录转座子的新类别，我们重点研究了以前未被表征的候选反转录转座子的基因类型。值得注意的是，我们检测到一组组蛋白基因通过反转录转位产生新的拷贝。在人类(图4A)和小鼠(图S4A)中，这些拷贝主要来自H3，而在果蝇中，它们主要来自H2B(图S4B)。例如，在人类中，H5通过逆转录产生了一个新的H3-5拷贝，并且观察到它比亲本拷贝更长(图4B)。为了进一步表征来自组蛋白基因的反转录转座子，我们对其长度进行了统计分析，发现它们的长度主要在1000 ~ 3000 nt之间(图4C)。
图4
新型的反转录转座子。人类组蛋白基因反转录转座子的系统发育树。B逆转录转座子的基因组浏览器可视化。顶部的反转录转座子的整个序列以不同的颜色显示;绿色表示两端的TSDs，蓝色表示基因体，红色表示poly(A)，其中poly(A)区域的核苷酸数量用括号表示，即A(n)。第一个轨道显示已知基因注释的信息。第二轨显示retroSeeker的识别结果。第三条轨道部分展示了retroSeeker的输入网络文件。C累积曲线显示来自人类组蛋白基因的反转录转座子的长度。D Circos图说明了与人类线粒体基因相关的反转录转座子的亲代关系。与线粒体基因相关的逆转录转座子，包括rRNA (E)、mRNA (F)和mRNA簇(G)，以及与vault RNA基因(H)相关的逆转录转座子的Genome Browser可视化。顶部逆转录转座子的整个序列以各种颜色显示;绿色表示两端的TSDs，蓝色表示基因体，红色表示poly(A)，其中poly(A)区域的核苷酸数量用括号表示，即A(n)。第一个轨道显示已知基因注释的信息。第二轨显示retroSeeker的识别结果。第三条轨道部分展示了retroSeeker的输入网络文件
我们的研究还发现，线粒体rna可能通过反转录转座整合到核基因组中(图4D和图S4C)。例如，在人类的4号染色体中发现了一个假定的线粒体rRNA反转录转位(图4E)。此外，一个预测的线粒体mRNA反转录转位嵌入在11号染色体的基因间区域(图4F)，一个包含三个亲本基因的线粒体mRNA反转录转位嵌入在X染色体的基因间区域(图4G)，具有标准的多聚(a)区和TSD序列。
最后，我们发现非编码库rna是创造一类全新的反转录转座子的创新起源(图4H)。例如，我们在人类中发现了具有完美TSDs和poly(a)区域的拱顶rna的反转录转座子(图4H)。总之，这些结果表明，组蛋白基因、线粒体rna和其他ncRNA基因可能通过共享的l1介导的逆转录机制，作为产生新型逆转录转座子的新来源。
1.1.5 trans图集criptio最活跃的反转录转座子
为了全面研究逆转录转座子的组织特异性模式，我们使用ENCODE联盟提供的公开RNA-seq数据集进一步分析了它们在19种正常组织中的表达谱(图5A)。我们使用tau分数算法来计算逆转录转座子的组织特异性，发现许多逆转录转座子在特定的人类组织中表现出不同的表达模式(图5B和C)。总的来说，我们在人类中鉴定了80,941个tau分数高于0.8的组织特异性逆转录转座子，主要位于人类基因组的重复元件、内含子和基因间区域(图5D)。例如，hsa-逆转录转座子-25779在甲状腺中的表达水平明显高于其他组织(图5E)。为了探索逆转录转座子的肿瘤特征，我们使用从癌症细胞系百科全书(Cancer Cell Line Encyclopedia, CCLE)中获得的16种癌症的RNA测序数据分析了逆转录转座子的表达。有趣的是，我们发现反转录转座子的一个子集在各种类型的癌症中表现出普遍的表达，没有任何癌症类型特异性模式(图S5A)。
图5
组织特异性反转录转座子图谱。三种类型的表达数据集集合的示意图。B累积曲线显示在人类中发现的所有反转录转座子的组织特异性Tau分数。C热图显示逆转录转座子在各种人体组织中的表达谱，使用ENCODE的总RNA-seq数据。在细胞中的表达水平被分类到相应的组织中。log2RPM:每百万读取的log2。D .已鉴定的组织特异性反转录转座子在注释基因类型中的分布。具有代表性的组织特异性反转录转座子hsa-retrotransposon-25779的E表达值。箱形图显示不同组织中的log2RPM (RPM，每百万读取数)值。每个框显示第一个四分位数、中位数和第三个四分位数
为了研究表达的反转录转座子在亚细胞RNA中的分布，我们分析了亚细胞表达数据，观察到大多数反转录转座子表现出核特异性定位模式，在核质和核核中的表达水平高于在细胞质中的表达水平(图S5B)。
为了探索逆转录转座子在其他代表性物种中的组织特异性，我们进一步分析了小鼠和苍蝇的表达数据。我们在小鼠(图S5C)和果蝇(图S5D)中检测到一个组织特异性逆转录转座子子集。值得注意的是，我们发现大量反转录转座子在小鼠中脑中高度表达(图S5C)。总之，我们的研究结果表明，大量的反转录转座子在某些组织中表现出特定的表达模式，而在癌症中表现出普遍的表达。
1.1.6 多物种反转录转座子的复杂进化模式
为了研究来自不同物种的反转录转座子的进化关系，我们对多个物种的反转录转座子得分进行了主成分分析(PCA)。值得注意的是，我们观察到逆转录转座子可以被分为五个分支，与已知物种的系统发育完全一致(图6A)。为了探索不同类型逆转录转座子的代表性进化特征，我们首先对Alu亚家族的组成进行了统计分析(图S2A)。我们发现AluS亚家族占据了最大的比例(60%，图S6A)，它被认为是Alu元素中第二古老的亚家族。具体来说，通过反转录转座子扩增，AluY(4400)、AluSx(2161)和AluSx1(2086)亚家族成为最成功的反转录转座子(图S6B)。在这里，我们展示了几个在Alu亚家族中表现出完美的tsd和典型多聚(A)区域的例子(图S6C, D和E)。
图6
反转录转座子的复杂进化模式。多物种反转录转座子得分的主成分分析(PCA)。B，物种特异性反转录转座子鉴定示意图。C人类反转录转座子的简化系统发育树。内部分支和根，指示物种的同源反转录转座子家族数量。树尖，每个物种的反转录转座子数。维恩图描述了在灵长类动物和人类中被识别出的独特和共有物品的总数。E小鼠反转录转座子的简化系统发育树。内部分支和根，指示物种的同源反转录转座子家族数量。树尖，每个物种的反转录转座子数。维恩图，描述了在啮齿类动物和小鼠中识别出的唯一和共享项目的总数。G果蝇反转录转座子的简化系统发育树。内部分支和根，指示物种的同源反转录转座子家族数量。树尖，每个物种的反转录转座子数。维恩图描述了在果蝇和苍蝇中被识别的独特和共有项目的总数
为了系统地探索物种特异性反转录转座子并进一步研究其进化模式，我们将候选反转录转座子的序列映射到其他物种的基因组中，要求序列相似性≥0.95(图6B)。我们发现，在大约25万年前灵长类动物分化之后，发生了相当大比例的反转录转位事件(7355至57,905)(图6C，左图)。
为了进一步探索反转录转座子在物种进化中的基因组分布，我们将反转录转座子位点分别定位到不同物种中已知的基因注释上。有趣的是，我们发现在整个物种进化过程中，重复元件的比例增加到约50%，而编码序列(CDSs)的比例减少到约5%(图6C，右图)。我们发现，人类中超过一半的逆转录转座子是灵长类特异性的(62.5%)，但只有一小部分是人类特异性的(1.98%，图6D)。有趣的是，我们还检测到可能起源于2亿多年前的反转录转座子(Myr)(图6C)。
为了探索其他模式生物的进化模式，我们将分析扩展到小鼠和苍蝇。我们观察到，大多数逆转录转位发生在小鼠分化后(53,114至121,073)(图6E，左图)。此外，重复元件的比例增加到约50%，而假基因的比例下降到约5%(图6E，右图)。我们的调查显示，62.5%的反转录转座子是啮齿类动物独有的，其中47.7%是在小鼠中特异性鉴定的(图6F)。有趣的是，几乎所有在果蝇中发现的反转录转座子都属于果蝇属(96.8%)(图6G和H)。总之，我们的研究结果表明，反转录转座子组成的动态重排增加了物种基因组的复杂性和多样性。
摘要。
1 介绍
2 讨论
3 方法
数据和材料的可用性
参考文献。
致谢。
作者信息
道德声明
# # # # #
通过将逆转录转座子事件的生物学模型(复制和粘贴)映射到比较基因组学数据(空白和填充)中，retroSeeker能够准确地识别任何物种基因组中的新逆转录转座子。我们发现，通过保守的识别和插入机制，分布在编码/非编码基因和基因间区域的反转录转座子多样性增加(图S2A, B和C)。我们观察到，在人类和小鼠中鉴定的大部分反转录转座子在反转录转座子的5 '端周围有几个共识序列元件(图2K, L, S2F和G)。重要的是，共识序列元件之一经常出现在反转录转座子插入位点上游10个核苷酸处，形成位点特异性TTAAAAN(10)基序(图2M)。有趣的是，我们没有在果蝇中检测到这些基序，这表明这种机制可能是在哺乳动物分化后出现的，并促进了逆转录转座子在高等生物中的扩增。虽然在家蚕R2逆转录转座子中也报道了位点特异性基序，但它们在体内仅靶向28S rRNA基因(Wilkinson et al. 2023)。与具有严格限制的R2反转录转座子相比，基于ORF2的反转录转座子在整个基因组中靶向更广泛的基因类型(图S2A和B)。这一发现表明ORF2蛋白可能是体内基因编辑的潜在工具。
了解新基因的起源对于解释新表型和生物多样性的起源和进化的遗传基础至关重要。与以往的研究集中于有限的mrna衍生的基因重复相反，我们系统地探索了与反转录转位相关的所有基因。我们的工作表明，通过反转录转位形成新基因有两种主要途径(图3B)。第一种途径涉及亲本基因的反向转位复制，随后发生突变或新功能化过程(图3C)。第二种途径涉及将这些逆转录拷贝插入到先前存在的基因中，从而产生将逆转录转座子与基因序列结合在一起的嵌合RNA产物(图3D)。与亲本基因主要参与核糖核蛋白复合物的生物发生和翻译(图3H)相比，我们观察到，通过反转录转位产生的新基因主要参与各种细胞发育和分化过程(图3I)，表明其功能谱的扩展。这一观察结果与最近强调孤儿逆转录基因如何在功能上取代亲本基因的工作是一致的(Carelli et al. 2016)。例如，人类通过反转录转位形成的新基因POU2F1(图S3G)是高等真核生物的转录调控因子，参与调控发育、分化、应激反应等过程(Hamashima et al. 2023)。
有人提出，逆转录转位活性ncrna的数量和多态性可能与活性Alu或L1元件相当(Weber 2006)。然而，由于缺乏ncRNA逆转录转座子的完整列表，这个问题在很大程度上是未知的。因此，我们首次对RNA物种的整个范围进行了全面的研究。重要的是，我们观察到vault rna在人类中具有反转录转位能力(图4H)，阐明了哺乳动物中ncRNA基因产生功能拷贝的保守机制。有趣的是，我们检测到线粒体(MT) rna可能嵌入核(NU)基因组(图4D和C)，产生核线粒体片段(numt)。numt存在于整个人类核基因组中，一些numt与疾病有关(Xue et al. 2023)。与之前提出的通过双链DNA断裂修复整合NUMT的机制(Wei et al. 2022)相反，我们的工作揭示了NUMT生成的一种新的可选途径。然而，mtDNA/mtRNA片段如何离开线粒体并转运到细胞核尚不清楚(Xue et al. 2023);因此，未来的研究需要进行更有效的验证。
反转录转座子曾被认为是由于缺乏转录元件而不表达的“基因组垃圾”。然而，最近的研究强调了它们在组织分化中的重要调节作用(Roller et al. 2021;Nam et al. 2023)。在这项研究中，我们系统地鉴定了一组具有组织特异性表达模式的反转录转座子，揭示了它们在基因表达调控中的作用。我们检测到相当比例的反转录转座子嵌入内含子或基因间区域(图5D)。因此，我们假设这些反转录转座子可能通过产生新的启动子、增强子或其他调控元件(如增强子rna)来塑造组织特异性基因调控模式，从而促进新的外显子或组织特异性。有趣的是，我们还观察到大量逆转录转座子在大脑中显著高表达(图S5C)，这与最近强调逆转录转座子是哺乳动物大脑进化的重要驱动因素的研究结果一致(Ferrari et al. 2021)。此外，反转录转座子活性的失调也可能导致神经系统疾病。因此，我们假设在其他组织中具有特异性的反转录转座子也可能在其他疾病中发挥重要作用。
总之，我们的研究提供了有价值的软件，揭示了人类、小鼠和苍蝇基因组中的反转录转座子多样性。本文所描述的发现有助于理解反转录转座子的生物发生机制，并介绍了数十万种新的反转录转座子，这些转座子的功能现在可以被研究。retroSeeker在其他种群和生物体中的进一步应用将加深我们对逆转录转座子复杂性的理解。
从UCSC genome Browser网站(Navarro Gonzalez et al. 2021)以网络文件格式下载了多物种的成对基因组比对，其中三个中心物种是人类、小鼠和苍蝇。有关net文件格式的详细说明，请参阅UCSC Genome Browser描述页面(https://genome.ucsc.edu/goldenPath/help/net.html)和net文件格式的下载页面(https://hgdownload.soe.ucsc.edu/goldenPath/hg38/vsMm10/)。此外，配对基因组比对文件(.net)也可以通过我们的开源管道(https://github.com/junhong-huang/retroSeeker/make_pairwise_alignment_pipeline.pl)生成，其中包括以下步骤:(1)利用lastz软件(https://github.com/lastz/lastz)以参数“-format=axt -ambiguous=iupac - action:target=multiple -strand=both -allocate:traceback=1.99G”生成axt格式的成对基因组比对;(2)利用axtChain软件(http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtChain)以参数“-linearGap=loose”将axt文件格式(.axt)转换为链格式(.chain);(3)通过chainNet软件(http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainNet)使用默认参数将链文件格式转换为net格式(.net)。
随后，我们使用retroSeeker，以net格式文件为输入，识别潜在的逆转录转座子。对同一物种基因组(如人-人网络)采用“自我模式”，对其他物种(如人-鼠网络)采用“不同物种模式”，参数均为默认。通过将逆转录转座子事件的生物学概念(复制和粘贴)映射到具有净格式(空白和填充)的比较基因组学数据中(Weber 2006)， retroSeeker可以识别任何物种基因组中的逆转录转座子。具体工作流程如下:(1)搜索查询基因组内的gap区域(即在参考基因组中存在但不在查询基因组中存在的区域)。(2)在上述空白区域(即在参考基因组和查询基因组中都存在的区域)内搜索填充区域。(3)在填充区的上游和下游端分别向内和向外延伸一定长度的核苷酸。(4)采用动态规划算法识别长度至少为7个核苷酸的TSDs。(5)搜索与3′-TSD相邻的poly(A)区域。(6)得到填充区域的反向互补序列，重复步骤4和步骤5。(7)使用评分系统评估poly(A)和TSD特征，并保留得分较高的反转录转座子候选人。为了确定高可信度的候选者，我们只纳入了TSD长度≥7和poly(A)长度≥5的反转录转座子进行进一步分析。
为了评估retroSeeker的敏感性和特异性，我们使用计算机生成的由60,000,000个随机a /T/C/G碱基组成的基因组进行了模拟研究。最初，我们通过随机将一个原始LINE1序列整合到基因组中来模拟LINE1反转录转座子的亲本。然后，对LINE1序列进行随机修饰，使其包含特定的反转录转座子特征，如poly(A)尾巴长度在5至20个核苷酸之间，靶位重复(TSDs)为7至20个核苷酸，碱基随机为A/T/C/G，插入方向为正链或负链。我们总共生成了1000个反转录转座子，每个转座子被随机整合到模拟基因组中。作为对照，我们还插入了1000个没有TSDs和poly(a)尾部的原始LINE1序列。最后，使用lastz对随机插入反转录转座子序列前后的全基因组进行两两比对，参数为“-format=axt -ambiguous=iupac - action:target=multiple -strand=both -allocate:traceback=1.99G”。此外，我们使用R软件包pROC对retroSeeker生成的反转录转座子评分进行ROC (receiver operating characteristic curve)分析(Robin et al. 2011)。这一分析使我们能够评估retroSeeker在区分反转录转座子和非反转录转座子序列方面的性能。
从UCSC genome Browser网站下载了人类(hg38)、小鼠(mm10)和果蝇(dm6)的基因组序列。人类和小鼠基因注释从GENCODE获得(人类版本39和小鼠版本25)(Frankish et al. 2021)。果蝇基因注释从UCSC (ensGeneV101)获取。通过RepeatMasker识别的重复元件从UCSC Genome Browser网站下载。使用BEDTools软件将所有反转录转座子与典型基因注释相交(Quinlan and Hall 2010)。使用trackBrowser工具可视化单个逆转录转座子(Zhang et al. 2023)，使用Circos可视化逆转录转座子的基因组分布及其与亲本基因的关系(Krzywinski et al. 2009)。
使用CD-HIT (Fu et al. 2012)对高度相似的序列进行聚类，参数设置为“-M 0 -c 0.80 -d 150 -T 30 -s 0.8 -A 0.8 -sc 1”。随后，在基序检测之前，只保留了一个有代表性的序列。我们使用HOMER (Heinz et al. 2010)的findMotifsGenome.pl脚本(参数为“-norevopp -noknown -rna -len 4,5,6,7,8,9 -p 20 -size given -dumpFasta”)和MEME(参数为“-rna -minw 5 -maxw 50 - allow -maxsize 0 -nmotifs 50 -brief 20,000 -p 10 -evt 5”)来重新识别motif。
使用ClustalW 2.1 (Chenna et al. 2003)使用默认参数对组蛋白基因相关的反转录转座子序列进行比对。使用FastTree 2.1 (Price et al. 2010)使用默认参数构建系统发育树，并使用R包ggtree (Yu 2020)绘制系统发育树。使用clusterProfiler (Yu et al. 2012)进行GO分析，只保留相似度小于0.7的GO项。
组织RNA-seq数据下载自ENCODE联盟(consortium, Moore et al. 2020)，癌症RNA-seq数据下载自CCLE (Barretina et al. 2012)。使用STAR (Dobin et al. 2013)软件将clean reads映射到参考基因组(hg38)，参数如下:-outSAMmultNmax -outFilterMultimapNmax 1 -genomeLoad NoSharedMemory -alignEndsType当地-outFilterType正常-outFilterMultimapScoreRange 0 -outFilterMismatchNmax 15 -outFilterMismatchNoverLmax 0.1 -outFilterScoreMin -outFilterScoreMinOverLread 0 -outFilterMatchNmin 18 -outFilterMatchNminOverLread 0.8 -alignIntronMin 5 -seedSearchStartLmax 15 -seedSearchStartLmaxOverLread 1 -seedSearchLmax 0 -alignTranscriptsPerReadNmax 20000 -alignWindowsPerReadNmax 20000 -seedMultimapNmax 20000-seedPerReadNmax 1000 -seedPerWindowNmax 100 -seedNoneLociPerWindow 20 -outSAMtype BAM Unsorted -outSAMmode Full -outSAMattributes All -outSAMunmapped None -outSAMorder Paired -outSAMprimaryFlag AllBestScore -outSAMreadID Standard -outReadsUnmapped Fastx -scoreGapNoncan -8 -scoreGapATAC -8 -scoreGapGCAG -4 -alignSJoverhangMin 15 -alignSJDBoverhangMin 5 - alignendsproude 150 ConcordantPair -scoreGenomicLengthLog2scale -1 -readFilesCommand zcat。然后利用feature - counts (Liao et al. 2014)计算逆转录转座子区域的表达，参数为“-F SAF -s 1 -T 16 -p”。
作为组织特异性的度量，我们使用Tau评分，这是现有方法中计算组织特异性的最佳选择(Kryuchkova-Mostacci and Robinson-Rechavi 2017)，该评分是基于log2 RNA-seq表达数据计算的。Tau的取值范围为0 ~ 1，其中0表示普遍表达，1表示特定表达。Tau值大于0.8的基因被认为是组织特异性的(Kryuchkova-Mostacci and Robinson-Rechavi 2016)。
发现逆转录转座子的一种新的计算方法。图S2。反转录转座子的特性。图S3。通过反转录产生新的基因。图S4。新型的反转录转座子。图S5。组织特异性反转录转座子图谱。图S6。反转录转座子的复杂进化模式。
本研究中使用的物种。
在人类中发现的反转录转座子。
小鼠中发现的反转录转座子。
在苍蝇中发现的反转录转座子。
ccDownload: /内容/ pdf / 10.1007 / s44307 - 023 - 00005 - 5. - pdf
