重复序列是真核动物基因组中的重要组成部分，例如人基因组中超过50%的序列都是重复序列 (Lander et al., 2001)，而在狗基因组中重复序列约占整个基因组的36.1%(Li et al., 2010)。重复序列主要分为两大类：串联重复序列 (tendam repeat) 和散在重复序列 (interpersed repeat)。其中，串联重复序列包括微卫星、小卫星序列等；散在重复序列又称可移动元件，或转座子元件 (transposable element，TE)。根据转座方式不同，TE可分为DNA转座子 (DNA transposon) 和逆转录转座子 (retrotransposon)，前者通过“剪切-粘贴”的机制，将DNA从基因组剪切后直接插入基因组中一个新的座位；后者则先将DNA转录出一段RNA，RNA经反转录酶重新合成DNA，再插入基因组某个新的位置，即通过“拷贝-粘贴”的方式在基因组中进行移动。常见的反转录转座子类别有长末端重复序列 (LTR)、长散布重复序列 (LINE)、SVA和短散布重复序列 (SINE) 等。大量转座子能在基因组中移动并进行扩增，是导致不同动物基因组结构多样性的重要原因；同时研究还发现它们也广泛介导机体内多种生物学过程，包括导致基因突变、调节基因表达、导致染色体断裂、引起mRNA的可变剪切以及多种肿瘤等疾病等 (Belancio et al., 2009；Copeland & Jenkins，2010；Lupski，2011；Ray & Batier，2011)。因此，系统研究基因组中的这些可移动元件对于阐明基因组的结构与功能具有重要意义。

然而受限于测序技术和组装技术等困难，鉴定和注释基因组中的转座子等重复序列一直是基因组学研究的一个挑战。由于目前基于二代测序技术产生的序列读长 (reads) 较短，而基因组序列采用k-mer (de Bruijn graph) 算法进行组装，高度相似的重复序列可能会被压缩到一起，从而影响对后续重复序列的识别。用于识别不同重复元件的方法或软件目前已超过30种，大体可将它们分为两大类：1) 序列比对法，基于特定结构来鉴定重复序列。该方法需要构建重复序列的文库，文库中包含很多来自不同物种某一重复序列的一致性序列 (Repbase数据库)，然后通过识别与已知重复序列相似的序列，从而对其进行分类。此方法最经典、最流行的软件是Greedier (Li et al., 2008) 和RepeatMasker (RM)(Smit et al., 2016)，Greedier在重复序列鉴定的敏感性上比RM稍高，但是重复序列的鉴定率只有RM的一半左右。此外，MaskerAid (Bedell et al., 2000) 和Censor (Kohany et al., 2006) 亦是基于该原理进行鉴定。2) 重复序列的从头预测法 (de novo)。该方法是利用重复序列或转座子自身的序列或结构特征，构建从头预测法或软件 (模型) 对序列进行识别，其优点在于能够根据TE自身的结构特征进行预测，不依赖于已有的转座子数据库，能够发现未知的TE。常见的从头预测法有Recon (Bao & Eddy，2002)、RepeatScout (RS)(Price et al., 2005)、Piler、LTR-finder、ReAS (Li et al., 2005) 等。

大熊猫Ailuropoda melanoleura是我国生物多样性保护的旗舰物种，其独特生物学和行为学特征一直是研究热点。为阐明其独特的遗传机制，2010年我国研究人员公布了一只雌性大熊猫的全基因组序列，发现大熊猫基因组与其他食肉目Carnivora动物相似，为2.25 G，其中约36.2%是由TE构成 (Li et al., 2010)，与狗基因组 (36.1%) 和猫基因组 (32.1%)(Pontius et al., 2007) 中TE的比例相似。然而大熊猫的2个Can-SINEs亚家族 (SINEC_b1和SINEC_b2) 的拷贝数比狗的高2.5倍，而其他一些Can-SINEs亚家族则显著低于狗，但相关文章只基于单一方法对大熊猫的TE做了初步统计，并未详细报道大熊猫基因组各种重复序列的类型，亚家族的数量、长度、密度、分化率等。大熊猫基因组重复序列的详细特征，尤其是那些在大熊猫与狗分化之后才插入大熊猫基因组的TE，它们究竟属于何种类型，如何塑造大熊猫基因组，以及它们对大熊猫的生长发育、繁殖代谢是否存在影响都不得而知。

对重复序列的鉴定与注释是研究其特征和功能的重要前提。本研究选择了2种常用的重复序列研究策略：基于序列比对的RM和从头预测的RS，分别对大熊猫基因组中的重复序列进行鉴定与注释，比较了2种方法的鉴定效果，为深入研究重复序列对大熊猫基因组结构和功能的影响提供重要的参考。

1 研究方法 1.1 基因组序列的下载

从GiGa (http://gigadb.org/) 数据库下载大熊猫全基因组序列，其组装水平为scaffold，总的序列 (scaffold) 数量为81 469，总长度为2 299 498 912 bp，N50约1.3 Mbp。

1.2 基于RM注释大熊猫重复序列

下载RepeatMasker (4.0.5)(Smit et al., 2016) 及其默认数据库Repbase (2015.8.7)(Smit et al., 2016) 进行本地安装。选择crossmatch作为其搜索比对算法，选择食肉目作为查询数据库类别，设定阈值为225，选择-low、-s参数，对整个大熊猫基因组序列进行重复序列重注释。

1.3 基于RS注释大熊猫重复序列

下载RepeatScout 1.0.5(http://bix.ucsd.edu/repeatscout/) 并进行本地安装。因RS无法处理大于1 G的序列文件，故将大熊猫基因组的序列分为4个部分，各部分分别包含490条、489条、805条和79 688条scaffolds。对过滤Perl脚本compare-out-to-gff.prl进行优化，降低循环次数以提高效率。选择k-mer值为17，过滤阈值为20，构建大熊猫重复序列的一致性序列。以RS重复序列为查询库、RM作为重复序列分类工具对大熊猫基因组进行注释，而后使用本地Perl脚本交叉注释，得到RS对重复序列分类注释的结果。

1.4 注释重复序列比较

使用本地Perl脚本统计2种方法重复序列数量及分布情况，R软件统计主要转座子的长度分布及数量，Excel绘图比较2种方法的差异。

2 结果 2.1 RS构建各种重复序列一致性序列的注释

基于现有数据库对RS注释的各种重复序列进行鉴定。RS总计构建出了3 400条重复序列的一致性序列，其中超过80%能够被已知的重复序列注释 (表 1)。在这些被注释的序列中，LINE占构建的一致性序列的比重最大 (47.23%)，LTR次之 (27.87%)，Small RNA和SINE较少，占比分别为19.00%和18.36%。而在注释到的LINE中，绝大部分都是LINE1，只有相对少的LINE2和L3/CR1。这表明了RS能够成功构建各种主要的重复类型，尤其是当一种重复序列的结构明显、数量众多时，其构建的该类重复序列的一致性序列也多，如LINE、LTR、SINE、Small RNA等。而DNA转座子在大熊猫基因组中的数量相对较少，结构并不明显，RS构建的核心序列也就相对较少，只有5.40%。此外还有20%重复序列的一致性序列不能被现有的重复序列数据库类型注释。

2.2 2种方法对大熊猫重复序列的注释统计

大熊猫总共81 647条scaffolds，采用RM与RS分别注释到3 474 463个、2 728 748个重复序列，分布于39 064条、25 035条scaffolds上。以RM作为重复序列分类工具对RS的结果进行了分类，发现RS的结果中能被注释的已知重复序列仅1 922 054个 (表 2)，即有约30%的重复序列能被RS识别，但无法注释到现有的已知重复序列类型。由于这部分重复序列与现有重复元件不同，主要针对注释的已知重复元件进行了比较。

RM鉴定大熊猫基因组中数量最多的重复元件是LINE，达到1 368 483个，占总TE的39.46%，LINE元件总长度占整个基因组的比例也最高，为21.13%。SINE元件数量仅次于LINE，共计1 259 550个，占TE总数的36.5%，其总长度占基因组的8.53%，其中包括MIR家族 (518 019个)。此外，RM还鉴定出413 355个LTR元件，380 310个DNA转座子和34 540个Small RNAs等 (表 3)。经RS鉴定的已知重复序列则包括1 018 913个LINE元件，其数量与RM鉴定结果基本相当；但RS仅鉴定了427 086个SINEs，远远少于RM鉴定的结果 (仅为其33.9%)，总长度仅占基因组的2.91%；此外RS注释了334 194个LTR元件和139 277个DNA转座子等 (表 4)。虽然2种方法均能够注释到主要的一些TE类型，但RS注释到的TE数量明显少于RM，尤其是SINE和DNA转座子的数量差异较大。

统计TE密度最高的前20条大熊猫scaffolds (图 1：A，B)，发现2种方法注释的结果差异很大，20条scaffolds各不相同。RS鉴定出的是scaffold2381(2.96/kbp)，20条scaffolds的TE平均密度为2.19/kbp；而RM注释出的是scaffold6473(6.95/kbp)，TE平均密度为3.53/kbp。从20条scaffolds的长度来看，RS鉴定较长的scaffold上TE密度较高，18条scaffolds都大于10 000 bp；相反，RM鉴定较短的scaffold上TE的密度高，密度最高的scaffold6473仅长719 bp，仅2条scaffolds大于10 000 bp。结果还显示4种类型的TE在各scaffold上的分布并不均匀。大多数scaffold上LINE的密度都高于其他类型，但在某些scaffold上，则表现为SINE或LTR元件的富集情况 (图 1：A，B)。此外，2种方法注释的LINE和LTR的密度差距不大，而RM注释的SINE和DNA转座子的密度分别是RS的3.8倍和5.1倍。

统计2种方法注释的大熊猫TE的分化率，即根据每个TE序列累积的突变位点计算其与一致性序列的差异 (图 2：A，B)。分化率越大，累积的突变越多，意味着该TE插入基因组的时间越长，是更古老的转座子。2种方法鉴定TE虽然在总数上存在差异，但对4类TE的分化率统计结果却相似，每种类型TE的分化率都有高有低，既有插入时间长的，也有新近插入基因组的。然而在分化率＜10%的区域，主要是LINE或SINE逆转座子，大部分的LTR和DNA转座子的分化率都集中分布在＞20%的区域，说明新近整合到大熊猫基因组的主要是一些LINE或SINE。

2.3 2种方法注释大熊猫4种TE亚家族的比较

在Repbase数据库中，食肉目动物的SINE元件包含33个亚家族，RM在大熊猫基因组中鉴定了所有的33个亚家族，而RS仅注释到其中的14个，前者为后者的2.36倍。RM鉴定＞100 000个拷贝的亚家族有5个，其中SINEC1_AME和SINEC1B_AME的数量均＞200 000，占总SINEs的18.9%和17.2%，是数量最丰富的SINE元件；拷贝数＜1 000的亚家族15个。相似的，RS鉴定的SINE亚家族中，SINEC1B_AME最多 (222 227个)，甚至超过RM的结果，占总SINEs的52%，其次为SINEC1_AME，占23.5%。与RM的结果相比，RS对SINEC1B_AME (17.2% vs. 52%) 及SINEC_B1(2.0% vs. 7.1%) 亚家族的鉴定更好，而对MIRB (15.4% vs. 0.58%) 和MIR (11.9% vs. 0.49%) 等亚家族的鉴定较差 (图 3：A)。

RM在大熊猫基因组中鉴定了LINE元件已知的所有143个亚家族，而RS仅注释到其中58个，但二者鉴定的LINEs总数相差不大。RM鉴定数量最多的都是L2(L2A、L2C和L2B)，三者共占总LINEs超过23.6%；其次是L1-1_AME亚家族；有47个LINE亚家族的拷贝数量都＜1 000。相反，RS对L2的注释较差，L2A仅有991个，占总数的0.097%，L2B也只有358个，没有L2C (而RM鉴定L2C＞100 000个)；但RS能较好地识别L1家族，数量排名前20的都是各种L1的亚家族，其中最多的是L1_CANID2(142 408个)，而RM仅鉴定有9 055个该亚家族拷贝；在RS识别的58个亚家族中，40个亚家族拷贝数＞1 000，它们都属于L1家族 (图 3：B)。

LTR元件共包含341个亚家族，是亚家族最丰富的TE类型。2种方法鉴定出的大熊猫基因组中LTR总数基本相当，RM鉴定了所有341个亚家族，而RS仅注释到其中的120个。RM鉴定最多的2个亚家族是MLT1D和MLT1B，前7名的都是MLT1家族的成员。341个亚家族中，仅有97个亚家族拷贝数＞1 000，其余都是数量少的小亚家族。RS鉴定数量最多的仍然是MLT1B (47 798个) 和MLT1D (44 050个)，分别是RM鉴定数量的3.21倍和2.77倍。然而RS对某些MLT1亚家族的鉴定效果却较差，如无法识别MLT1K、MLT1L、MLT1T亚家族，而RM鉴定这3个亚家族的数量都＞10 000。在RS注释的120个LTR亚家族中，拷贝数＞1 000的有52个，其余68个亚家族的拷贝数都很少 (图 3：C)。

RM在大熊猫基因组中鉴定了已知DNA转座子的所有288个亚家族，而RS仅注释到其中的28个，仅9.7%的亚家族被成功注释，是4种TE中比例最低的。RM鉴定最多的亚家族是MER5A (31 040个) 和MER5B (21 562个)，拷贝数＜1 000的小亚家族193个。与RM相比，RS对某些MER亚家族的鉴定效果好，MER20、MER5A、MER5A1、MER3都超过10 000个拷贝，有的比RM的结果还高。但RS对某些MER亚家族却无法识别，如MER5B、MER103C拷贝数都为0，而RM鉴定二者分别有21 562个、7 717个拷贝，此外，它也无法识别MAMTIP2B亚家族 (图 3：D)。

2.4 2种方法注释4种TE亚家族长度的比较

为比较2种方法注释大熊猫4类TE的长度，统计了每个亚家族内所有拷贝的平均长度 (可发邮件向作者索取)。总的来说，2种方法注释TE的长度与其本身长度具有一致性，若这类TE本身很长，则2种方法注释到其长度也相对较长，反之亦然。LINEs平均长度最长，SINEs最短。同类TE的长度差异也很大，最长的LINE亚家族是L1_CANID2(1 238.0 bp)，最短为X7D_LINE (65.6 bp)。2种方法均鉴定到SINE元件中平均长度最长的是SINEC1B_AME和SINEC1_AME亚家族，而它们也是鉴定数量最多的SINE亚家族。而2种方法都注释LTR22_AME亚家族是大熊猫基因组中较长的LTR逆转座子，但是其拷贝数却较少。

RM注释的各种TE平均长度均长于RS，也就是说RM注释的结果更接近TE的全长。RS注释了14个SINE亚家族、58个LINE亚家族、120个LTR亚家族和28个DNA转座子亚家族的平均长度分别为96.4 bp、170.5 bp、132.7 bp和119.5 bp；而RM鉴定的同样的TE亚家族平均长度分别为133.8 bp、479.8 bp、349.3 bp和252.3 bp (图 4)，分别为前者的1.39倍、2.81倍、2.63倍和2.11倍。而对长度越完整的转座子的鉴定亦越准确，如LTR22_AME亚家族，2种方法统计其平均长度较接近 (443.0 bp和401.3 bp)，鉴定其拷贝数也相似 (分别为717和643)。而RS鉴定各MIR亚家族的长度远低于RM (如MIR亚家族长度仅为后者的50%；MIRB仅为后者的61.1%)，因而RS鉴定MIR家族的数量远远少于RM。对L2家族的鉴定结果亦是如此，RS鉴定L2A、L2B的平均长度仅为RM的27%、21%，而RM鉴定的数量为RS的140倍、230倍。

3 讨论

此前关于食肉目动物TE的了解，绝大多数来自模式动物——狗或猫。比较家犬和贵宾犬基因组发现，两者包含大量仍然具有活跃转座能力的CanSINE元件 (Wang & Kirkness，2005)，它们在不同基因组中的转座活动可能导致皮肤的不同花纹或某些疾病 (Clark et al., 2006；Walters-Conte et al., 2011)。而在猫亚目Feliformia中，SINE元件至少可分为2大类，其中一类在不同猫科Felidae动物中也仍然非常活跃 (Walters-Conte et al., 2014)。本研究系统研究了大熊猫基因组中的重复序列，尤其是其中的TE，详细阐明了TE的组成、类型、数量、亚家族、长度分布和分化率等信息。大熊猫基因组在TE的组成与类型上与狗和猫相似，都包括4大类TE，即SINE、LINE、LTR和DNA转座子。Li等 (2010)曾推测大熊猫基因组存在约70 Mbp的TE序列，它们的分化率＜10%，可能是大熊猫特异的一些TE。这与本研究结果一致，2种鉴定方法都注释了分化率＜10%的转座子，它们都只属于SINE和LINE。深入研究这些仍然活跃的转座元件，将有利于阐明它们对大熊猫疾病、基因组稳定性和多样性的影响和作用机制。

由于不同动物的TE类型及丰富度差异极大，利用现有的重复序列注释方法研究非模式动物的重复序列非常困难。如灵长目Primates动物数量最大的Alu元件、鸟类中最丰富的CR1元件，在食肉目动物中几乎没有或极少，因此很难找到一个统一有效的重复序列研究方法。同时，现有的重复序列数据库也无法识别非模式动物中尚未被认识的新的重复类型，所以基于现有数据库去注释非模式动物时往往会存在许多问题。另一方面，基于从头预测法的方式对非模式动物重复序列进行注释，有可能发现新的TE类型，但这种方法可能存在大量的假阳性，即产生一些非TE的序列 (Bergman & Jenkins，2007)。基于此，本研究综合采用了基于数据库的方法和从头注释的方法，更有利于全面地了解大熊猫基因组中转座子的特征。结果显示，RM与RS都能鉴定大熊猫基因组中几种主要的TE类型，但它们在TE亚家族的类型、数量和平均长度上存在明显差异。RM注释到的TE数量在绝大部分亚家族中均多于RS，在某些亚家族中较为接近，而在一些亚家族中则少于RS，而RS注释的TE的亚家族类型及平均长度均小于RM。利用RS (Fernando et al., 2012) 仅注释了大熊猫1条scaffold292上的TE，其结果也远少于RM注释的数量，这与本研究结果类似。这表明了RS所采用的k-mer打断多序列比对重构重复序列识别算法可能存在一定的局限性，难以重构出所有的转座子序列。但我们发现这种算法在对某些TE亚家族上的识别能力要强于RM，如结构十分明显的LINE1、LTR逆转座子，而对DNA转座子和某些SINE元件则非常差。其次，RS还存在一个内存限制的问题，不能够同时处理所有大熊猫基因组的序列，这可能也是造成其无法识别所有大熊猫转座子的原因。此外，RS有较为严格的对于重复序列拷贝数的阈值要求，低于这个阈值的序列会被过滤，这也可能造成一些低拷贝的TE亚家族被过滤掉而没能保留下来进行分类。最后，由于组装的困难，大熊猫基因组中重复序列会产生很多碎片，因而难以被识别，而从头预测法对序列结构的完整性及计算资源的要求又比较高，这可能是导致RS的鉴定不及RM的重要原因。虽然在已知TE类型的鉴定上，RS的表现不及RM，但本研究发现RS构建的大熊猫重复序列的一致性序列中，有20%并不属于现有的重复序列类型，排除部分假阳性的结果，其中也可能包含一些大熊猫特有的TE类型。因此通过整合多种算法优化TE识别率，深入研究潜在的新的重复序列，将有可能全面揭示大熊猫基因组在重复序列上与其他食肉目动物的重要差异。