红毛丹基因组中MITEs转座子的鉴定与演化分析

马誉生，张文萍，明瑞光

（1.福建农林大学生命科学学院，福州 350002；2.福建农林大学海峡联合研究院基因组与生物技术中心，福州 350002；3.伊利诺伊大学厄巴纳香槟分校植物生物学系，美国 伊利诺伊州厄巴纳 61801）

转座子（transposon）是动植物基因组的重要组成成分[1]，是基因组演化过程中不可忽视的一部分[2]，它们可以通过复制或移位从基因组的一个位点转移至另一个位点[3]，从而造成基因组重排，并对基因组大小、染色体结构和附近基因的表达水平等有着不同程度的影响[4]。转座子的含量和种类在不同的生物基因组中存在较大的差异，根据其转座的不同机制，大致可分为ClassⅠ和ClassⅡ两大类[5]。ClassⅠ型转座子以RNA为中间媒介进行转座，可以自身合成逆转录酶然后将其本身反转录为cDNA，以“复制-粘贴”的形式整合到另一个位点，由于其逆转录的特性又被称为逆转座子（retrotransposon）[6]。ClassⅡ型转座子以DNA为中间媒介进行转座，与ClassⅠ型不同的是此类转座子在转座酶的作用下从原位点直接切除并重新整合至新位点，以“剪切-粘贴”的形式进行转座，因此也被称为DNA转座子（DNA transposon）[7]。

微型反向重复转座元件MITEs属于ClassⅡ型DNA转座子中的TIR类转座子，不具备编码转座酶的能力，属于非自主型转座元件（non-autonomous transposable elements）[5]。此类基因组元件首先在植物中发现[8]，并在后续的研究中发现广泛存在于各种生物体基因组中，包括脊椎动物[9]，无脊椎动物[10]，真菌[11]以及病毒[12]。通常MITEs转座子具有以下典型特征：①转座子两端为末端反向互补的重复序列（terminal inverted repeat,TIR）；②在转座子TIR外两侧为具有正向重复的靶位点重复序列（target site duplication,TSD）；③中间为无编码能力的序列，该序列中AT含量较高，有形成二级结构的趋势；④完整片段长度较短，一般在50～800 bp之间[13]。MITEs转座子通过借助与其TIR序列高度相似的对应自主转座子编码的转座酶进行转座，但其拷贝数通常高于其对应的自主转座子[5]。

MITEs转座子在植物基因组中的频繁转座以及高拷贝数等特性对其宿主基因组结构多样性（如等位基因多态性）有着至关重要的影响，甚至影响宿主相关基因的表达从而影响表型[14]。这种现象在桑葚（Morus notabilis）[15]、葡萄（Vitis viniferaL.）[16]、胡萝卜（Daucus carota）[17]和小麦（Triticum-Aegilopsgroup）[18]等多种植物基因组中都有报道，说明了MITEs的存在及转座活动对宿主影响的普遍性。植物MITEs转座子数据库P-MITE于2014年建立并投入使用[19]，该数据库整理了一些植物基因组中MITEs转座子的初步鉴定和分类结果，是植物基因组MITEs研究过程中一个良好的参考平台，由该数据库统计结果（表1）可知，总体上MITEs转座子在各植物基因组中具有较高的拷贝数和含量。

表1 MITEs转座子在植物基因组中的含量[19]Table 1 MITEs content in plant genomes[19]

红毛丹（Nephelium lappaceumLinn）属于无患子科热带水果，其果肉、种子和树干等都具有很高的应用价值，如果肉可做成果酱，树干可用于建筑材料，种仁油可制成肥皂和蜡烛，因此具有极大的商业价值[20]。红毛丹高质量基因组的公布，为红毛丹功能基因组学及基因组演化方面的研究提供了良好的数据基础[21]，而转座子存在于几乎所有的真核生物基因组中，并对其宿主基因组演化、基因表达调控及表型可塑性有着不可忽视的影响[13,22-23]，在红毛丹中还没有对MITEs转座子进行过系统的鉴定和分析，由此本研究利用红毛丹全基因组数据来鉴定、分类并注释红毛丹基因组中的MITEs转座子，以期完善红毛丹基因组的注释并促进其转座子的研究，为红毛丹基因组演化和相关基因功能方面的分析提供关键的数据依据。

1 材料和方法

1.1 红毛丹基因组中MITEs转座子的鉴定与注释

本实验室前期已经完成红毛丹基因组组装及相关基因注释[21]，选用MITE-Tracker软件[24]对红毛丹基因组中的MITEs转座子进行搜寻，设置参数为‘-tsd_min_len 2-tsd_max_len 10-mite_min_len 50-mite_max_len 650’，MITE-Tracker对基因组中具有MITEs转座子结构且满足参数条件的序列提取出来并聚类为不同的家族，由于目前没有软件可以精确辨认所有MITEs转座子的边界，因此后续需要提取出每个家族成员的序列，对每个家族运用MUSCLE-3.7[25]软件进行多序列比对，对边界有缺失的序列扩展前后50 bp，再提取出对应的序列重新进行多序列比对，运用BioEdit-7.0.9.0[26]软件对所得的各家族序列TIR末端进行人工矫正，并按照所总结的MITEs转座子结构汇总表对各个家族进行超家族分类注释（表2），最后运用卷积神经网络深度学习软件DeepTE[27]对注释结果及MITEs序列进行修改验证。运用DAMBE-7软件[28]计算每个MITEs超家族的一致序列（consensus sequence）。完成转座子注释后，计算鉴定出的MITEs序列中腺嘌呤脱氧核糖核酸和胸腺嘧啶脱氧核糖核酸（A和T）的比例，并统计每个超家族在红毛丹基因组中所占的比例。

表2 MITEs转座子超家族分类标准Table 2 Classification criteria for MITEs superfamilies

1.2 MITEs转座子各家族进化树的构建

对各家族的一致序列进行多序列比对后，利用MEGA-X[29]软件构建各个家族代表序列的neighbor-joining进化树，校验Bootstrap参数设置为1 000次，用iTOL（https://itol.embl.de/）在线软件对进化树运行结果进行注释及可视化。

1.3 MITEs转座子插入时间的估算

选取K2P模型（Kimura 2 parameter distances）计算每个MITEs拷贝的插入时间[30]，从而了解MITEs转座子在红毛丹基因组中的扩增情况。首先使用R语言程序包ape（http://ape-package.ird.fr/），结合公式K=-1/2ln(1-2p-q)-1/4ln(1-2q)（q代表颠换位点所占比例，p代表转换位点所占比例）计算每个MITEs完整拷贝与相应超家族一致序列（假定祖先序列）之间的遗传距离，然后运用公式T=K/2r（K为完整拷贝MITEs序列与一致序列之间的遗传距离，r为核苷酸位点替换速率）[31]计算每个转座子的插入时间，在本研究中采用1.38×10-8作为核苷酸位点替换速率。

1.4 MITEs转座子插入位置及密度分布分析

利用软件MITE-Tracker所产生的MITEs转座子注释文件，结合红毛丹基因组注释文件中各元素的长度、位置和正反链等信息，统计每个转座子在基因组中的位置，主要对以下位置的插入进行统计：基因5′端（5 kb以内）、基因3′端（5 kb以内）和基因（内含子和外显子），若MITEs序列任意位点的坐标在上述区域坐标之内，则统计一次插入，最后将所有转座子的插入位置信息及插入次数进行整合。获得红毛丹基因和MITEs转座子在基因组中的坐标信息后，利用R语言程序包gtrellis（https://www.statmethods.net/advgraphs/trellis.html）绘制基因和转座子在红毛丹各染色体上的密度分布图，通过比较分析MITEs转座子在红毛丹染色体上的插入分布模式。

1.5 MITEs转座子在基因附近的插入偏好性分析

提取红毛丹基因组注释文件中每个基因的坐标和距离红毛丹基因5′端和3′端5 kb以内的MITEs的相关注释信息，分别划分为10个区段（每个区段长500 bp），共计20个区段，按照与基因末端坐标由远到近的距离依次统计各区段中MITEs转座子的插入次数，对数据进行可视化从而分析转座子的分布情况与基因距离之间的关系。

1.6 MITEs转座子插入的基因功能富集分析

使用R语言程序包ClusterProfile V3.8（https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html）对MITEs转座子邻近的基因进行GO富集分析，将得到的结果按分子功能（MF）、生物过程（BP）和细胞组成（CC）3部分分类，取最显著的前20个term进行展示比较，然后对MITEs转座子邻近的基因进KEGG富集分析，取最显著的前20个term进行展示比较。

2 结果与分析

2.1 红毛丹基因组中MITEs转座子的鉴定及分类

本研究共鉴定出MITEs转座子完整拷贝588条，全长325 kb，大约占红毛丹全基因组总长（328 Mb）的0.10%，占红毛丹转座子总长度（131.95 Mb）的0.25%，与其他植物基因组中MITEs的占比相比，其含量相对较低。对边界不完整的MITEs拷贝进行人工矫正（图1），根据MITEs转座子结构汇总表的分类标准对软件MITE-Traker所生成的所有家族进行分类。经反复比对及RepeatMasker中重复序列数据库和DeepTE软件训练模型进行深度学习验证，得到红毛丹基因组中MITEs转座子的超家族：PIF/Harbinger、Mutator、Ac-mMITE、hAT和Ginger，其余无法分类的转座子归为Unknown（UN）。各个超家族在红毛丹基因组MITEs转座子总量中所占的比例差异较大（表3），这可能与不同MITEs转座子所对应的自主转座子的拷贝量和活性相关，其中最大的两个超家族为PIF/Harbinger和Mutator，分别有362条和108条完整拷贝，在红毛丹基因组MITEs中的占比分别为64.87%和19.35%，最小的超家族为Ac-mMITE，只有10条完整拷贝，占比为1.80%。

表3 红毛丹基因组中MITEs超家族统计Table 3 Statistics of superfamilies of MITEs in rambutan genome

图1 MITEs转座子手工矫正TIR序列对比图Figure 1 Comparison plots for correction of MITEs TIR sequence

2.2 红毛丹基因组中MITEs转座子的AT含量

红毛丹基因组中MITEs在总体上AT碱基含量的占比为70.63%，其中AT碱基含量比例最高的超家族为Ac-mMITE（78.22%），其次为PIF/Harbinger超家族（76.00%），因此这两种超家族结构稳定性较低，形成二级结构的趋势较高（表4），有利于转座活动。AT碱基数含量较低的超家族为Ginger（37.98%），远远低于总体水平，因此这类超家族形成二级结构的趋势较低。

表4 红毛丹基因组中MITEs转座子AT碱基含量统计Table 4 Statistics on AT bases content of MITEs in rambutan genome

2.3 红毛丹基因组中MITEs的插入时间

红毛丹基因组中MITEs转座子共经历了5次转座“爆发”（图2），其中最早的一次“爆发”大约发生于17～18百万年之前，规模最小，参与的超家族仅有hAT和PIF/Harbinger两种；最近一次的转座“爆发”则发生于1～2百万年之前，MITEs插入数量仅次于规模最大的一次，有将近80个MITEs转座子参与其中，包括Ginger、hAT、Mutator和PIF/Harbinger超家族；规模最大的一次转座“爆发”发生于4～5百万年之前，超过80个MITEs转座子参与，虽然此次转座“爆发”事件规模最大，但参与其中的超家族组成却较为简单，只有Mutator和PIF/Harbinger两种超家族。值得注意的是，在转座子的各个插入时期中几乎每个时期都有超家族PIF/Harbinger的参与，并且在各个时期MITEs的含量上几乎占据主导地位，而Ac-mMITE只参与到少数几个时期的转座事件中且所占比例极小，这与前文中对各个转座子含量分析的结果相一致。

图2 红毛丹基因组中MITEs转座子插入时间Figure 2 The insertion time of MITEs in rambutan genome

2.4 红毛丹基因组MITEs各家族演化树分析

在系统发育演化树中（图3），隶属于同一个超家族的各个亚家族一般具有很近的演化关系，总体上各家族具有较好的聚类效果，这从侧面反映了本研究中所应用的MITEs转座子分类方法的准确率较高。从演化树的分支数来看，PIF/Harbinger超家族的亚家族数最多，主要分为4个类群，这与该超家族作为转座“爆发”事件的主要参与者经历了多轮扩增有关，其次为Mutator，主要分为2个类群，亚家族最少的超家族为Ac-mMITE和Ginger，这些结果也与它们在红毛丹基因组中含量的高低顺序及转座“爆发”事件参与度相一致。未被鉴定出的4个超家族（unknown）中F20、F24、F28和F29分别与PIF/Harbinger、Mutaor、Ginger和hAT超家族具有较近的演化关系，但由于在基因组演化过程中这几个超家族典型的结构完整性已经缺失，无法根据其TIR和TSD序列判断其超家族的类型。

图3 MITEs转座子各家族进化树Figure 3 Phylogenetic tree of MITEs families

2.5 MITEs转座子插入位置及插入密度分布分析

MITEs转座子在红毛丹基因组各个位置上的分布显示出了较大的差异，在基因5′端（5′-flank）和基因3′端（3′-flank）的插入数量最多，分别为136个和137个，说明MITEs转座子的插入主要分布在红毛丹基因两端的位置，且两端插入数量基本持平（表5）。红毛丹基因（gene）内部插入的MITEs总数为35个，其中插入内含子（intron）的MITEs数为35个，插入外显子（exon）的MITEs数为2个，其中有2个MITEs转座子插入区域横跨相邻的外显子和内含子区段。分析红毛丹基因组中MITEs转座子和基因的分布密度（图4），发现MITEs转座子在红毛丹各条染色体上的总体插入比较分散，每条染色体上都有一定数量的插入，但是密度相对较低。MITEs转座子在基因分布较多的区域密度相对较高，说明在红毛丹基因组中该类转座子具有插入到基因密度较高区域的趋势。

表5 红毛丹基因组中MITEs转座子插入分布Table 5 Insertion distribution of MITEs in rambutan genome

图4 基因和MITEs转座子插入密度分布的比较Figure 4 Comparison of density distribution between genes and MITEs

2.6 MITEs转座子在基因附近插入偏好性分析

由红毛丹基因组中MITEs转座子在基因附近插入位置的分布（图5）可知，MITEs转座子在基因两端分布较多。在距离红毛丹基因5′端5 kb范围内，-500～0区段MITEs转座子插入次数只有12次，随着与基因的距离越来越远，在-2 000～-1 501区段插入次数达到峰值，有34次转座子的插入，随后随着距离增加插入次数逐渐减少。在距离红毛丹基因3′端5 kb范围内，0～500区段MITEs转座子插入次数为0，随着与基因距离越来越远次数逐渐增加，在3 501～4 000区段MITEs插入次数达到峰值，有27次转座子的插入，随后插入数逐渐减少。综上，在基因两侧5 kb范围以内，MITEs转座子插入偏好的情况为，距离基因两侧最近的区段，MITEs转座子的插入数最少，随着距离增加在某一区段MITEs转座子的插入偏好会达到峰值，然后又继续减少，在总体趋势上呈现出5 kb以内从离基因较远的某一点开始，距离基因越近，MITEs插入偏好性越低的分布模式，说明距离基因较近的区域，MITEs的转座活动受到了明显的抑制。

图5 MITEs转座子在基因附近的分布情况Figure 5 Distribution of MITEs near the genes

2.7 MITEs转座子关联基因的功能富集分析

在MITEs插入到内部的基因中，有4个含有GO编号（图6），基因编号（表6）分别为：Nl01g01800、Nl01g10720、Nl01g15070和Nl10g00860，这些基因主要参与硫代葡萄糖苷的代谢及生物合成过程（GO:0019757和GO:0019758）、葡糖异硫氰酸盐的代谢及生物合成过程（GO:0019760和GO:0019761）和S-糖苷的代谢及生物合成过程（GO:0016143和GO:0016144）。分子功能方面，这些基因与各种蛋白复合物的结合相关（GO:0005488、GO:0005515、GO:0008017、GO:0008092、GO:0015631、GO:0032403、GO:0044877和GO:0051011），其次为酶活性功能相关基因（GO:0003824、GO:0016853、GO:0016866和GO:0050486），而这些功能对于相关活性物质的代谢及生物合成过程至关重要。在KEGG富集结果中（图6），这些基因（表6）主要富集在次级代谢物生物合成、角质、软木脂蜡质的生物合成、N端多糖的生物合成和硫代葡萄糖苷的生物合成等方面，这些代谢相关通路与上述GO生物过程和分子功能富集结果基本一致。

图6 MITEs转座子插入内部基因功能富集分析结果Figure 6 Gene functional enrichment of inside inserted gene

在MITEs插入到5′端5kb范围内的基因中，有2个含有GO编号（图7），基因编号（表6）分别为：N101g14070和N101g01800，这些基因主要参与生物 体 生 殖（GO:0000003、GO:0022414和 GO:0044702）及发育（GO:0032502、GO:0044767、GO:0009790和GO:0010154等）的相关过程。分子功能方面，这些基因与各种复合物的结合相关（GO:0000166、GO:0005515、GO:0005524和 GO:0008092等），这些功能是生物体生殖发育过程中必不可少的关键要素。在KEGG富集结果中（图7），这些基因（表6）主要富集在次级代谢物的生物合成、类苯基丙烷的生物合成、氮代谢和氨基酸的生物合成等方面，是生物生殖和发育过程中必不可少的环节。

表6 MITEs关联基因NR注释结果Table 6 NR annotation of MITEs inserted genes

图7 MITEs转座子插入5'端侧翼区的基因功能富集分析结果Figure 7 Gene functional enrichment of 5'-flank inserted gene

3 讨论

转座子最初由巴巴拉·迈克林托克于20世纪40年代在玉米（Zea mays）的第9条染色体中发现[32]，然而在很长一段时间里转座子这类重复序列被认为是“垃圾基因”，直到近二十年来随着全基因组测序技术的蓬勃发展为各个物种中转座子的全面鉴定、分类及注释带来了先决条件，越来越多的转座子结构、功能和应用方面的谜题被逐渐揭开[33-34]，这时人们才发现转座子对染色体结构[35]、基因组重排[36]、基因组大小[37]、新基因的形成[38]以及基因表达调控[39]等方面都有着较为显著的影响。

本研究中，MITEs属于非自主型转座子，必须依赖与其同源的自主转座子编码的转座酶进行移位，再插入到基因组的其他位点，因此自主转座子的活性及编码序列完整性决定了其同源MITEs转座子的转座效率及拷贝量。由于转座子受到选择压力，根据转座子的发展史，自主转座元件的活性及结构完整性在基因组演化过程中可能会随着时间发生消退，随后其对应的非自主转座元件拷贝数也开始逐渐减少[40]，红毛丹中MITEs的含量相对较低，这可能与其基因组中自主DNA转座子的活性[5]和其在双子叶基因组三倍化事件之后没有再发生过另外的全基因组复制事件有关[21]。在所有MITEs超家族中，PIF/Harbinger的相对含量占据主导地位且在各个MITEs的转座“爆发”时期贡献度最高，而AcmMITE相对含量最低且在MITEs的转座“爆发”的各个时期贡献度最低，因此可以推测PIF/Harbinger的自主转座子在近期转座子扩增的过程中仍具有较高的活性和结构完整性，而超家族Ac-mMITE的扩增活性最低，其大部分同源的自主转座子可能已经失去活性及结构完整性。另外，MITEs转座“爆发”事件在一定程度上可以反映植物基因组演化历史中的重要事件，并且在宿主基因组可塑性和转录应答的改变上可能有一定的贡献[41]，比如Chen J.等在水稻的重组近交系中研究由转座元件爆发所产生的基因组多样性现象时发现Ping/mPing转座子会影响其宿主基因组的多态性[35]，说明转座子对基因组的演化具有重要意义。由此可见，MITEs转座子在红毛丹基因组中的转座活动必然会对基因组结构和演化活动造成一定的影响。

在MITEs插入内部的基因功能富集分析中，MITEs插入到某些通路相关基因的内部会改变该基因的结构，从而最终影响到通路本身，甚至改变宿主的某些性状，Wei L.等在水稻的赤霉素及油菜素甾醇合成途径中发现，如果该通路的关键基因有MITEs的插入会影响其旗叶夹角和水稻株高等相关性状[42]。在蔬菜和水果中有大量以葡萄糖苷形式存在的非挥发性风味前体物质，是蔬菜水果增香的重要大分子[43]，而MITEs转座子在这些化合物通路相关基因的插入可能会改变葡萄糖苷、葡糖异硫氰酸盐和S-糖苷的代谢及生物合成相关基因的编码序列，最终对红毛丹相关风味等重要农艺性状造成影响。在MITEs插入到5′端基因的功能富集分析中，有两个基因主要与生物体的生殖发育等过程相关，而基因5′端是基因顺式调控元件所在的区域，MITEs转座子在这些区域的插入可能会对红毛丹生殖发育相关的表型造成影响，如Y.Lee等在玉米中研究一种转座子监督机制保护雄性玉米植株生育能力的过程中发现，转座子的转座活动会增加基因组的不稳定性，从而可能产生雄性不育的表型[44]。在MITEs插入位置分析中，插入到基因内部的MITEs比例很低，其中大多数MITEs插入到了内含子当中，会在被插入基因的结构上产生一定程度的改变，甚至影响可变剪切，少数MITEs插入到了外显子中，在这些基因的功能方面可能会造成直接影响。

MITEs在红毛丹基因组中插入密度与基因分布密度具有较高的相关性，而且大量MITEs插入到了红毛丹基因的5′端和3′端附近，这种大量在基因附近的转座活动可能会在基因表达调控网络的水平上对基因功能产生影响，最终可能改变宿主表型，产生表型可塑性的改变，如E.Butelli等在血橙的表型研究中发现当Copia-like反转座子插入Ruby基因编码区上游并受到一定寒冷诱导时，会使其产生红色果肉的表型[45]。同样，红毛丹中MITEs转座子在这些区域的转座事件可能在转录和后转录水平上对其基因的表达调控有一定的影响，至于影响水平的高低，是否会改变相关的性状，需要后续对其进行更加详细的分析和实验验证，Chen L.等学者在菠萝基因组及其驯化过程的研究中发现，MITEs转座子在菠萝基因附近的插入位置具有高度多样性，推测这一现象在菠萝驯化过程中可能会通过体细胞突变作为菠萝新性状形成的主要来源之一[46]。随着越来越多MITEs转座子的鉴定及相关分析的细化，MITEs在其宿主基因组的演化及相关性状的多样性上所扮演的角色会逐渐明朗，对这些信息及相关效应的充分利用可以为其宿主基因组演化及性状改良等方面的研究带来更多的途径。