Tc1/Mariner轉座子超家族在褶皺臂尾輪蟲中的進化和表達分析?

范正杰，王玉玨，張全啟，2，3??

(1.中國海洋大學海洋生物遺傳學與育種教育部重點實驗室，山東 青島 266003；2.青島海洋科學與技術試點國家實驗室 海洋漁業科學與食物產出過程功能實驗室，山東 青島 266237；3.中國海洋大學三亞海洋研究院 海南省熱帶水產種質重點實驗室，海南 三亞 572000)

轉座子(Transposons)，又被稱為轉座元件(Transposable elements, TEs)，是一段能通過非同源重組在宿主基因組內或基因組間進行移動的散在分布的重復序列[1]。真核生物的轉座子按轉座機制可分為兩大類：Ⅰ類元件即逆轉錄元件，是通過逆轉錄酶使用“復制-粘貼”的機制進行轉座的RNA轉座子；Ⅱ類元件DNA轉座子，則是使用基于DNA的“剪切-粘貼”轉座模式直接將轉錄子從原始位點整合至目標位點的轉座子，其中最大的亞綱是兩端有末端反向重復序列(Terminal inverted repeat，TIR)結構的TIR轉座子，包含著分布廣泛且種類眾多的超家族，例如Tc1/Mariner超家族。Tc1/Mariner轉座子超家族可能是自然界中分布最廣的轉座子超家族，在輪蟲、真菌、植物、魚類和哺乳動物等多種生物中都有其存在，然而其絕大多數是失活突變[2-4]，且它因Tc1和Mariner合稱而得名[5]。其中Tc1轉座子于秀麗隱桿線蟲(Caenorhabditiselegans)基因組研究時被發現[6]，Mariner是最小的自主DNA轉座子之一，以其水平轉移的傾向而聞名，目前已經發現其在動物、植物、真菌中均廣泛分布，即使在蛭形輪蟲這類具有高轉座抑制的物種中，Mariner轉座子也通過其強大的不依賴宿主的特性而大量的存在與繁殖[4]。

Tc1/Mariner轉座子的長度在1～5 kb之間，編碼282～345個氨基酸的轉座酶，其兩側有2個TIR，且TIR長度在17～1 100 bp之間[3,7]。來自不同Tc1/Mariner轉座子家族的轉座酶蛋白序列并不十分相似，但都具有兩個特征結構域：一個包含螺旋-轉角-螺旋(Helix-Turn-Helix, HTH)氨基末端，用于識別和結合TIR；另一個包含DDD/E羧基末端催化基序結構域，其第一和第二個天冬氨酸(D)殘基之間至少有92個氨基酸，而第二個天冬氨酸(D)殘基和第三個天冬氨酸(D)殘基或谷氨酸殘基(E)之間的距離是變量。目前Tc1/Mariner超家族大致可根據此數目分為7個主要亞家族，分別為DD34E(Tc1)、DD34D(mariner)、DDxD(pogo)、DD37D(maT)、DD37E、DD39D和DD41D(rosa)。

Tc1/Mariner轉座子不依賴宿主的特征來執行轉座過程，因此其存在不局限于一個特定的宿主。事實上，已經發現了許多Tc1/Mariner在不同宿主之間水平轉移的情況，例如在褐帶卷蛾中發現的一個Mariner轉座子通過水平轉移方式轉移到它的寄生體后，也在其它寄主和寄生體間發現了Mariner轉座水平轉移的現象，表明寄生—寄主關系可能在轉座子水平轉移中占重要地位[8]。該現象存在于海洋甲殼類動物之間[9]、不同目昆蟲之間[10-11]，甚至存在于不同門的生物之間，如人類和寄生線蟲之間的轉移[12]。然而，目前還不清楚轉座子是如何侵入新的基因組的。參與這種水平轉移的潛在載體是外部寄生者，如螨蟲(可能是果蠅中P元件水平轉移的載體[13])？還是內部寄生者，如病毒[14]？目前針對Tc1/Mariner的水平轉移在昆蟲中報道較多[15]。

轉座子在基因組中的分布并非是隨機的，部分轉座子與一些功能元件有密切的聯系[16-17]，轉座子不僅影響了基因組的結構，還參與了基因表達調控的過程，例如在人類和小鼠中發現Alu插入的數量與基因差異表達存在相關性[18]。在大鼠和小鼠中發現，長末端重復序列(Long terminal repeated，LTR)和長散在重復序列(Long interspersed nuclear elements，LINE)與基因表達正相關[19]。此外，轉座子插入到基因或基因側翼區域也可能帶來基因功能突變[20]。

褶皺臂尾輪蟲(Brachionusplicatilis)在分類學上隸屬于輪蟲動物門(Rotifera)，單巢綱(Monogononta)、游泳目(Ploimida)、臂尾輪蟲科(Brachionidae)、臂尾輪蟲屬(Brachionus)，是目前唯一能在海水養殖中實現規模化養殖的輪蟲，具有抗逆性強、營養豐富、繁殖快和游動慢等特點，因此該物種成為海水魚蝦蟹類幼體重要的開口餌料[21]。同時褶皺臂尾輪蟲具有典型的孤雌生殖世代和有性生殖世代交替的兼性繁殖生活史，會在外界環境變化，如溫度、營養、種群密度和pH溶解氧劇變時進入混交世代。先前的轉座子對于輪蟲門基因組影響的研究大多關注含量及結構[22]，而本研究基于褶皺臂尾輪蟲基因組數據，針對褶皺臂尾輪蟲基因附近的轉座子家族進行鑒定和富集，并針對其Tc1/Mariner超家族進行了進化和表達分析，以從功能的視角去探索Tc1/Mariner超家族在褶皺臂尾輪蟲基因組中扮演的角色。

1 材料與方法

1.1 褶皺臂尾輪蟲基因組重復序列注釋

褶皺臂尾輪蟲基因組數據(Bioproject:PRJNA 719948)和基因注釋信息由本實驗室分析獲得。重復序列注釋方法分為同源序列比對和從頭預測兩類。同源序列比對方法基于RepBase[23]，使用RepeatMasker[24]和Repeatproteinmask[25]軟件識別與已知重復序列相似的序列。從頭預測使用LTR_FINDER[26]，Piler[27]，RepeatScout[28]，RepeatModeler[29]等軟件：首先建立從頭測序重復序列庫；再通過Repeatmasker 軟件預測。此外，使用TRF[30](http://tandem.bu.edu/trf/trf.html)尋找基因組中串聯重復序列(Tandem repeat)。重復序列的注釋均使用軟件默認參數。將RepeatModeler、RepeatScout、Piler、LTR_finder 軟件預測出來的結果結合RepBase 庫采用Uclust[31]的軟件(遵從80-80-80 原則)進行整合獲得最終注釋結果。

1.2 轉座子及鄰近基因表達分析

褶皺臂尾輪蟲轉錄組數據(Bioproject: PRJNA 720972)由本實驗室測序分析獲得，取自8個發育階段：非混交胚胎、非混交雌性、雄性、雄性胚胎、攜帶雄性胚胎的混交雌性、攜帶休眠卵的混交雌性、新產休眠卵、休眠三個月后的休眠卵(見圖1)。使用Genomic-Ranges軟件[32]及rtracklayer軟件[33]對轉座子及基因結構進行了統計，并且以基因兩側的5 000 bp的長度分別作為基因的上游和下游部分。轉錄組數據進行reads過濾，并使用STAR軟件[34]進行比對得到SAM文件，再使用Samtools軟件[35]轉換為按名字排序的BAM文件并用于后續的分析。將基因注釋、轉座子注釋文件和各個時期轉錄組進行比對結果，使用TEtranscripts軟件[36]進行counts統計，并使用R語言軟件中的DEseq2軟件[37]進行標準化及處理，使用dplyr、factoextra[38]、FactoMineR軟件[39]進行平行組及各個組之間的主成分分析以確保質量，最后用TBtools軟件[40]來得到各個發育時期的所有轉座子及基因的FPKM表達量。規定只有在各個發育時期FPKM≥3時認為該轉座子有表達，進行篩選后將有表達的轉座子使用R語言軟件中的clusterProfiler軟件[41]進行富集分析以獲得轉座子家族，對富集到的轉座子家族與臨近基因的各世代表達進行Person相關性分析，以絕對值0.5作為相關的系數標準，使用R語言軟件中的ggplot2[38]進行繪圖。

圖1 褶皺臂尾輪蟲生活史模式圖

1.3 Tc1/Mariner轉座子超家族的表達和進化分析

利用shell腳本提取褶皺臂尾輪蟲基因組中的Tc1/Mariner轉座子超家族注釋信息，并分別統計其家族種類、拷貝數和表達量，最后對臨近基因的功能進行分析。

為了獲得基于Tc1/Mariner轉座子完整轉座酶，從NCBI及Repbase中下載了544條Tc1/Mariner轉座酶序列，以褶皺臂尾輪蟲基因組為目標庫，使用tBlastN尋找轉座酶序列(相似度>30%、長度>80、e值<10-5)，對每個拷貝的兩側分別延長1 000 bp，使用Seqkit[42]的subseq功能基于bed格式提取完整fasta序列，使用TBtools中的ORF_Prediction獲取序列的ORF序列，并保留每個序列預測出的長度大于300 bp的氨基酸序列，在MuscleX軟件中使用對齊(即使用默認設置)，手動進行篩選，保留有完整的DDE/D結構域的轉座酶序列，并基于DDE/D結構域確認其家族，在PHYRE2網站(http://www.sbg.bio.ic.ac.uk/phyre2)上預測各家族轉座酶的空間結構。

使用已知的具有完整DDE結構域的轉座酶家族(來自于果蠅、蛭形輪蟲、家蠶等)作為參考，使用muscle軟件[43]進行對齊，用ProtTest軟件[44]選擇最優氨基酸代替模型，基于該模型使用RAxML軟件[45]構建系統進化樹,bootstrap設為1 000。對于獲得的系統樹，使用EvolView軟件[46]進行美化。

利用Tc1/Mariner各個亞家族的注釋信息，使用TEtranscripts軟件處理結果來得到各個時期的表達量。用GenomicRanges軟件及rtracklayer軟件尋找各家族附近基因，并基于基因功能注釋得到Tc1/Mariner各亞家族附近基因的功能富集結果。

2 結果

2.1 褶皺臂尾輪蟲基因附近轉座子家族富集分析

主成分分析結果顯示(見圖2)，8個組24個轉錄組聚類情況良好，可用于后續分析。褶皺臂尾輪蟲基因組中共鑒定到419 244個轉座子(不含簡單重復序列)，從中篩選得到了92 225個轉座子在褶皺臂尾輪蟲的生命周期中的不同階段發生了表達。共富集到10個轉座子家族中，其中Ⅰ轉座子家族有3個：ERV4、Pao和tRNA-Deu；Ⅱ轉座子家族有7個：包括隸屬于Tc1/Mariner超家族的Mariner、Fot1和Tigger3個亞家族，以及不屬于Tc1/Mariner超家族的hAT1、hAT、Academ和Sola。將這10個轉座子家族的轉座子拷貝表達量與鄰近基因的表達量進行相關性分析，發現多數轉座子的表達模式與附近基因的表達模式呈正相關，負相關者較少，其中呈相關性的Mariner家族的轉座子以分布在基因下游為主(見圖3a)，而Fot1家族的轉座子以mRNA分布居多(見圖3b),Tigger家族的轉座子則分布較為均勻(見圖3c)。

(AE: 非混交胚胎Amictic embryos；AMF: 非混交雌性Amictic females；M: 雄性Males；ME: 雄性胚胎Male embryos；M: 攜帶雄性胚胎的混交雌性Mictic females with male embryos；MFRE: 攜帶休眠卵的混交雌性Mictic females with resting eggs；RE: 新產的休眠卵Resting eggs；RED: 休眠三個月后的休眠卵Resting eggs after three months’ dormancy)

((a)Mariner家族Mariner family；(b)Fot1家族Fot1 family；(c)Tigger家族Tigger family。第一列為Pearson相關性系數，第二列-log10(p值)，第三列轉座子位置，第四列為轉座子FPKM表達量，第五列為轉座子附近基因FPKM表達量。1st column: correlation; 2nd column:-log10(p-value);3rd column: TE position; 4th column: TE FPKM expression level; 5th column: nearby genes’ FPKM expression level.AE: 非混交胚胎Amictic embryos；AMF: 非混交雌性Amictic females；M: 雄性Males；ME: 雄性胚胎Male embryos；M: 攜帶雄性胚胎的混交雌性Mictic females with male embryos；MFRE: 攜帶休眠卵的混交雌性Mictic females with resting eggs；RE: 新產的休眠卵Resting eggs；RED: 休眠三個月后的休眠卵Resting eggs after three months’ dormancy)

2.2 Tc1/Mariner的鑒定和轉座酶結構及系統發育分析

褶皺臂尾輪蟲中共發現了23 378個Tc1/Mariner轉座子拷貝，并分為7個亞家族，其中拷貝數目較多的3個家族為Tc1、Pogo和Mariner家族(見表1)。

統計了7個亞家族在基因內部、基因上下游及基因間區的分布情況(見表1)。所有亞家族表現出了相似的分布特征，即均在基因間區分布最多。除了Trigger亞家族(48.19%)和Tc1亞家族(47.73%)，其余5個亞家族在基因間區的轉座子分布占比都超過了50%。Sagan亞家族在基因間區的分布占比最高(65.69%)。在基因內部分布的情況中可以發現，相較于其他6個亞家族，Sagan亞家族在基因內部的分布最少(僅占2.93%)。

表1 Tc1/Mariner轉座子7個亞家族的位置分布及拷貝數

基于同源比對的方法，本文作者在褶皺臂尾輪蟲基因組中共發現了29條完整的Tc1/Mariner轉座子序列。為了確定褶皺臂尾輪蟲中Tc1/Mariner轉座子的系統發育關系，基于轉座酶序列使用Repbase和NCBI中下載的近源物種中的Tc1/Mariner轉座子基因序列進行了系統發育分析。結果發現褶皺臂尾輪蟲中的Tc1/Mariner轉座子可以分為3個不同的亞家族即Tc1、Tc2和Pogo，其中大多數為Pogo轉座子(見圖4)。系統發育分析結果表明，褶皺臂尾輪蟲Tc1亞家族與家蠶的Tc1亞家族關系較近，而與果蠅中的Tc1家族關系較遠；Tc2亞家族與秀麗隱桿線蟲的Tc2聚為姐妹群；Pogo亞家族與真渦蟲Pogo亞家族關系較近。總體來看，這與褶臂尾輪蟲Tc1/Mariner轉座子系統發育關系與物種之間的親緣關系基本一致。

(不同的背景顏色代表不同的Tc1/Mariner轉座子亞家族，紅色分類名代表褶皺臂尾輪蟲中鑒定到的序列。Clade colors denote different subfamilies of Tc1/Mariner transposons; Taxon names in red indicate sequences from B. plicatilis.)

然后對29條完整轉座酶進行了結構預測和序列比對。Tc1/Mariner轉座酶的N端為DNA結合區域，C端為催化區域，空間結構預測結果(見圖5)表明褶皺臂尾輪蟲的Tc1/Mariner轉座酶均在Phyre2中高置信命中了c3hosA(轉座酶)結構，其N端NA結合區域含有兩段螺旋-轉角-螺旋(HTH，helix-turn-helix)結構，而在轉座酶的C端均具有DDE/D基序作為DNA剪切轉座活性的催化區域。序列對比結果顯示Tc1/Mariner超家族轉座酶催化區域的第二個天冬氨酸(D)和第三個天冬氨酸(D)之間，或第二個天冬氨酸(D)和谷氨酸(E)之間有30～35個氨基酸，基于其數目不同可以分為不同的家族，分別是Tc1家族(DD34E，見圖6a)、Tc2(DD35D，見圖6b)家族和Pogo(DDxD，見圖6c)家族。

圖5 Tc1/Mariner轉座酶的空間結構預測

2.3 Tc1/Mariner轉座子超家族的表達模式及功能

Tc1/Mariner轉座子7個亞家族在各個發育階段均有表達，并且具有相似的表達模式，但表達量并非恒定一致(見圖7)，在孤雌生殖階段的非混交胚胎和非混交雌性中表達量偏低，在有性生殖階段中，表達量較高，其中轉座子在雄性階段表達量最高，在攜帶雄性胚胎的需精雌性和攜帶休眠卵的混交雌性中表達量較低，在休眠卵(包括新產和休眠三個月)中的表達量較高。

(AE: 非混交胚胎Amictic embryos；AMF: 非混交雌性Amictic females；M: 雄性Males；ME: 雄性胚胎Male embryos；MFM: 攜帶雄性胚胎的混交雌性Mictic females with male embryos；MFRE: 攜帶休眠卵的混交雌性Mictic females with resting eggs；RE: 新產的休眠卵Resting eggs；RED: 休眠三個月后的休眠卵Resting eggs after three months’ dormancy.)

為了了解各個亞家族臨近基因的功能，本文作者對每個亞家族臨近的基因進行GO和KEGG富集分析。

GO功能富集(見圖8)發現了大量離子相關功能，包括離子結合和離子轉運。其中離子結合功能包括鋅離子結合和鎂離子結合等，而離子轉運涉包括陰離子轉運等。同時也發現了細胞周期調控功能，包括細胞周期功能DNA修復功能和多細胞生物發育的調節功能。此外還發現應對刺激的功能，包括在應對外界壓力和生物過程的調控的功能，以及損傷的反應。

圖8 篩選出Tc1/Mariner臨近基因的GO功能富集

KEGG通路富集(見圖9)中發現的細胞周期通路和環境信息處理的通路有大量的富集。其中大量富集于細胞周期的通路包括減數分裂和衰老通路。大量富集于環境信息處理的通路包括鞘脂信號通路、Notch、Jak-STAT、HIF1、TGF-beta及Wnt信號通路等。

圖9 篩選出的Tc1/Mariner臨近基因的KEGG信號通路富集

總之，Tc1/Mariner轉座子超家族附近基因中有大量涉及離子相關、應對刺激、環境信息處理及生物發育調節的功能和通路。

3 討論

通過統計褶皺臂尾輪蟲基因組中的轉座子的分布，發現分布在基因內含子區域的轉座子較多。相較于基因編碼區域，內含子區域的轉座子受到的選擇壓力更小,從而獲得相較于其他區域更高的存活機會[47]，同時內含子區域的轉座子插入可能影響mRNA的可變剪切及轉錄產物的正確加工。但是目前還沒有足夠的證據證明褶皺臂尾輪蟲中轉座子的插入位置偏倚是隨機插入后受到選擇壓力的選擇結果，還是在轉座子本身的偏好性帶來的[48]。

轉座子在基因組中的分布與基因表達的各個方面相關，包括表達水平、轉錄本多樣性和基因調控因子的活性[49]。目前，關于轉座子家族作為整體對基因表達的影響的研究還較少，有研究報道LINE1家族在小鼠原始生殖干細胞和著床前胚胎中RNA含量豐富并且其敲除會抑制原始生殖干細胞更新[50]。本研究通過使用生物信息學手段分析了轉座子表達與輪蟲各個生命階段的相關關系，在發育的各個階段都鑒定到了表達的轉座子及富集的轉座子家族。共富集到了10個轉座子家族，其中DNA類型的轉座子含優勢，其中有Mariner、Fot1和Tigger3個亞家族均來自于Tc1/Mariner超家族。同時發現轉座子的表達與臨近基因的表達以正相關為主。

已有文獻報道脊椎動物中Tc1/Mariner轉座子具有多樣性，其豐度各有差異[51]。Tc1/Mariner轉座子超家族在哺乳動物和鳥類中表現出非常低的多樣性和豐度[52-53]，與哺乳類和鳥類中不同，早前的研究在硬骨魚的基因組中觀察到了Tc1/Mariner轉座子的高多樣性和高豐度[54]。例如轉座子可以插入到基因組中的新的調控元件、外顯子及內含子中介導基因融合和基因沉默[55]。而在Tc1/Mariner超家族中，許多轉座子例如Tc1亞家族不僅可以參與基因表達調控，還可以促進新基因的產生，進而使物種適應新環境[56]。本研究發現Tc1/Mariner轉座子在褶皺臂尾輪蟲基因組中表現出多樣性，共鑒定出7個Tc1/Mariner轉座子亞家族。根據其催化區域第二個天冬氨酸(D)和第三個天冬氨酸(D)之間，或第二個天冬氨酸(D)和谷氨酸(E)之間的氨基酸數目將其鑒定為Tc1(DD34E)、Tc2(DD35D)、Mariner(DD34D)、Pogo(DDxD)、Sagan(DD30D)、Tigger(DD32D，DD36D)和Fot1(DD30D)。然而，在自然界中發現的絕大多數Tc1/Mariner轉座子被認為是有缺陷的，因為它們在轉座酶基因的編碼區域內包含了移碼、插入/缺失和過早終止密碼子，因此很少有轉座子(如Passport和Thm3)被認為是功能活性元件[57-58]。在哺乳動物和鳥類中，所有Tc1/Mariner轉座子都是化石轉座子且結構有缺陷[52-53,59]。在兩棲動物中，Tc1/Mariner超家族的活動似乎也已經滅絕，盡管Tc1/Mariner超家族占青蛙基因組的很大比例(約5%)，但大多數元件都是古老的，并且沒有獨立的家族被鑒定出來[60]。爬行動物中Tc1/Mariner多樣性也很低，Tc1/Mariner超家族是蜥蜴DNA轉座子中第二多的超家族，但只有一個獨立的家族被鑒定出來[61]。在褶皺臂尾輪蟲中，Tc1/Mariner轉座子具有多樣性，這可能是因為在淡水生態系統中，洪水和干旱等頻繁的壓力可以加快轉座活動，從而有助于宿主產生新的變異以適應變化的環境[25, 62-63]。

Tc1/Mariner轉座酶中最后兩個殘基之間的氨基酸數量在30～36之間。除去缺陷拷貝中的突變，大多數DD34E(Tc1)和DD34D(Mariner)元件包含典型的34個氨基酸間距(DD34E和DD34D)[64]。在本研究中，本文作者發現褶皺臂尾輪蟲的DD34E(Tc1)中的大部分元件都表現出典型的DD34E 結構。在DDxD(Pogo)亞家族在不同物種中可以觀察到距離的變化，例如真菌的Flipper、Pot2和Fot1中的DD35D[65-66]、果蠅的DD30D[25,67]、人類的DD32D(Tigger1)[68]，而在硬骨魚中，所有完整的轉座子元件都顯示出非常保守的DD35D 結構域[69]。本研究鑒定出的3個亞家族的29條完整轉座酶序列也均具有良好的保守性，包括Tc1-DD34E、Tc2-DD35D和Pog-DDxD，這一結果表明，硬骨魚類中Tc1/Mariner轉座酶的活性可能有其特殊性。

對褶皺臂尾輪蟲的Tc1/Mariner轉座子的各個亞家族的各發育階段的表達情況和臨近基因功能進行了分析，發現各個亞家族的表達模式基本相似，均在雄性發育階段表達量最高，在新產休眠卵和休眠三個月后的休眠卵的表達量較高，這2個時期均為輪蟲發育狀態發生變化的特殊時期，這期間轉座子處于活躍狀態，這為后續研究轉座子在輪蟲發育中的作用提供了重要的材料和線索。在對其臨近基因功能分析中，各個亞家族的GO功能諸多涉及到多細胞生物發育的調節、離子通道、細胞周期、對刺激反應的調控、DNA修復等；KEGG分析也發現與環境信息處理、生物發育調節有關。這些都可能與環境適應有關，而在輪蟲基因組中，DNA修復和離子轉運屬于輪蟲適應環境比較獨特的功能，涉及到神經調節、滲透壓調節、細胞周期調節、無性世代有害變異修復和休眠卵滯育等一系列的生物學過程。顯然，輪蟲基因組對環境的適應是綜合的調節，不局限于基因層面，還包含了轉座子表達對基因表達的調控。因此，有關Tc1/Mariner轉座子具體的生物學功能，需要進一步的實驗去驗證。

本研究通過對褶皺臂尾輪蟲Tc1/Mariner轉座子超家族在基因組中的分布、進化和表達進行系統分析，這為系統地認知Tc1/Mariner轉座子并從功能角度理解其對褶皺臂尾輪蟲基因組的作用提供了新的線索。