劉震 張樹林 史瑩慧 彭仁海

（安陽工學(xué)院生物與食品工程學(xué)院，安陽 455000）

長末端重復(fù)序列（Long terminal repeat，LTR）反轉(zhuǎn)錄轉(zhuǎn)座子是真核生物基因組中普遍存在的一類可移動(dòng)的DNA序列，它們以RNA為媒介，通過“復(fù)制-粘貼”的形式在基因組中不斷自我擴(kuò)增。LTR反轉(zhuǎn)錄轉(zhuǎn)座子的插入和刪除對基因組的進(jìn)化和物種的環(huán)境適應(yīng)能力都具有重要的意義［1］，此外，LTR反轉(zhuǎn)錄轉(zhuǎn)座子、基因組重復(fù)和多倍化是導(dǎo)致植物基因組擴(kuò)增和收縮的三個(gè)主要原因［2，3］。LTR反轉(zhuǎn)錄轉(zhuǎn)座子根據(jù)內(nèi)部編碼基因（gag，pol和int）排列順序的不同，分為Copia和Gypsy兩類主要的超家族［4］。LTR末端不編碼蛋白質(zhì)，但包含轉(zhuǎn)錄的起始信號(hào)和終止信號(hào)，內(nèi)部編碼區(qū)gag基因編碼的蛋白負(fù)責(zé)反轉(zhuǎn)錄轉(zhuǎn)座子RNA的成熟和包裝，pol基因編碼反轉(zhuǎn)錄酶和 RNAse H，INT（Integrase）基因編碼整合酶［4］。

LTR反轉(zhuǎn)錄轉(zhuǎn)座子的插入除了可以導(dǎo)致基因組的膨脹外，更重要的是影響插入位點(diǎn)及其相鄰基因的表達(dá)［5，6］。LTR反轉(zhuǎn)錄轉(zhuǎn)座子能夠完成自我轉(zhuǎn)錄是因?yàn)槠浔旧砗修D(zhuǎn)錄所需的調(diào)控元件。當(dāng)LTR反轉(zhuǎn)錄轉(zhuǎn)座子插入到基因編碼區(qū)就可能導(dǎo)致該基因轉(zhuǎn)錄成不完整RNA序列，進(jìn)而不能被翻譯成完整的肽鏈，或者該RNA失去其調(diào)控能力；當(dāng)LTR反轉(zhuǎn)錄轉(zhuǎn)座子插入到基因附近區(qū)域時(shí)，其序列內(nèi)的調(diào)控元件將會(huì)發(fā)揮作用，并影響附近基因的表達(dá)［7］。

棉纖維是紡織工業(yè)中天然纖維的主要來源。中國種植的四倍體棉花主要包括陸地棉和海島棉兩個(gè)栽培種，陸地棉豐產(chǎn)性好，而海島棉不僅纖維品質(zhì)優(yōu)良，而且是鹽堿地的主要栽培作物。隨著國內(nèi)紡織技術(shù)的發(fā)展和人們生活水平的提高，海島棉的需求量也逐年增加。

四倍體海島棉的基因組為2.57 G，大約69.11%為重復(fù)序列，其中A亞組重復(fù)序列占73.5%（1 098 Mb），D 亞組重復(fù)序列占 63.5%（541.6 Mb）［8，9］，而A基因組的亞洲棉與D基因組的雷蒙德氏棉中轉(zhuǎn)座子的含量則分別為57.0%和68.5%［10］。在四倍體棉種轉(zhuǎn)座子起源的相關(guān)研究中，發(fā)現(xiàn)更多的LTR反轉(zhuǎn)錄轉(zhuǎn)座子起源于A基因組，或者說，四倍體A亞組中有更多的轉(zhuǎn)座子拷貝［11，12］。此外，雷蒙德氏棉的轉(zhuǎn)座子數(shù)據(jù)庫也已經(jīng)公布［13］。這些數(shù)據(jù)為深入研究海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子提供了絕佳的機(jī)遇。本研究首先綜合多種不同的方法挖掘了海島棉基因組中的LTR反轉(zhuǎn)錄轉(zhuǎn)座子，然后對這些轉(zhuǎn)座子進(jìn)行了家族分類、周邊基因的功能富集、數(shù)量分布和統(tǒng)計(jì)分析。本研究對海島棉基因功能分析和基因組進(jìn)化有重要的參考價(jià)值。

1 材料與方法

1.1 材料

海島棉的基因組序列、基因注釋和GO注釋文 件 均 從 COTTONGEN（https://www.cottongen.org/）下載。

1.2 方法

1.2.1 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子的挖掘與分類 分別通過依據(jù)LTR反轉(zhuǎn)錄轉(zhuǎn)座子結(jié)構(gòu)特征的工具LTR_STRUC［14］和 LTRharvest［15］；依據(jù) LTR 反轉(zhuǎn)錄轉(zhuǎn)座子重復(fù)特征的工具PILER［16］；以及綜合性工具RepeatModeler［17］搜尋海島棉基因組中的LTR反轉(zhuǎn)錄轉(zhuǎn)座子。使用REPCLASS［18］軟件將上述結(jié)果序列歸類到相應(yīng)的超家族，并將相同的超家族合并，之后再與已知重復(fù)序列數(shù)據(jù)庫Repbase［19］進(jìn)行進(jìn)一步的合并。利用CD-HIT［20］去除合并結(jié)果中的冗余序列，得到海島棉特異的LTR反轉(zhuǎn)錄轉(zhuǎn)座子序列庫，最后利用RepeatMasker［21］注釋海島棉基因組中的LTR反轉(zhuǎn)錄轉(zhuǎn)座子，由同一參考序列注釋到的一組序列被認(rèn)為是一個(gè)家族。分析過程中，要求LTR反轉(zhuǎn)錄轉(zhuǎn)座子序列最短為80 bp，每個(gè)LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族在基因組中有100次以上的重復(fù)拷貝。

1.2.2 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子的數(shù)量與分布 通過Perl腳本從RepeatMasker結(jié)果文件中收集LTR反轉(zhuǎn)錄轉(zhuǎn)座子在海島棉基因組的數(shù)量和位置，并利用gff注釋文件的數(shù)據(jù)收集基因的數(shù)量和位置。統(tǒng)計(jì)分析染色體中每100 kb范圍內(nèi)的LTR反轉(zhuǎn)錄轉(zhuǎn)座子與基因的數(shù)量，并通過Circos［22］繪制分布圖。

1.2.3 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊基因的GO富集分析 查找海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子上、下游20 kb范圍內(nèi)的基因，利用基因組GO注釋文件確定這些基因的GO注釋條目，并使用WEGO［23］（http://wego.genomics.org.cn/）進(jìn)行富集分析。

2 結(jié)果

2.1 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子的數(shù)量分析

圖1 海島棉A亞組和D亞組共有和特異的LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族

海島棉為四倍體棉種，A亞組和D亞組各包含13條染色體。數(shù)據(jù)結(jié)果（圖1）表明海島棉基因組中共包含2 018個(gè)100拷貝以上的LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族，其中1 930個(gè)家族共同存在于A亞組和D亞組中，84個(gè)家族為A亞組特異，4個(gè)家族為D亞組特異。A亞組共包含274 360個(gè)LTR反轉(zhuǎn)錄轉(zhuǎn)座子拷貝，D亞組則包含209 415個(gè)拷貝，因此，LTR反轉(zhuǎn)錄轉(zhuǎn)座子在A亞組中的拷貝要比D亞組多，這一特征對于高拷貝數(shù)的LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族表現(xiàn)的更為明顯（圖2）。此外，從圖2還可以大致看出，一個(gè)LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族在較大的染色體上有較多的拷貝數(shù)，在較小的染色體上有較少的拷貝數(shù)，LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族特異分布在少數(shù)染色體上的情況并不明顯。

圖2 不同拷貝數(shù)的LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族在海島棉各染色體上的分布

海島棉基因組每Mb序列平均包含242.21拷貝的LTR反轉(zhuǎn)錄轉(zhuǎn)座子，通過計(jì)算發(fā)現(xiàn)海島棉基因組中每Mb染色體包含的LTR反轉(zhuǎn)錄轉(zhuǎn)座子的拷貝數(shù)與染色體大小的皮爾森相關(guān)系數(shù)為-0.52，表現(xiàn)出中度負(fù)相關(guān)的關(guān)系（表1）。

2.2 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子的分布特征

海島棉基因組中的LTR反轉(zhuǎn)錄轉(zhuǎn)座子整體在A亞組染色體的后端部分分布較少（4號(hào)染色體除外），而在D亞組染色體則相對分布均勻（圖3）。Copia超家族在染色體的起始端有一個(gè)明顯的波峰，這一特點(diǎn)A亞組和D亞組類似，同時(shí)，Copia超家族在A亞組染色體的后端部分分布較少，而D亞組染色體則沒有該特征。Gypsy超家族整體在染色體中部分布多，兩端分布少，進(jìn)一步比較A亞組和D亞組可以發(fā)現(xiàn)，Gypsy超家族在A亞組后半部分分布的更少。因此，A亞組染色體后端LTR反轉(zhuǎn)錄轉(zhuǎn)座子分布較少是這兩類主要超家族共同的分布特征。

比較海島棉基因組中基因與LTR反轉(zhuǎn)錄轉(zhuǎn)座子的分布發(fā)現(xiàn)Gypsy超家族與基因的分布呈近似反比關(guān)系，而Copia超家族與基因則沒有明顯的數(shù)量分布關(guān)系。

表1 海島棉各染色體的大小及LTR反轉(zhuǎn)錄轉(zhuǎn)座子的分布密度

圖3 LTR反轉(zhuǎn)錄轉(zhuǎn)座子在海島棉基因組上的分布

2.3 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊基因的富集分析

分別對海島棉A亞組與D亞組LTR反轉(zhuǎn)錄轉(zhuǎn)座子上下游20 kb范圍內(nèi)的基因進(jìn)行了GO富集分析。結(jié)果共涉及GO分類體系中的細(xì)胞組件、分子功能和生物過程3個(gè)大類別中的9、12和19個(gè)小類別。A亞組和D亞組LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊分別有47 979和55 880個(gè)基因具有GO注釋。從圖4中可以看出兩類超家族的富集情況基本相同，在細(xì)胞組件中，涉及較多的條目依次是細(xì)胞（Cell GO：0005623）、細(xì) 胞 組 分（Cell part GO：0044464）、細(xì)胞器（Organelle GO：0043226）和大分子復(fù)合物（Macromolecular complex GO：0032991） 等；在分子功能方面，主要富集在結(jié)合活性（Binding GO：0005488）和催化活性（Catalytic activity GO：0003824）等類別中；而在生物學(xué)過程中，涉及較多的條目依次是代謝過程（Metabolic process GO：0008152）、細(xì) 胞 過 程（Cellular process GO：0009987）、生 物 調(diào) 節(jié)（Biological regulation GO：0065007）、定位（Localization GO：0051179）、建立定位（Establishment of localization GO：0051234）和色素（Pigmentation GO：0043473）等。將二倍體亞洲棉（A組）、雷蒙德氏棉（D組）與四倍體海島棉（AD組）進(jìn)行比較，發(fā)現(xiàn)LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊基因的GO富集情況類似，只是在基因數(shù)量和百分比方面有差別（數(shù)據(jù)未發(fā)表）。

圖4 海島棉LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊基因的GO富集分析

3 討論

LTR反轉(zhuǎn)錄轉(zhuǎn)座子是植物基因組的重要成分，是推動(dòng)基因組大小變異和進(jìn)化的重要因素［2］。精確而完整的LTR反轉(zhuǎn)錄轉(zhuǎn)座子注釋對研究基因組大小變異和進(jìn)化具有非常重要的意義。從基因組中挖掘轉(zhuǎn)座子序列的算法主要有三類：依據(jù)轉(zhuǎn)座子的結(jié)構(gòu)特征、依據(jù)轉(zhuǎn)座子在基因組中的重復(fù)特征和依據(jù)已知轉(zhuǎn)座子序列進(jìn)行同源搜索，每種方法都有各自的優(yōu)勢和缺陷［24］。本文首先綜合使用前兩類算法的軟件挖掘了海島棉基因組中的轉(zhuǎn)座子序列，再進(jìn)一步將這些序列合并、去冗余，構(gòu)建出一個(gè)海島棉特異轉(zhuǎn)座子序列數(shù)據(jù)庫。最后依據(jù)該庫通過RepeatMasker軟件用同源搜索的方法注釋海島棉基因組中的轉(zhuǎn)座子序列，從而獲得了海島棉基因組中非常完整的轉(zhuǎn)座子信息。進(jìn)一步的數(shù)據(jù)分析發(fā)現(xiàn)海島棉基因組中每Mb染色體包含的LTR反轉(zhuǎn)錄轉(zhuǎn)座子的拷貝數(shù)與染色體的大小具有一定的負(fù)相關(guān)性，但這一特征在其他物種中是否存在還需要進(jìn)一步的研究。

LTR反轉(zhuǎn)錄轉(zhuǎn)座子在海島棉A亞組和D亞組的分布曲線有較大的差別（圖3），而在相同亞組內(nèi)部的各染色體上則具有類似的分布曲線（A亞組4號(hào)染色體除外）。高拷貝數(shù)的LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族在海島棉A亞組和D亞組的分布特征也具有一定的差異。此外，在轉(zhuǎn)座子活性研究中發(fā)現(xiàn)一個(gè)Copia類轉(zhuǎn)座子僅插入到海島棉A亞組的HD1基因中，而在D亞組的HD1基因中則沒有插入［25］。這些數(shù)據(jù)都表明轉(zhuǎn)座子在四倍體海島棉A亞組和D亞組中并不是完全相同的。然而，LTR反轉(zhuǎn)錄轉(zhuǎn)座子周邊基因的富集分析則表明A、D亞組之間非常類似。研究通過比較四倍體海島棉A、D亞組LTR反轉(zhuǎn)錄轉(zhuǎn)座子的數(shù)量與分布特征使我們對其有了更好的認(rèn)識(shí)，這將為海島棉基因組研究提供數(shù)據(jù)支持。

4 結(jié)論

本研究結(jié)果表明，絕大多數(shù)LTR反轉(zhuǎn)錄轉(zhuǎn)座子家族被海島棉A亞組和D亞組共同擁有，同時(shí)，兩個(gè)亞組也分別存在少數(shù)特異家族。海島棉染色體的大小與LTR反轉(zhuǎn)錄轉(zhuǎn)座子的數(shù)量有關(guān)。此外，本研究也發(fā)現(xiàn)在海島棉基因組中，Gypsy超家族分布較多的位置基因分布較少，但Copia超家族的分布則與基因沒有明顯的關(guān)系。

［1］Oliver KR, McComb JA, Greene WK. Transposable elements：powerful contributors to angiosperm evolution and diversity［J］.Genome Biol Evol, 2013, 5（10）：1886-1901.

［2］Bennetzen JL. Transposable element contributions to plant gene and genome evolution［J］. Plant Mol Biol, 2000, 42（1）：251-269.

［3］Vitte C, Panaud O. LTR retrotransposons and flowering plant genome size：emergence of the increase/decrease model［J］. Cytogenet Genome Res, 2005, 110（1-4）：91-107.

［4］Wicker T, Sabot F, Hua-Van A, et al. A unified classification system for eukaryotic transposable elements［J］. Nat Rev Genet, 2007, 8（12）：973-982.

［5］Kobayashi S, Goto-Yamamoto N, Hirochika H. Retrotransposoninduced mutations in grape skin color［J］. Science, 2004, 304（5673）：982.

［6］Mirouze M, Reinders J, Bucher E, et al. Selective epigenetic control of retrotransposition inArabidopsis［J］. Nature, 2009, 461（7262）：427-430.

［7］Domingues DS, Cruz GM, Metcalfe CJ, et al. Analysis of plant LTR-retrotransposons at the fine-scale family level reveals individual molecular patterns［J］. BMC Genomics, 2012, 13（1）：137.

［8］Liu X, Zhao B, Zheng HJ, et al.Gossypium barbadensegenome sequence provides insight into the evolution of extra-long staple fiber and specialized metabolites［J］. Sci Rep, 2015, 5：14139.

［9］Yuan D, Tang Z, Wang M, et al. The genome sequence of Sea-Island cotton（Gossypium barbadense）provides insights into the allopolyploidization and development of superior spinnable fibres［J］. Sci Rep, 2015, 5 ：17662.

［10］Wang K, Huang G, Zhu Y. Transposable elements play an important role during cotton genome evolution and fiber cell development［J］. Sci China Life Sci, 2016, 59（2）：112-121.

［11］Hu G, Hawkins JS, Grover CE, et al. The history and disposition of transposable elements in polyploidGossypium［J］. Genome,2010, 53（8）：599-607.

［12］Hawkins JS, Kim H, Nason JD, et al. Differential lineage-specific amplification of transposable elements is responsible for genome size variation inGossypium［J］. Genome Res, 2006, 16（10）：1252-1261.

［13］Xu Z, Liu J, Ni W, et al. GrTEdb：the first web-based database of transposable elements in cotton（Gossypium raimondii）［J］.Database（Oxford）, 2017, 2017（1）.

［14］McCarthy EM, McDonald JF. LTR_STRUC ：a novel search and identification program for LTR retrotransposons［J］.Bioinformatics, 2003, 19（3）：362-367.

［15］Ellinghaus D, Kurtz S, Willhoeft U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons［J］. BMC Bioinformatics, 2008, 9（1）：1-14.

［16］Edgar RC, Myers EW. PILER ：identification and classification of genomic repeats［J］. Bioinformatics, 2005, 21 Suppl 1 ：i152-i158.

［17］Huda A, Jordan IK. Analysis of transposable element sequences using CENSOR and RepeatMasker［J］. Methods Mol Biol, 2009,537（537）：323-336.

［18］Feschotte C, Keswani U, Ranganathan N, et al. Exploring repetitive DNA landscapes using REPCLASS, a tool that automates the classification of transposable elements in eukaryotic genomes［J］.Genome Biol Evol, 2009, 1（1）：205-220.

［19］Bao W, Kojima KK, Kohany O. Repbase Update, a database of repetitive elements in eukaryotic genomes［J］. Mob DNA, 2015,6（1）：11.

［20］Fu L, Niu B, Zhu Z, et al. CD-HIT：accelerated for clustering the next-generation sequencing data［J］. Bioinformatics, 2012, 28（23）：3150-3152.

［21］Tarailo-Graovac M, Chen N. Using RepeatMasker to identify repetitive elements in genomic sequences［J］. Curr Protoc Bioinformatics, 2009, Chapter 4：4-10.

［22］Naquin D, D’Aubenton-Carafa Y, Thermes C, et al. CIRCUS：a package for Circos display of structural genome variations from paired-end and mate-pair sequencing data［J］. BMC Bioinformatics, 2014, 15（1）：198.

［23］Ye J, Fang L, Zheng H, et al. WEGO ：a web tool for plotting GO annotations［J］. Nucleic Acids Res, 2006, 34（Web Server issue）：W293-W297.

［24］Lerat E. Identifying repeats and transposable elements in sequenced genomes：how to find your way through the dense forest of programs［J］. Heredity（Edinb）, 2010, 104（6）：520-533.

［25］Cao Y, Jiang Y, Ding M, et al. Molecular characterization of a transcriptionally active Ty1/copia-like retrotransposon inGossypium［J］. Plant Cell Rep, 2015, 34（6）：1037-1047.