無芒隱子草全基因組水平全長LTR 反轉錄轉座子鑒定及其中斷基因分析

王藝蒙，馬甜甜，歐陽子鳳，張吉宇

（草地農業生態系統國家重點實驗室，蘭州大學草地農業科技學院，甘肅 蘭州 730020）

無芒隱子草（Cleistogenes songorica）屬于禾本科隱子草屬，是荒漠草地旱生種，可以在年降水量100 mm 的西北內陸地區良好生長，對穩定荒漠草原生態系統具有重要作用［1］。其根系可以改善土壤團聚結構，增強土壤的凝聚力，在干旱半干旱地區具有抗水土流失能力，無芒隱子草還具有適口性好、耐踐踏、青綠期長等多種優良特性，是優良的飼用禾草和草坪草［2］。其為異源四倍體（2n=4X=40），“騰格里”無芒隱子草是2016 年由我國草品種審定委員會審定登記的野生栽培品種（登記號：499），抗旱性極強，在干旱半干旱地區的生態治理、園林綠化和草原恢復中具有非常重要的應用價值［3］。

轉座子（transposon）為植物基因組中含有的大量可移動的遺傳因子，根據轉座過程可分為DNA 轉座子和反轉錄轉座子兩類［4］。反轉錄轉座子（retrotransposon）是一類存在于植物基因組中的特殊的DNA 序列，以一種“復制粘貼”的轉座機制來實現快速的自我增殖，其最早在哺乳動物和酵母基因組中被發現［5］。

相比于其他轉座子，全長LTR 反轉錄轉座子（long terminal repeat retrotransposon，LTR-RTs）有其獨特的結構特征和轉座機制，最顯著的特征是在轉座子的5′和3′端出現一對同向高度同源序列，長度從幾百bp 到約5000 bp 不等，這兩段標志序列被稱為 TG-CA 盒（TG-CA box）［6］。在緊鄰 5′LTR 的下游區域，有一個和該物種中某類tRNA 序列的3′末端反向互補，長度為十多個bp 的序列，該序列是反轉錄的起始引物的結合位點，稱為引物結合位點 PBS（primer binding site）［7］。類似地，在 3′LTR 的上游附近，也有十多 bp 的序列負責引導 DNA 鏈的復制，稱該區域為PPT（poly purine trait）。完整的LTR-RTs 結構具有兩個LTR、LTR 兩側有一對二核苷酸回文基序，以及兩端的5 bp 靶位點重復TSD（target site domain）［8］。除此之外，根據pol 基因區的不同分為Copia超家族和Gypsy超家族，在Copia超家族中，反轉錄酶位點在整合酶位點的前面，而在Gypsy超家族中則剛好相反［9］。

研究表明，許多活性轉座子插入了宿主基因并中斷其功能，因此活性轉座子可以應用于鑒定功能基因的研究，例如在水稻（Oryza sativa）基因組中發現的Tos17、RIRE7，激活后的LTR-RTs 還可以作為植物轉基因技術中的基因轉移載體［10］，例如煙草（Nicotiana tabacum）中發現的Tto1，所以近年來的研究均關注具有潛在活性的全長LTR-RTs［11］。

本課題組于近期完成“騰格里”無芒隱子草基因組測序，基因組大小為543 M［12］，本研究基于無芒隱子草的全基因組水平，篩選了潛在活性的全長LTR-RTs，對這些序列進行了系統進化樹的構建和插入時間的分析，對全長LTR-RTs 及被中斷的基因的位置關系進行分析，對被中斷基因進行GO 功能注釋和表達量分析。這有助于了解無芒隱子草的進化過程，為全長LTR 反轉錄轉座子在牧草中的后續研究奠定了基礎［13］。

1 材料與方法

1.1 全長LTR-RTs 篩選

使用 Repeat Masker 和 Repeat Protein Mask 軟件（http：//www.repeatmasker.org/cgi-bin/WEB Repeat Mask?er）基于 RepBase 庫（http：//www.girinst.org/repbase）［14］分別注釋基因組得到的部分候選 LTR-RTs。使用 LTRFinder［15］和 Repeat Modeler 對基因組中的重復序列進行De novo預測［16］，再以此為庫利用軟件 Repeat Masker 得到另一部分候選LTR-RTs。整合以上結果去冗余即為從無芒隱子草全基因組中預測得到的全部LTR-RTs。使用 RepClass 將獲得的 LTR-RTs 分類為Gypsy、Copia和 Other。

完整的LTR-RTs 長度在85～5000 bp，包含兩端LTR 區域，以及Gypsy和Copia超家族元件典型結構域中3個結構域的序列，結構域包括Gag（種屬特異抗原）、RT（反轉錄酶Reverse transcriptase）、RH（RNA 酶RNase H）、Int（整合酶 Integrase）、PR（編碼蛋白酶 Protease）［17］?；谏鲜鲈瓌t使用 LTR-FINDER［14］和 LTR-retriever［18］分兩步對無芒隱子草全基因組進行LTR-RTs 篩選，得到候選全長LTR-RTs，使用BLAST（the basic local align?ment search tool）與候選全長LTR-RTs 的兩個LTR 之間的編碼區域進行比對來確定結構域。

1.2 全長LTR-RTs 插入時間分析

通過LTR 反轉錄轉座子兩端LTR 序列同源性與宿主的堿基替換速率估算插入時間，當LTR 反轉錄轉座子進行轉座時，LTR 序列由同一模板DNA 合成，因此起初兩端LTR 序列同源性為100%［19］，但是隨著植物基因組的變化，LTR 序列會發生不同程度的突變，導致兩端LTR 序列同源性下降，所以根據LTR-RTs 兩端LTR 序列同源性估算插入時間。使用MEGA 7 比對LTR-RTs 的兩端LTR 序列獲得序列分化度（K），通過公式T=K/2r估測插入時間，r為 LTR 序列的平均替換率，K為序列分化度，本研究中使用 1.3×10?8替代率·bp?1·a?1［20］。

1.3 全長LTR-RTs 進化分析

從 Gypsy Database（http：//www.gydb.org/index.php/Main_Page）中下載已報道的Gypsy和Copia超家族序列的RT 氨基酸序列［21］（表1），與使用NCBI 的保守結構域數據庫預測全長LTR-RTs 的RT 氨基酸序列進行比對，利用 clusterX 和 FastTree 聚類和構建系統進化樹，使用在線軟件 iTOL（https：//itol.embl.de/）可視化［22］。

表1 用作系統發育分析的參考全長LTR-RTs 序列Table 1 Reference full-length LTR-RTs sequences for phylogenetic analysis

1.4 全長LTR-RTs 中斷基因功能注釋

基于無芒隱子草基因和全長LTR-RTs 的染色體起止位置，統計被中斷的基因。對用GO（gene ontology）進行注釋的不同數據進行比對，并以直觀的柱狀圖形式呈現數據結果的差異，序列比對篩選被中斷基因的同源基因，進行基因功能注釋，找尋可能具有功能的全長LTR-RTs［23］。

1.5 全長LTR-RTs 中斷基因的表達量分析

對無芒隱子草進行干旱脅迫處理，取對照、輕度脅迫、中度脅迫、重度脅迫、復水處理的地上地下部樣品進行RNA-seq 測序，參考Yan 等［24］的測序數據，對被中斷基因進行表達量分析。

1.6 數據統計與分析

使用 LTR-retriever 等進行 LTR-RTs 篩選，用 RepClass 將獲得的 LTR-RTs 分類，利用 clusterX 和 Fasttree 聚類和構建系統進化樹，用百邁克云平臺進行可視化。通過在線數據分析平臺Omic Share Tools 分析被中斷基因在干旱脅迫下的表達量。

2 結果與分析

2.1 無芒隱子草全基因組LTR-RTs

從無芒隱子草全基因組中共鑒定出了299079 個LTR-RTs。其中199460 個被分類為Gypsy超家族，81865 個被分類為Copia超家族，剩余17754 個被統稱為Other。這些LTR-RTs 長度為0.001 Kb 到44.72 Kb 不等，平均長度為0.78 Kb?？傞L度為144136979 bp，占無芒隱子草全基因組大小的26.54%。

用LTR-FINDER 對無芒隱子草基因組進行預測，共得到4689 個候選全長LTR 反轉錄轉座子，進一步用LTR-retriever 篩選得到1170 個高質量候選全長LTR-RTs。將這些反轉錄轉座子與數據庫進行對比后，產生了845 個結構全長 LTR-RTs，發現 410 條屬于Copia類型，435 條屬于Gypsy類型。從預測結果來看，Copia類型反轉錄轉座子占所有預測反轉錄轉座子的48.52%，略低于Gypsy類型的51.48%。有794 個全長LTR-RTs 可以被成功定位到染色體上，20 條染色體上的Copia∶Gypsy比率為0.65～2.50（表2）。盡管每條染色體上的全長LTRRTs 數量不等，但每千萬堿基LTR-RTs 密度沒有明顯差異。

表2 鑒定出的全長LTR-RTs 在無芒隱子草基因組20 條染色體上的分布情況Table 2 Distribution of the identified full-length LTR-RTs on the 20 chromosomes of the C.songorica genome

2.2 全長LTR-RTs 的分類

對預測出的435 條Copia超家族成員RT 結構域氨基酸序列與已鑒別的18 條其他物種的RT 結構域氨基酸序列進行聚類分析（圖1）。并對預測出的408 條Gypsy超家族和已鑒別的27 條其他物種的RT 結構域氨基酸序列進行聚類分析（圖2），參考Llorens 等［25］關于LTR-RTs 系統發育的分析（表1），聚類結果顯示無芒隱子草中存在Copia超家族中的Sire、Oryco、Retrofit、Tork4 個支系。在本研究中，Oryco是Copia超家族中最重要的支系，有138 個成員，最小的支系是Sire有39 個成員，另有104 條序列未分類。Gypsy超家族的RT 序列的6 個支系分別為Athila、Tat、Reina、CRM、Galadriel、Del/Tekay，而聚類結果顯示無芒隱子草中存在 5 個支系，僅Galadriel支系沒有成員。其中最大的支系是Tat有182 個成員，最小的支系是Athila僅有1 個成員，另有105 條序列未分類。對于構建的系統進化樹，雖然Copia和Gypsy兩大超家族的是以基因中和功能域的順序來區分，但是從序列相似度上也可以區分兩大家族［26］。

圖1 無芒隱子草全長LTR-RTs Copia 超家族逆轉錄酶（RT）氨基酸序列聚類樹Fig. 1 Phylogenetic tree of the reverse transcriptase（RT）in Copia retrotransposons of C. songorica

圖2 無芒隱子草全長LTR-RTs Gypsy 超家族逆轉錄酶（RT）氨基酸序列聚類樹Fig. 2 Phylogenetic tree of the reverse transcriptase（RT）in Gypsy retrotransposons of C. songorica

2.3 全長LTR-RTs 插入時間分析

對無芒隱子草全長LTR-RTs 篩選各階段的序列進行了插入時間的估測（圖3A），并對Gypsy和Copia超家族成員的插入時間進行比較分析（圖3B）。發現大部分LTR-RTs 的插入時間發生在9 百萬年間，而LTR-RTs 的轉座時間在1.0～1.5 百萬年達到最高值［27］。在Copia和Gypsy超家族中，插入時間超過5 萬年的反轉錄轉座子的比例都很小。

圖3 無芒隱子草全長LTR-RTs 插入時間分析Fig.3 Analysis of insertion time of full-length LTR-RTs in C. songorica

對每個超家族支系的插入時間的評估也顯示了全長LTR-RTs 的差異（表3）。在無芒隱子草基因組中鑒定的成員數大于10 的8 個支系和2 個未分類組中，只有Gypsy超家族中的支系CRM和未分類支系在0.5 百萬年間出現插入高峰，其他支系在更古老的時代表現出激增，例如Copia超家族中大部分的Retrofit在0.5～1.5 百萬年間插入（68%），大部分的Tork支系在2～4 百萬年間插入（75%），大部分的Oryco支系在0.5～2.5 百萬年間插入（80%）。

表3 無芒隱子草中Copia 和Gypsy 超家族成員在不同支系中的插入時間估計Table 3 Estimated insertion times of different lineages of full-length Copia and Gypsy superfamily members in C.songorica

由于LTR-RTs 插入基因序列中可能導致基因突變，造成功能缺陷或改變，并且這種突變是可遺傳的，因此對全長LTR-RTs 及被中斷的基因的位置關系進行分析。無芒隱子草中被全長LTR-RTs 中斷的基因有183 個，全長LTR-RTs 和基因的位置關系初步歸類如圖4，兩者之間的位置關系分為全長LTR-RTs 與被中斷基因在LTR區間產生交互（圖4A）、全長LTR-RTs 與中斷基因在編碼區產生交互（圖4B）、全長LTR-RTs 位于被中斷基因的內含子區間內（圖4C）、全長LTR-RTs 位于被中斷基因的外顯子區間內（圖4D）、被中斷基因位于全長LTR-RTs的LTR 區間內（圖4E）、被中斷基因位于全長LTR-RTs 的編碼區內（圖4F）以及全長LTR-RTs 與被中斷基因沒有交互等類型。

圖4 全長LTR-RTs 及被全長LTR-RTs 中斷的基因的關系Fig.4 Relationship between full-length LTR-RTs and genes interrupted by full-length LTR-RTs

無芒隱子草中被全長LTR-RTs 中斷的基因為183 個，GO 功能注釋到的被中斷基因有145 個。被中斷的基因在生物學功能中主要行使新陳代謝功能（metabolic process）、細胞學功能（cellular process）和單一有機體功能（single-organism process）。在細胞組分中主要行使細胞部分（cell part）、細胞器（organelle）和細胞膜部分（mem?brane）等。在分子功能方面，主要行使結合（binding）和催化作用（catalytic activity）等功能（圖5）。

圖5 被LTR-RTs 中斷基因的GO 注釋Fig.5 GO functional of genes interrupted by LTR-RTs in C. songorica

被中斷基因的KEGG 富集分析表明，基因顯著富集在剪切復合體（spliceosome）、內吞作用（endocytosis）、內質網中的蛋白質合成（protein processing in endoplasmic reticulum）以及植物激素信號傳導（plant hormone signal transduction）。結果表明，被LTR-RTs 中斷的基因主要參與蛋白質翻譯、化合物的代謝與修飾、信號傳導等途徑（圖6）。

圖6 無芒隱子草中被LTR-RTs 中斷基因的KEGG 富集分析Fig.6 KEGG enrichment analysis of genes interrupted by LTR-RTs in C. songorica

2.4 全長LTR-RTs 中斷基因的表達量分析

對不同干旱脅迫處理下的被全長LTR-RTs 中斷的183 個基因進行表達量分析（圖7），發現大部分基因都不表達，從中篩選出表達量較高的10 個被中斷基因（表4），這10 個被中斷的基因均為全長LTR-RTs 位于被中斷基因內的類型，其中一個是全長LTR-RTs 位于被中斷基因的外顯子區間內，其余都是全長LTR-RTs 位于內含子區間內。篩選出的基因都能被GO 功能注釋到，且主要行使激活脫落酸的信號通路、光周期、刺激應答和氧化還原酶活性等功能。

圖7 對無芒隱子草進行干旱脅迫后被中斷基因表達量熱圖Fig. 7 Gene expression heatmap of interrupted gene after drought stress in C. songorica

表4 無芒隱子草中10 個被中斷基因的GO 功能和在干旱脅迫處理下的表達量Table 4 GO functions of 10 interrupted genes and expression under drought stress and rehydration in C.songorica（FPKM）

3 討論

全長LTR 反轉錄轉座子有其獨特的結構特征和轉座機制，最顯著的特征是在轉座子的5′和3′端出現一對同向高度同源序列［28］。通過 LTR-FINDER 和 LTR retriever 分兩步對無芒隱子草全基因組進行LTRRTs 篩選，LTR-FINDER 能夠精確而迅速地在全基因組規模的序列中對LTR-RTs 作出準確預測，在很大程度上解決了計算時間問題。但由于LTR-FIND?ER 是從結構出發的，也可能帶來許多的假陽性結果［29］，所 以 進 一 步 用 LTR retriever 篩 選 候 選 全 長LTR-RTs。再通過NCBI 的BLAST 數據庫比對確認結構域，僅有845 個反轉錄轉座子包含兩端LTR 區域，以及Gypsy和Copia超家族元件典型結構域（Gag、RT、RH、Int、PR）中 3 個結構域的序列，而非全長LTR-RTs 共有29079 個，說明無芒隱子草中大量LTR-RTs 在進化過程中丟失了內部元件，失去活性［30］。

鑒于LTR-RTs 的進化速度很高，特別是當分析中包括不完整的LTR-RTs 時，對全長LTR-RTs 的界定是一個巨大的挑戰［31］。在本試驗中鑒定產生了845個全長 LTR-RTs，發現 410 條屬于Copia類型，435 條屬于Gypsy類型，Copia和Gypsy元件的平均比率為1.17。來自同一家族的兩個超家族之間的全長元件在無芒隱子草中分布相似［32］。為了闡明這些序列的系統發育關系，使用序列中的RT 結構域的氨基酸序列構建系統進化樹［33］。本研究參考Llorens 等［25］關于LTR-RTs 系統發育的分析，根據RT 序列的進化樹分析，將無芒隱子草基因組中發現的845 個全長LTR-RTs 分為不同的支系，這些支系也廣泛存在于水稻、小麥（Triticum aesti?vum）和擬南芥中［33］。除此之外，LTR-RTs 的數量和類型在無芒隱子草20 對染色體上分布相對均勻，沒有出現某一超家族或在某條染色體上的過度擴張和消除現象，這與小麥基因組中的LTR-RTs 的分布類似［34］。

LTR-RTs 之間的核苷酸同一性有助于估計其插入時間，因為單個LTR 反轉錄轉座子的兩個LTR 在整合后通常在核苷酸序列水平上是相同的［35］。事實上，學者們普遍認為，最近插入的LTR-RTs 通常在結構上更完整，而插入時間較長的包含更高比例的截短LTR-RTs 和單獨的LTR［35］。因此插入時間越近，具有潛在活性的可能性越大。在這項研究中發現兩個LTR 之間有很高的相似性，因此具有潛在的活性。從LTR-RTs 插入時間可知，無芒隱子草基因組在近期發生大量LTR-RTs 的插入，插入時間在1.0～1.5 百萬年達到最高值，且超過75%的全長序列都是在 3.0 百萬年內插入的［36］。

對禾本科Copia家族的LTR 反轉錄轉座子系統發育關系的研究表明［37］，Copia超家族起源非常古老，用60%蛋白相似度定義RT 的家族，發現很多RT 家族都可以在禾本科內共有。同樣的，在Gypsy家族上的研究結果也表明［38］，Gypsy超家族的一些主要成員在單子葉和雙子葉植物都共有，這說明Gypsy超家族的起源也非常古老。在不同的超家族和支系中，已確定的LTR-RTs 的擴增時間有很大差異，例如無芒隱子草中最近50 年的Gypsy超家族的活性高于Copia超家族。在這個意義上，本研究觀察到不同支系之間轉座活性的時間模式不同，例如Ret?rofit/Ale，Tork/GMR和Tat近期最活躍。相反，盡管最近插入了幾個Del/Tekay支系的成員，但活動的高峰出現在1.0～2.5 百萬年前。這表明LTR-RTs 在無芒隱子草基因組中的插入是波浪式的，而不是連續的，在其他物種中也發現了 LTR-RTs 的類似行為［39］。

LTR-RTs 與基因組的大小變化和進化有著密切的聯系，LTR-RTs 插入基因序列中可能導致基因突變，造成功能缺陷或改變，并且這種突變是可遺傳的［39］。LTR-RTs 的插入可以致使基因失活，例如葡萄的紫色表皮控制基因VvmybA1 被命名為Gret1 的Gypsy超家族LTR-RTs 插入，導致基因轉錄失活，表皮變為白色［40］。本試驗中無芒隱子草中篩選出的被全長LTR-RTs 中斷的基因為183 個，被GO 功能注釋到的基因有145 個，未被注釋到的基因長度都較短，這可能是因為短序列可能很難通過BLAST 得到注釋信息。為了闡明LTR-RTs 與基因有著密切的聯系，對全長LTR-RTs 和基因的位置關系進行大致分類，對被LTR-RTs 中斷的基因進行了GO 注釋和富集分析，結果顯示被中斷的基因主要作用于轉錄因子的激活，這可能表明全長LTR-RTs 對插入的基因具有偏好性與選擇性［41］，在水稻中也觀察到類似的現象［42］。

植物基因組中大量的轉座子在轉座和轉錄上均是失活的，只有少數反轉錄轉座子在生物或非生物因素脅迫下被激活，逆境能激活LTR-RTs 的轉錄，其原因可能是LTR-RTs 對宿主基因組的一種適應，其機制是通過重組獲得類似于寄主相關逆境誘導基因的啟動子，從而使LTR-RTs 在逆境條件下發生轉座。有研究表明毛竹（Phyl?lostachys heterocyclaMitford cv.Pubescens）基因組中非編碼基因存在著多個具備轉座能力的轉座子，生物和非生物脅迫會影響LTR-RTs 的活性表達［43］。在本試驗中，經RNA-seq 測序發現，僅有少數被全長LTR-RTs 中斷的基因在干旱脅迫下差異性表達，推測是這些基因在干旱脅迫下相關防御機制發揮作用，相關的全長LTR-RTs 也開始活躍起來，從而使中斷基因的表達量上調。在差異性表達的中斷基因中篩選出被GO 功能注釋到的10 個基因，發現其主要行使激活脫落酸的信號通路、光周期、刺激應答和氧化還原酶活性等功能，推測是脅迫誘導相關基因的表達，從而使LTR-RTs 在逆境條件下發生轉座。不僅是干旱脅迫，在組培、外源激素、化學誘變劑和其他逆境脅迫等處理下的LTR-RTs 中也存在差異表達，說明影響LTR-RTs 表達活性的條件具有多樣性［44］。還有研究表明，在擬南芥、小麥、煙草和玉米等植物中已克隆出具有活性的反轉錄轉座子，這些LTR-RTs 在基因內部插入時會導致基因表達和基因表達產物結構發生變化［45］，進一步說明了全長LTR-RTs 對插入的基因具有偏好性與選擇性，且LTR-RTs 可能與植物的防御機制有關。

4 結論

綜上所述，本研究基于無芒隱子草全基因組序列，篩選出具有潛在活性的全長LTR 反轉錄轉座子845 個，其中有410 個屬于Gypsy超家族，435 個屬于Copia超家族。對這些基因組序列進行分析，得出全長LTR 反轉錄轉座子的插入高峰期為1.0～1.5 百萬年間，其大量轉座發生在4 百萬年內。對全長LTR-RTs 及被中斷的基因的位置關系進行分析，以及對被中斷基因進行GO 功能注釋，經RNA-seq 表達量分析發現，僅有少數被全長LTR-RTs中斷的基因在干旱脅迫下差異性表達。這為全長LTR 反轉錄轉座子與中斷基因在無芒隱子草及牧草中的后續研究奠定了基礎。