吳盈萍，古力娜·巴克，李海英，趙曉鈺，馬 晨，姚瑩瑩，曹 妍，梅志勇，潘 璐

（新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，新疆烏魯木齊 830000）

新疆伊犁鵝（Yili Goose）是新疆優(yōu)良、珍稀的地方禽種，是我國(guó)唯一由灰雁（Anser anser）馴化而來的中小型家鵝品種[1]。新疆伊犁鵝的馴化歷史較短、馴化程度較低，其具有一定的短距離飛行能力，耐粗飼、抗病力強(qiáng)、抗寒耐熱且肉質(zhì)優(yōu)良，于2000 年8 月被農(nóng)業(yè)部（現(xiàn)為“農(nóng)業(yè)農(nóng)村部”，下同）列入第一批國(guó)家級(jí)畜禽品種資源保護(hù)名錄[2]。但新疆伊犁鵝的產(chǎn)蛋量較低，成年母鵝平均年產(chǎn)蛋量?jī)H為8～12 枚[3]，同時(shí)其就巢性很強(qiáng)，嚴(yán)重制約了地方鵝品種產(chǎn)蛋性能的提高。

轉(zhuǎn)錄組在研究畜禽生長(zhǎng)發(fā)育及生理活動(dòng)中發(fā)揮重要作用，也是表型關(guān)聯(lián)研究的重要方法。隨著轉(zhuǎn)錄組測(cè)序技術(shù)的快速發(fā)展，RNA-seq 技術(shù)已廣泛應(yīng)用于檢測(cè)新的轉(zhuǎn)錄本[4]、差異基因篩選、可變剪接[5-7]、基因結(jié)構(gòu)優(yōu)化[8-9]以及SNP 篩查[10-11]等相關(guān)研究。張蕾[12]等人利用轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)高郵鴨雙黃蛋高、低產(chǎn)組卵巢組織進(jìn)行測(cè)序，對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)控和注釋，篩選與高郵鴨雙黃蛋產(chǎn)蛋性能相關(guān)的SNP 位點(diǎn)，在所覆蓋的基因組數(shù)據(jù)中，檢測(cè)到位于基因外顯子區(qū)域的108 260～116 478 個(gè)SNP 位點(diǎn)，轉(zhuǎn)換類型總數(shù)極顯著高于顛換類型總數(shù)（P＜0.01）。韓坤鵬[13]等通過轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)京海黃雞卵巢進(jìn)行測(cè)序分析，基于所選參考基因組序列，共發(fā)掘4 431 個(gè)新基因，其中1 809 個(gè)新基因得到注釋。本研究應(yīng)用RNA-seq 技術(shù)對(duì)高、低產(chǎn)蛋量伊犁鵝性腺軸組織進(jìn)行SNP 查找、可變剪切篩查、基因結(jié)構(gòu)優(yōu)化以及新轉(zhuǎn)錄本預(yù)測(cè)，為進(jìn)一步解析產(chǎn)蛋量調(diào)控的分子機(jī)制提供理論依據(jù)。

1 材料與方法

1.1 樣品采集 本研究所用實(shí)驗(yàn)動(dòng)物由新疆額敏縣恒鑫實(shí)業(yè)有限公司提供。根據(jù)系譜資料的連續(xù)產(chǎn)蛋量記錄以及本次記錄的全期產(chǎn)蛋量，在體重相近、飼養(yǎng)條件相同的3 歲齡新疆伊犁鵝母鵝中篩選產(chǎn)蛋量具有極顯著差異的高、低產(chǎn)母鵝各4 只。高產(chǎn)母鵝的產(chǎn)蛋量分別為：14、16、18、19 個(gè)，低產(chǎn)母鵝的產(chǎn)蛋量分別為：5、6、7、5 個(gè)。屠宰后快速采集新疆伊犁鵝下丘腦、垂體、卵巢組織，用PBS 沖洗1～2 次之后，立即分裝放入含有RNA 保存液的凍存管中，做好標(biāo)記，樣品編號(hào)分別為X01-04（低產(chǎn)母鵝下丘腦）、X21-24（高產(chǎn)母鵝下丘腦）、C01-04（低產(chǎn)母鵝垂體）、C21-24（高產(chǎn)母鵝垂體）、L01-04（低產(chǎn)母鵝卵巢），L21-24（高產(chǎn)母鵝卵巢），4℃靜置過夜，隨后再放到-80℃冰箱保存，用于組織樣總RNA 提取。

1.2 組織樣總RNA 提取 采用Trizol 法分別提取24 個(gè)樣品的RNA，并分別用瓊脂糖凝膠電泳、Nanodrop、Qubit、Agilent 2100 4 種方法對(duì)其質(zhì)量進(jìn)行檢測(cè)以保證使用合格的樣品進(jìn)行轉(zhuǎn)錄組測(cè)序。

1.3 cDNA 文庫(kù)的構(gòu)建及Illumina 測(cè)序 樣品檢測(cè)合格后進(jìn)行cDNA 文庫(kù)構(gòu)建，構(gòu)建完成后先使用Qubit 2.0進(jìn)行初步定量，稀釋文庫(kù)至1 ng/μL，隨后使用Agilent 2100 對(duì)文庫(kù)的插入片段長(zhǎng)度進(jìn)行檢測(cè)，插入片段長(zhǎng)度符合預(yù)期后，使用Q-PCR 方法對(duì)文庫(kù)的有效濃度進(jìn)行準(zhǔn)確定量（文庫(kù)有效濃度＞2 nM），以保證文庫(kù)質(zhì)量。庫(kù)檢合格后，進(jìn)行Illumina HiSeq 2500 高通量測(cè)序平臺(tái)測(cè)序，測(cè)序讀長(zhǎng)為雙末端測(cè)序150 bp。

1.4 測(cè)序數(shù)據(jù)分析 高通量測(cè)序得到的原始圖像數(shù)據(jù)文件經(jīng)CASAVA 堿基識(shí)別（Base Calling）分析得到原始測(cè)序序列（Raw reads），對(duì)Raw reads 進(jìn)行數(shù)據(jù)過濾，得到高質(zhì)量的過濾后測(cè)序數(shù)據(jù)（Clean reads），后續(xù)分析都基于Clean reads。使用序列高效比對(duì)軟件HISAT2.0.4 對(duì)各樣品測(cè)序數(shù)據(jù)與參考基因組（https://www.ncbi.nlm.nih.gov/assembly/GCF_000002315.6）進(jìn)行序列比對(duì)，從而獲取Clean reads 在參考基因組上的位置及測(cè)序樣品特有的序列信息。

1.5 基因結(jié)構(gòu)分析 SNP 分析：首先通過samtools 和picard-tools 等工具對(duì)比對(duì)結(jié)果進(jìn)行染色體坐標(biāo)排序、reads 去重復(fù)等處理，然后通過變異檢測(cè)軟件samtools分別查找NGS 數(shù)據(jù)與參考序列區(qū)別（SNP calling），最后進(jìn)行過濾，得到SNP 信息。

可變剪切分析：以每個(gè)進(jìn)行差異可變剪切分析的比較組為單位，本研究首先統(tǒng)計(jì)發(fā)生的可變剪切事件的種類及數(shù)量，然后分別計(jì)算每類可變剪切事件表達(dá)量，最后對(duì)每類可變剪切事件進(jìn)行差異分析。

新轉(zhuǎn)錄本預(yù)測(cè)及基因結(jié)構(gòu)優(yōu)化：將所有測(cè)序reads數(shù)據(jù)的基因組定位結(jié)果放到一起，用Cufflinks 進(jìn)行組裝，然后用Cuffcompare 和已知的基因模型進(jìn)行比較，可以發(fā)現(xiàn)新基因（相對(duì)于原有基因注釋文件）和已知基因新的外顯子區(qū)域，以及對(duì)已知基因的起始和終止位置進(jìn)行優(yōu)化。

2 結(jié)果

2.1 RNA 質(zhì)量檢測(cè)及測(cè)序數(shù)據(jù)質(zhì)量檢測(cè)

2.1.1 RNA 質(zhì)量檢測(cè) 由圖1 可知，RNA 分子28 s 和18 s 條帶清晰，且亮度比例約為2:1，表明RNA 質(zhì)量好、完整性好、RNA 降解率低。

圖1 總RNA 瓊脂糖凝膠電泳檢測(cè)結(jié)果

2.1.2 測(cè)序數(shù)據(jù)及其質(zhì)量控制 由表1 可知，通過測(cè)序，高、低產(chǎn)新疆伊犁鵝分別獲得了579 563 136 個(gè)reads和596 933 010 個(gè)reads；總堿基數(shù)為176.47 G，每個(gè)樣品均獲得了5.98 G 以上的堿基，Q20 在95.10%以上，Q30 在88.57%以上，GC 含量在48.28%～50.57%之間，說明轉(zhuǎn)錄組測(cè)序得到的結(jié)果質(zhì)量較高，滿足后續(xù)結(jié)果分析的需求。

表1 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估

2.2 基因結(jié)構(gòu)分析

2.2.1 SNP 分析 由表2 可知，在24 個(gè)樣品上分別得到208 788～539 612 個(gè)可能的SNP 位點(diǎn)。在各樣品中，大多數(shù)堿基取代均為轉(zhuǎn)換大于顛換，其中轉(zhuǎn)換型SNP 所占百分比為72.22%～73.09%；顛換型SNP 所占百分比為26.9%～27.78%。同時(shí)，24 個(gè)樣品的雜合SNP 位點(diǎn)所占百分比為17.83%～35.75%，純合比例為64.25%～82.17%。

表2 SNP 位點(diǎn)統(tǒng)計(jì)

2.2.2 可變剪切事件分類和數(shù)量統(tǒng)計(jì) 性腺軸組織樣品的各種可變剪切方式的統(tǒng)計(jì)結(jié)果如圖2 所示。在定量過程中，rMATs 采取了2 種定量方式，分別為：Junction Count only、Reads on target and junction counts。本文只對(duì)其中5 類可變剪切事件包括外顯子跳躍（SE:Skipped exon）、外顯子選擇性跳躍（MXE:Mutually exclusive exon）、第一個(gè)外顯子可變剪切（A5SS:Alternative 5'splice site）、最后一個(gè)外顯子可變剪切（A3SS:Alternative 3'splice site）和內(nèi)含子滯留（RI:Retained intron）進(jìn)行了鑒定。

圖2 可變剪切分類和數(shù)量統(tǒng)計(jì)

在樣品X2 vs X0（高產(chǎn)組下丘腦vs 低產(chǎn)組下丘腦）、C2 vs C0（高產(chǎn)組垂體vs 低產(chǎn)組垂體）和L2 vs L0（高產(chǎn)組卵巢vs 低產(chǎn)組卵巢）中同時(shí)使用Junction Counts 和reads on target 檢測(cè)到的可變剪切事件分別為20 531、21 634、22 044 個(gè)，表明許多基因都包含著多種可變剪切方式，從而造成基因在表達(dá)水平上的變化，由于基因剪切方式的不同，使得生物能夠表達(dá)多種蛋白，進(jìn)而引起表型的多樣性。在5 類可變剪切事件中，X2 vs X0、C2 vs C0 和L2 vs L0 均為SE 事件比例最大，MXE 次之，RI 事件比例最小。

2.2.3 已注釋基因結(jié)構(gòu)優(yōu)化 由表3 可知，本研究共有11 383 個(gè)基因的起始位置和終止位置得到了優(yōu)化，其中正鏈5 484 條，負(fù)鏈5 899 條。

表3 部分已知基因結(jié)構(gòu)優(yōu)化

2.3 新轉(zhuǎn)錄本的發(fā)掘及注釋

2.3.1 新轉(zhuǎn)錄本的發(fā)掘 本研究共發(fā)掘新基因3 119 個(gè)，獲得注釋的基因902 個(gè)。本實(shí)驗(yàn)垂體和卵巢組織差異表達(dá)基因篩選過程中將Q value＜0.05 作為篩選標(biāo)準(zhǔn)，篩選差異基因。當(dāng)Q value＜0.05 基因作為下丘腦差異表達(dá)閾值時(shí)，篩選到的差異新基因?yàn)?。因此，將下丘腦差異表達(dá)基因的篩選閾值調(diào)整為P value值＜0.005。差異表達(dá)新基因共64 個(gè)，其中49 個(gè)新基因表達(dá)上調(diào)，15個(gè)新基因表達(dá)下調(diào)（表4）。

表4 差異表達(dá)新基因數(shù)目統(tǒng)計(jì)表

新轉(zhuǎn)錄本在染色體上的分布如圖3 所示，新轉(zhuǎn)錄本分布前三的染色體分別為NW_013185656.1 染色體343條、NW_013185751.1 染色體321 條、NW_013185661.1染色體314 條。

圖3 轉(zhuǎn)錄本在染色體上的分布

2.3.2 新轉(zhuǎn)錄本的注釋 利用Blast2Go 軟件對(duì)篩選到的新轉(zhuǎn)錄本進(jìn)行GO 富集分析，3 119 個(gè)新轉(zhuǎn)錄本注釋到細(xì)胞組成、生物學(xué)過程和分子功能。新轉(zhuǎn)錄本共富集到2 500 條GO Term，細(xì)胞組成占60.96%，生物學(xué)過程占26.16%，分子功能占12.88.%。其中，顯著富集的有多生物過程、DNA 拓?fù)洚悩?gòu)酶活性和DNA 拓?fù)渥兓▓D4）。在細(xì)胞組成方面主要富集到其他生物細(xì)胞、細(xì)胞外區(qū)域、線粒體膜部分、線粒體內(nèi)膜等；在生物學(xué)過程方面主要富集到調(diào)節(jié)免疫系統(tǒng)過程、ATP 生物合成過程、交配行為、繁殖行為、繁殖過程等；在分子功能方面主要富集到通道調(diào)節(jié)器活性、激素活性、結(jié)構(gòu)分子活性、雌激素受體活性等。

利用KOBAS（2.0）軟件對(duì)KEGG 注釋通路進(jìn)行分析，3 119 個(gè)新轉(zhuǎn)錄本注釋到49 條代謝通路。代謝途徑見圖5，顯著富集的有泛素介導(dǎo)的蛋白質(zhì)水解和核苷酸切除修復(fù)途徑。其中，核苷酸切除修復(fù)途徑富集因子最大。除顯著富集途徑外，富集排名前十的途徑有：RNA 轉(zhuǎn)運(yùn)、精氨酸和脯氨酸代謝、組氨酸代謝、氨基酸的生物合成、氨?；?tRNA 的生物合成、β-丙氨酸代謝、鞘脂代謝、糖酵解/糖異生、抗壞血酸鹽和藻酸鹽代謝、核糖體。

圖5 性腺軸組織差異基因KEGG 富集散點(diǎn)圖

3 討 論

3.1 基因結(jié)構(gòu)分析 SNP 是指在基因組上由單個(gè)核苷酸變異形成的遺傳標(biāo)記，其數(shù)量很多，多態(tài)性豐富。從理論上看，每個(gè)SNP 位點(diǎn)都可以有4 種不同的變異形式，但實(shí)際上發(fā)生的只有2 種，即轉(zhuǎn)換和顛換。鄭杰[14]等研究顯示，犏牛新鮮囊胚和凍融囊胚兩樣本中轉(zhuǎn)換型SNP 比例都明顯高于顛換型，轉(zhuǎn)換/ 顛換比值分別為2.57 和2.45。黃育雯[15]對(duì)安格斯牛下丘腦和垂體的每個(gè)樣品的SNP 轉(zhuǎn)換和顛換進(jìn)行計(jì)數(shù)，結(jié)果顯示，各樣品堿基轉(zhuǎn)換的可能性均超過70%，遠(yuǎn)遠(yuǎn)超過了顛換，同時(shí)，雜合SNP 位點(diǎn)超過37.68%。本研究在性腺軸組織24 個(gè)樣品中共獲得1 048 575 個(gè)SNP 位點(diǎn)，其中轉(zhuǎn)換型約占72.22%～73.09%，顛換型約占26.9%～27.78%，轉(zhuǎn)換與顛換的比值在2.60～2.72 之間。在高產(chǎn)伊犁鵝的垂體和卵巢組織中SNP 數(shù)目明顯高于低產(chǎn)伊犁鵝，提示SNP 在伊犁鵝產(chǎn)蛋量方面起重要的調(diào)控作用。

圖4 性腺軸組織差異基因GO 富集柱狀圖

可變剪切是調(diào)節(jié)基因表達(dá)和蛋白質(zhì)組多樣性的主要因素之一，在轉(zhuǎn)錄組中大量存在[16]。真核生物中存在著大量的可變剪切，動(dòng)物細(xì)胞中也存在著組織特異性[17]。李慧鋒[18]等對(duì)白來航母雞首次產(chǎn)蛋前后的肝臟組織進(jìn)行分析，得到可變剪切事件在產(chǎn)蛋前后分別為14 518個(gè)和14 683 個(gè)，其中外顯子跳躍事件最多，5' 可變剪切位點(diǎn)事件最少。Li 等[19]報(bào)道，20 周齡青年雞和30周齡產(chǎn)蛋雞肝臟的可變剪切類型主要為內(nèi)含子滯留（RI 32%）和外顯子跳躍（SE 24%）。本研究發(fā)現(xiàn)，在5類可變剪切事件中，X2 vs X0、C2 vs C0 和L2 vs L0均為SE 事件比例最大，MXE 次之，RI 事件比例最小，表明高、低產(chǎn)新疆伊犁鵝性腺軸組織的可變剪切類型中，SE 事件為主要的剪切方式，是導(dǎo)致蛋白質(zhì)形成差異的主要機(jī)制，也可能是表現(xiàn)產(chǎn)蛋量差異的主要機(jī)制。

此外，通過RNA-seq 越來越多的新轉(zhuǎn)錄本被發(fā)現(xiàn)，顧麗紅[20]等通過轉(zhuǎn)錄組測(cè)序技術(shù)，對(duì)180 日齡文昌雞和隱性白羽克洛雞胸肌進(jìn)行測(cè)序，在4 個(gè)樣品中發(fā)現(xiàn)了289～339 個(gè)新轉(zhuǎn)錄本。葉保國(guó)[21]等對(duì)北京鴨6 個(gè)組織的RNA-seq 數(shù)據(jù)進(jìn)行新轉(zhuǎn)錄本分析，預(yù)測(cè)得到了8 141個(gè)新轉(zhuǎn)錄本，單個(gè)樣品可以預(yù)測(cè)得到3 636～4 918 個(gè)新轉(zhuǎn)錄本。裴星朝[22]對(duì)番鴨下丘腦、垂體和卵巢3 種組織進(jìn)行轉(zhuǎn)錄組深度測(cè)序，對(duì)基因結(jié)構(gòu)優(yōu)化分析，發(fā)掘了1 427 個(gè)新基因。本研究顯示，轉(zhuǎn)錄組測(cè)序結(jié)果中新轉(zhuǎn)錄本的表達(dá)譜分析揭示了這些新轉(zhuǎn)錄本也參與了伊犁鵝的繁殖過程。本研究發(fā)現(xiàn)在性腺軸組織中總共有3 119個(gè)新轉(zhuǎn)錄本，獲得注釋的基因902 個(gè)。差異表達(dá)新基因64 個(gè)，其中49 個(gè)差異新基因表達(dá)上調(diào)，15 個(gè)差異新基因表達(dá)下調(diào)。

3.2 功能分類分析 為揭示新轉(zhuǎn)錄本的生理功能，本實(shí)驗(yàn)對(duì)篩選到的新轉(zhuǎn)錄本進(jìn)行了GO 富集分析，其中細(xì)胞組成所占比例最大。此結(jié)果表明，參與伊犁鵝繁殖生理活動(dòng)的分子組分廣泛分布于細(xì)胞組分，這提示細(xì)胞組分在繁殖生理活動(dòng)中扮演著重要角色。此外，參與伊犁鵝性腺軸生理活動(dòng)的新轉(zhuǎn)錄本顯著富集在多生物過程、DNA 拓?fù)洚悩?gòu)酶活性和DNA 拓?fù)渥兓?，說明這些生物學(xué)過程對(duì)于伊犁鵝繁殖性能具有重要調(diào)控作用。

KEGG 信號(hào)通路分析表明，這些新轉(zhuǎn)錄本共涉及到49 條代謝通路，其中泛素介導(dǎo)的蛋白水解和核苷酸切除修復(fù)途徑通路顯著富集。泛素-蛋白酶體途徑是一種依賴ATP 進(jìn)行的、具有高度特異性和選擇性的蛋白質(zhì)降解途徑，主要負(fù)責(zé)降解細(xì)胞內(nèi)超過80% 的正?；虍惓５鞍踪|(zhì)，泛素介導(dǎo)的蛋白質(zhì)降解系統(tǒng)涉及細(xì)胞的多種重要生理功能，參與了基因轉(zhuǎn)錄、蛋白質(zhì)翻譯、信號(hào)傳導(dǎo)、細(xì)胞周期控制以及生長(zhǎng)發(fā)育等幾乎所有的生命活動(dòng)過程[23-26]。核苷酸切除修復(fù)為DNA 修復(fù)途徑之一，且核苷酸切除修復(fù)基因具有高度多態(tài)性。修復(fù)基因突變后可能引起DNA 修復(fù)能力的改變，從而導(dǎo)致機(jī)體免疫系統(tǒng)的改變[27]。生理上，DNA 修復(fù)過程需要多種酶作用，與修復(fù)蛋白表達(dá)水平密切相關(guān)，而后者受轉(zhuǎn)錄水平調(diào)控。因此，DNA 修復(fù)基因轉(zhuǎn)錄水平可以反映機(jī)體受損時(shí)細(xì)胞DNA 修復(fù)受損的能力[28]。新轉(zhuǎn)錄本富集到的信號(hào)通路分析只能為性腺軸生物學(xué)功能的信號(hào)調(diào)控途徑提供一個(gè)方向，但是新疆伊犁鵝繁殖性能的具體調(diào)控過程還需要對(duì)涉及到相關(guān)通路的對(duì)應(yīng)基因進(jìn)行更深入的研究。

4 結(jié) 論

在高、低產(chǎn)蛋量新疆伊犁鵝性腺軸中共篩選出1 048 575個(gè)SNP位點(diǎn)，64 209個(gè)可變剪切，優(yōu)化了11 384個(gè)基因，篩選到3 119 個(gè)新轉(zhuǎn)錄本，獲得注釋的轉(zhuǎn)錄本902 個(gè)。這些新轉(zhuǎn)錄本的發(fā)現(xiàn)和功能得以注釋為進(jìn)一步探索挖掘新基因奠定基礎(chǔ)，也為伊犁鵝產(chǎn)蛋量相關(guān)基因的篩選提供科學(xué)依據(jù)。