類烏齊牦牛全基因組RNA編輯位點預測及剖析

王嘉博, 舒 濤,柴志欣, 王吉坤, 王 會,武志娟, 唐 友,鐘金城*, 姬秋梅

(1.西南民族大學 青藏高原動物遺傳資源保護與利用四川省、教育部重點實驗室，成都 610041；2.吉林農(nóng)業(yè)科技學院，長春 132309； 3.大麥、牦牛種質(zhì)資源與遺傳改良國家重點實驗室，拉薩 850000)

牦牛是高原特有的大型動物，盡管其生存環(huán)境惡劣，但仍然為高原地區(qū)的人們提供著生活必需品以及主要的勞役[1]。類烏齊牦牛是西藏東部地區(qū)特有的牦牛資源，具有優(yōu)質(zhì)牛肉風味以及極強的高原適應性。其肉顏色鮮紅，彈性致密，具有獨特的腥膻味，口味濃郁，營養(yǎng)價值高[2]。類烏齊牦牛在海拔3 700～4 500 m之間的高原、高山乃至河谷地帶均有分布。隨著現(xiàn)代分子生物學的發(fā)展，牦牛遺傳資源的挖掘、高原生態(tài)和高寒地區(qū)遺傳環(huán)境互作等研究逐漸成為熱點[3-5]。從基因水平解析牦牛組織分化、特異性差異表達以及各種器官與高原環(huán)境互作等問題的研究逐漸深入[3,6-7]，而如何在多個組織中找到基因間表達的異同以及探索器官特異性RNA編輯產(chǎn)物的分類還未見報道。

RNA編輯位點(RNA editing site, REs)是自然界中廣泛存在的一種對關(guān)鍵基因表達修飾的位點[8]，這種現(xiàn)象在哺乳動物中尤為普遍，目前，已有的研究表明，人、羊、牛等動物均存在這種現(xiàn)象[9-12]。RNA編輯位點是“一因多效”的重要體現(xiàn)，其可以高效地利用同一基因在不同環(huán)境或不同組織中針對表達產(chǎn)物進行修飾，從而影響生物體對于環(huán)境的反射和適應[13-15]。牦牛不同組織及不同發(fā)育階段的表達產(chǎn)物具有十分明顯的差異[16]，這是由于牦牛生存在復雜多變的高原地區(qū)[17]。探索牦牛不同組織間的RNA編輯位點，并通過生物信息手段進行分類和歸納，可以幫助理解牦牛在高原地區(qū)RNA多種形態(tài)編輯的遺傳規(guī)律，并解析牦牛不同組織中特有的編輯方式和類型，為深入研究牦牛遺傳資源與高原環(huán)境互作提供節(jié)點基因。

鑒于牛胚胎發(fā)育過程中組織分化和發(fā)育的差異，結(jié)合類烏齊牦牛組織在高原適應性中的作用[18]，本研究以類烏齊牦牛大腦、小腦、臀部脂肪和肌肉組織為研究對象，通過生物信息學對組織間RNA編輯位點進行探測和分類，最終通過風險評估預測這些位點對于組織生長分化的作用以及牦牛對高原惡劣環(huán)境適應性的反饋能力。

1 材料與方法

1.1 樣品采集

本試驗以西藏自治區(qū)類烏齊市3頭4.5歲成年、健康的類烏齊母牦牛作為試驗材料，符合動物福利屠宰后，采集其大腦(DN)、小腦(XN)、臀部肌肉(TJ)和臀部脂肪(TZ)，經(jīng)DEPC水沖洗干凈后，迅速將其包裹在錫箔紙中并置于液氮中保存，以用于后續(xù)試驗和分析。

1.2 類烏齊牦牛RNA-seq測序

共計3個個體，4個組織樣的RNA文庫建立，用TRIzol?Plus RNA Purification Kit (Invitrogen, USA)從樣本中提取總RNA，經(jīng)過瓊脂糖凝膠電泳回收18～30 bp的片段。建庫之前使用Epicentre Ribo-zeroTM rRNA Removal Kit剔除rRNA[19](Epicentre,美國)。通過反轉(zhuǎn)錄獲得總RNA的cDNA序列，經(jīng)過片段化處理和PCR擴增，利用Illumina HiSeqTM 4000 測序平臺進行高通量測序。

1.3 RNA編輯位點預測

最終得到3頭類烏齊牦牛在大腦、小腦、臀部肌肉及脂肪總RNA的Fasta格式測序文件，通過去掉比對rRNA庫中的序列和必要測序質(zhì)量控制過程，其中質(zhì)控采用Fastp軟件依照默認的參數(shù)進行。過濾后的干凈數(shù)據(jù)與Ensemble上最新的牦牛參考基因組(Bosgru_v3.0)進行比對，生成比對后bam文件并進行序列排序。將排好序的bam文件使用SPRINT[20]和JACUSA(v1.3.0)[21]進行REs的預測，SPRINT利用編輯位點變異基因型與周邊基因型變異程度進行區(qū)段化處理，利用區(qū)段內(nèi)連鎖關(guān)系排除假陽性的位點，而JACUSA可以利用重復樣本進行多重比對，相對而言假陰性較低。其中，SPRINT的命令行sprint_from_bam函數(shù)包進行，獲得3個個體4個組織的REs估計位點，利用自行編輯的R語言腳本對3個個體之間相同組織REs進行篩選，獲得每個組織特異性REs位點；JACUSA采用Github上推薦的命令，對同一組織3個個體進行一次直接估計REs。這種多樣本重復驗證估計的方法是JACUSA特有的功能，以此降低假陽性的發(fā)生。兩種方法均使用表達序列大于9進行過濾(參考SPRINT軟件使用手冊和參考文獻[20]中的過濾標準)，并比較兩種方法在4個組織中預測的結(jié)果，取交集作為候選REs位點。

1.4 高風險變異位點的預測

結(jié)合牦牛參考基因組fasta文件和cDNA、ncRNA以及生成蛋白質(zhì)序列文件，利用snpEff[22]軟件自行建立牦牛3號參考基因組的注釋數(shù)據(jù)庫，代碼為java -jar snpEff.jar build -gtf22 -v BosGru3，對兩款軟件的預測REs進行注釋，分析單位點的變異在基因組中與已知基因的相互關(guān)系，所屬RNA類型以及編輯變異在后續(xù)功能是否存在風險，使用風險等級評估這些位點。風險等級分為High、Moderate、Low以及Modifier。

1.5 利用RDDs驗證REs

RNA-DNA different sites(RRDs)利用表達譜數(shù)據(jù)反轉(zhuǎn)錄的cDNA序列與基因組中DNA序列進行直接對比確定REs的真實性，其原理和操作步驟參考梁浩等[8]的研究。選擇風險等級為中級的兩個基因進行試驗驗證，通過Primer(5.0)軟件設(shè)計特異性引物(擎科生物公司)，反轉(zhuǎn)錄表達產(chǎn)物成cDNA， 采用一代測序手段進行cDNA和DNA的測序及序列比對。在DNAMAN(V6，2005 Lynnon)上實現(xiàn)序列的比對，比較兩個基因在4個組織中的cDNA和DNA基因序列間的差異。

2 結(jié) 果

2.1 測序質(zhì)量鑒定及數(shù)據(jù)比對

經(jīng)過對下機數(shù)據(jù)的質(zhì)量控制以及篩選，總計獲得將近50G的原始數(shù)據(jù)(raw data)。經(jīng)過fastp軟件過濾后的數(shù)據(jù)GC含量在50%左右，Q20數(shù)據(jù)在98%以上，Q30數(shù)據(jù)在95%以上。通過BWA軟件對比牦牛第三版參考基因組，獲得12個樣本組織的bam文件，每個樣本的比對率超過90%。說明樣本與參考群體高度相似，比對reads可以用來進行下游分析(表1)。

表1 過濾后數(shù)據(jù)與參考基因組比對結(jié)果Table 1 The mapping result of data after filtering with reference genome

2.2 兩種方法預測牦牛REs

通過SPRINT和JACUSA兩款REs軟件對類烏齊牦牛腦組織和肌肉脂肪組織中REs進行預測。SPRINT在大腦、小腦、臀肌、臀脂中分別報告了28 864、 17 523、1 653和52 474個估計REs；JACUSA在大腦、小腦、臀肌、臀脂中分別報告了91 222、83 242、 61 293和171 630個估計REs(圖1)。將這些REs整合統(tǒng)計發(fā)現(xiàn)共計24 784 個編輯事件，在4個組織中均發(fā)現(xiàn)的編輯事件為4 015個。JACUSA預測的REs數(shù)量明顯比SPRINT多，這可能是由于JACUSA采用3個個體進行統(tǒng)一分析，只要有一個樣本出現(xiàn)REs, JACUSA就會報告這個位點。本研究目標是獲得組織特異性REs，因此選取兩款REs預測的交集部分并去除在所有組織中均檢測出來的REs(排除基因組雜合子的干擾)作為候選REs，按照REs富集在各染色體中的數(shù)量分布(圖2)，在臀脂中檢測出的REs最多，其次是大腦組織，然后是小腦組織，最后在臀肌組織中檢測出的REs最少。

A、B、C、D. 大腦、臀部肌肉、臀部脂肪和小腦組織中兩款軟件REs預測位點的共有交集文氏圖。圖中圓圈的大小和所包含REs數(shù)量成正比A, B, C, D. Indicates intersection Venn plot of REs predicted by two softwares in the DN, TJ, TZ and XN tissues. The circle size shows the number of REs圖1 對比兩種軟件預測REs共有部分的文氏圖Fig.1 The intersection Venn plot of REs predicted by comparing between JACUSA and SPRINT

圖2 探測的REs在牦牛染色體上的分布情況Fig.2 The distribution of REs in the whole genome

通過比對4種組織中REs的結(jié)果，對比SPRINT和JACUSA兩個軟件的預測結(jié)果，其在不同組織中的交集作為后續(xù)REs鑒定和下游分析的樣本(圖1)。其中圖1A、B、C、D分別是大腦、臀部肌肉、臀部脂肪和小腦組織中兩種軟件預測REs的文氏圖。

2.3 REs類型分布

針對檢測出的候選REs，按其在基因中的功能、變異產(chǎn)生的位置、編輯類型進行分型統(tǒng)計，結(jié)果見圖3。在所有REs編輯類型中，主要以A-G和T-C類型為主，廣泛存在嘌呤之間的轉(zhuǎn)換和嘧啶之間的轉(zhuǎn)換，而僅在臀脂組織中檢測A-T的編輯類型(圖3A)。在統(tǒng)計編輯產(chǎn)生的位置和所在RNA類型時，由于位點數(shù)據(jù)差距過大，因此本研究采用lg(REs數(shù)量+1)的形式對結(jié)果進行轉(zhuǎn)換。從圖3B中可以看出，多數(shù)REs位于基因下游區(qū)域、基因上游區(qū)域、基因間和內(nèi)含子區(qū)域(lg(RES+1)>2)。盡管4個組織中REs探測數(shù)目各不相同，但其變異類型、分布以及所屬RNA種類的比例相近。通過注釋觀察到，即使很多REs所屬RNA類型未知，但依然有很多REs所屬的RNA是編碼蛋白的mRNA或假基因(圖3C)。這些結(jié)果表明，REs與各組織中蛋白質(zhì)編碼活動具有密切關(guān)系。

圖3 4種組織中檢測的REs編輯類型、變異所在區(qū)域以及所屬RNA類型Fig.3 The editing types, location and RNA types of REs in 4 tissues

2.4 REs風險等級評估

對注釋后的REs統(tǒng)計編輯位點變異風險，共計獲得64個存在風險的REs。其中在大腦組織中發(fā)現(xiàn)一個高風險REs，這個REs可以導致SON基因提前終止翻譯蛋白質(zhì)。在大腦、小腦、臀脂和臀肌中分別發(fā)現(xiàn)13、8、9和2個中度風險REs，還有7、6、14和4個低風險REs(表2)。這個高度風險的REs位于牦牛1號染色體上的SON基因，具體變異位點在37 202 992 bp上，在Ensemble基因編號為ENSBGRG00000000139，SON是重要的DNA和RNA結(jié)合蛋白。

表2 4個組織中REs風險檢測結(jié)果Table 2 The prediction risk of REs in 4 tissues

對4個組織中檢測的非Modifier類型REs所在的基因進行GO富集分析，大腦中的REs相關(guān)基因主要富集在DNA代謝負調(diào)節(jié)條目中，小腦中的REs相關(guān)基因主要富集在神經(jīng)營養(yǎng)因子TRK正調(diào)節(jié)、蛋白多聚泛素化負調(diào)節(jié)以及腦源性神經(jīng)刺激受體條目中，臀脂中REs相關(guān)基因主要富集在細胞胚胎發(fā)生組分條目中，臀肌中REs相關(guān)基因主要富集在乙醇代謝條目中(圖4)。

A、B、C、D. 大腦、臀部肌肉、臀部脂肪和小腦組織中REs所在基因的GO聚類結(jié)果。圓圈的大小代表富集基因的多少，顏色越靠近紅色代表其P值越小A, B, C, D. The GO enrichment result of the REs in the DN, TJ, TZ and XN. The circle size indicates the number of enriched genes, the more red color show the fewer P value圖4 4種組織中REs所在基因的GO分析結(jié)果Fig.4 The GO analysis of REs located in genes in 4 tissues

2.5 REs的驗證

通過對兩個具有中度風險的REs進行DNA和cDNA克隆測序驗證，結(jié)果表明，EVA1C在3個類烏齊個體的小腦組織均發(fā)現(xiàn)了A-G的REs(圖5)。

A、B. EVA1C和SON基因克隆測序驗證結(jié)果A, B. The EVA1C and SON genes validation result圖5 通過DNA、cDNA克隆測序驗證REsFig.5 The REs validation by cloning and sequencing the gene DNA and cDNA

3 討 論

RNA編輯在生物界廣泛存在，其主要作用是高效利用同一段基因在不同時間或者不同組織中差異編輯，生成特異性的表達產(chǎn)物和翻譯產(chǎn)物。這種“一因多效”的機制可以增加生物多樣性，并與細胞生長、發(fā)育和分化密切相關(guān)。有關(guān)雞、豬、羊等畜禽的REs探測有大量報道[11,14,23-25]，但是在牦牛多組織中探測REs并未見報道。牦牛由于獨特的生存環(huán)境，其基因的利用以及與環(huán)境互作的關(guān)系更加復雜，在RNA水平的編輯可能更加獨特[4,26-27]。

本研究在4個組織樣本中均發(fā)現(xiàn)了特異性REs位點，其分布于除了性染色體外的所有常染色體。本研究還發(fā)現(xiàn)，臀部肌肉中REs數(shù)量最少，但在染色體中的分布并沒有相對集中。大腦組織中REs分布最多的是8號染色體，小腦中REs分布最集中的也是8號染色體，而在臀部肌肉和脂肪組織中REs分布主要集中在3、4、5號染色體，這說明REs分布與組織功能存在可能的關(guān)聯(lián)。

在諸多哺乳動物以及禽類的REs探測研究中，A-I廣泛存在于動物的生長發(fā)育過程中[28-29]，A-I是由雙鏈RNA腺苷脫氨酶催化腺苷A，C6位置的氨基水解脫掉氨基后形成次黃苷I，而次黃苷I在反轉(zhuǎn)錄成CDS完成測序時被識別成鳥苷G，因此,這個過程也被稱為A-G(A轉(zhuǎn)G)，本研究發(fā)現(xiàn)的編輯風險等級最高的位點就屬于這種編輯類型。A-I類型的RNA編輯方式在哺乳動物腦組織中普遍存在，有研究在小鼠中探測RNA編輯位點發(fā)現(xiàn)，這些編輯位點并不隨年齡和個體健康狀況改變，這些編輯位點似乎和遺傳有一定關(guān)系[11]。2012年，Ekdahl等[30]報道，A-I類型的RNA編輯方式對于動物腦發(fā)育的影響很大。2016年，Behm和?hman[24]對RNA編輯方式在哺乳動物大腦中的影響進行了深入研究發(fā)現(xiàn)，通過擴大編輯后亞型的多樣性，為蛋白質(zhì)提供了靈活的編輯，這為神經(jīng)元功能對環(huán)境的感知和適應起著重要作用，這種類型的編輯多集中在突觸傳遞的關(guān)鍵介質(zhì)中，對神經(jīng)元信號的整體影響顯著，對大腦發(fā)育至關(guān)重要[11,30]。

相對而言，本研究發(fā)現(xiàn)的主要編輯類型C-T是一種C-U的編輯過程，主要指胞苷脫氨基后形成尿苷，在反轉(zhuǎn)錄成CDS完成測序時被識別成胸腺苷T[31-32]。盡管本研究發(fā)現(xiàn)，A-G和C-T是牦牛RNA編輯的主要兩種方式，但其所占比與其他研究結(jié)果并不一致，本研究中，C-T類型的編輯屬于優(yōu)勢編輯類型，這可能是由于牦牛組織分化過程中存在更多與環(huán)境的互作。2019年，Zhang等[25]在對綿羊腎和脾組織RNA編輯位點的研究時發(fā)現(xiàn)，其主要編輯類型為A-G和C-T，與本研究結(jié)果相似，可能是由于同為草食動物的原因，同樣，這些RNA編輯位點所在的位置也主要集中在內(nèi)含子、基因間區(qū)域和基因下游區(qū)域，這些結(jié)果也與本研究相似。

通過GO分析，對比4個組織中REs所在基因的功能發(fā)現(xiàn)，大腦組織中REs所在基因(4個)主要富集在DNA分子代謝進程中(P<0.01)，這可能是由于牦牛大腦為應對復雜的外部環(huán)境產(chǎn)生的調(diào)整DNA分子代謝的特殊途徑；臀部脂肪組織中REs所在基因主要富集在TRK營養(yǎng)性受體信號通路、多聚泛素化蛋白調(diào)節(jié)以及腦源性神經(jīng)營養(yǎng)細胞受體通路(P<0.01)；臀部肌肉組織中REs(>30個)所在基因主要富集在胚胎細胞組成分化通路上，這表明，肌肉組織中REs參與了肌肉細胞分化的重要過程；而小腦組織中REs所在基因主要富集在乙醇代謝途徑上(P<0.01)。這些GO分析結(jié)果從REs所在基因功能的富集上表明，不同分化組織中REs所參與功能往往與其組織特定化功能有關(guān)[32]。

本研究發(fā)現(xiàn)的高風險RNA編輯位點位于牦牛SON基因上，SON是一種核蛋白，參與多種細胞過程，包括轉(zhuǎn)錄、信使RNA (pre-messenger RNA, mRNA)剪接和細胞周期調(diào)控[33-34]。近年來，研究者對SON在核組織和pre-mRNA剪接中的作用以及這些活動在維持細胞健康方面的影響進行了探究[35-36]。此外，SON在干細胞以及癌癥、流感和肝炎等多種疾病的發(fā)病過程中發(fā)揮了關(guān)鍵作用。本研究發(fā)現(xiàn)，在腦和肌肉脂肪組織之間，SON基因存在高風險編輯事件，這說明牦牛不同組織中利用SON基因來調(diào)控腦和肌肉脂肪組織中細胞健康以及諸多環(huán)境互作方面的基因表達。同樣，本研究發(fā)現(xiàn)的其他64個存在風險的REs，盡管通過風險評估不是最高，但在基因組表達產(chǎn)物調(diào)控中具有非常重要的生物學意義，是基因組應對復雜環(huán)境的“一因多用”。這種高效的“一因多用”機制如何在生物體內(nèi)具體調(diào)節(jié)，還有待進一步深入研究和探索。

4 結(jié) 論

本研究以類烏齊牦牛大腦、小腦、臀部脂肪及肌肉組織為試驗材料，通過高通量測序以及RNA-seq分析，對4個組織中相同基因的不同RNA編輯位點進行預測并對不同編輯類型進行分析。類烏齊牦牛大腦、小腦組織中的REs編輯過程和類型較為相近，而脂肪和肌肉組織中的REs編輯過程和類型比較相近。鑒定的一個高風險編輯位點使其所在基因的翻譯提前終止。這些發(fā)現(xiàn)將有助于研究牦牛各組織分化及發(fā)育過程中基因調(diào)控的分子規(guī)律。