基于RNA-seq 對復合益生菌作用于肉雞回腸中新轉錄本預測及可變剪接、SNP 分析

曹威榮， 王燕飛， 張若男， 賈 浩， 劉 璇， 張利環

（1.朔州職業技術學院，山西 朔州 036002；2.山西農業大學生命科學學院，山西 晉中 030801）

黃羽肉雞抗病力和對劣質飼料的適應性強，易飼養，具有較高的繁殖性能，肉質風味與口感也較佳。我國黃羽肉雞種類繁多，其中黃麻肉雞體型圓潤、整體發育勻稱、雞皮較薄、皮下脂肪適中，肉中富含較高蛋白和維生素A，具早熟易肥、肉質緊實滑嫩等特點，因此在全國被廣泛養殖（阿卜杜熱合曼·達尼亞爾，2018）。 食用動物攝取的抗生素，以及產生的食物中的抗生素殘留， 已被公認為人類耐藥性迅速增加的主要原因之一 （Davies 等，2010）。 濫用抗生素飼喂食用動物，會直接選擇對抗生素有抗藥性的微生物， 并把食用動物的系統轉變成抗生素抗藥性基因的貯存庫（Claesson 等，2012）。因此，減少食用抗生素的攝入量，消除食用動物農業中的抗生素殘留， 已成為食品安全和公共衛生的首要任務。

目前，在養雞業中，飼料中抗生素的替代品包括纖維降解酶、益生菌、益生元、共生元和植物生物元素等（Sopkova 等，2017）。 其中益生菌生產成本低，在不同種類的宿主動物中應用范圍廣，具有優勢（Gaggia 等，2010）。益生菌被定義為在機體內定植一段時間后可改善腸道微生態并對宿主動物產生積極影響的活性微生物 （Sanders 等，2008）。干酪乳桿菌、 嗜酸乳桿菌和雙歧桿菌一起被稱為“健康三益菌”。 乳酸菌是一種非致病性和毒素性的革蘭氏陽性發酵菌，無芽孢與鞭毛，菌落白色粗糙，區別于其他厭氧菌，大多數可以在氧氣存在下生長。雙歧桿菌是定居于人胃腸道的微生物之一，雖不屬于乳酸菌， 但可以產生乳酸。 在飲用含有1×109CFU 的干酪乳桿菌溶液時， 益生菌因其特性對胃液、 水解酶和膽汁酸具有高耐受性并可在胃腸道中存活并停留3 d（Radicioni 等，2019）。 乳酸菌在腸道中因可分泌諸如細菌素的抗菌物質而具有抗炎性和抑菌特性，產生的酸可降低腸道pH并防止產氣莢膜梭菌的增殖，有利于腸道健康。

RNA-seq 數據的高分辨率不僅可以探索基因表達， 還可以預測動物中的可變剪接事件（Li等，2019）。 RNA-seq 可以反映出mRNA、smallRNA、noncodingRNA 等或者其中一些的表達水平，有助于查看突變/SNP 和基因表達隨時間的變化（Maher 等，2009），還可以查看miRNA、tRNA 和核糖體分析等（Ingolia 等，2012）。SNP 是一種基于核苷酸變異性的有效遺傳標記， 在動物遺傳學研究和育種中有著廣泛的應用（Li 等，2019）。 可變剪接是真核生物中的一種正?，F象， 其帶來了重要的蛋白質多樣性并允許基因產生不同的mRNA轉錄本， 這些轉錄本被翻譯成不同的蛋白質（Li等，2020）。

本試驗對肉雞回腸進行了新轉錄本的發掘和注釋、可變剪接和SNP 分析，以期為益生菌飼喂肉雞新基因的挖掘提供數據基礎， 進而更好的探究其分子機制。

1 材料與方法

1.1 試驗動物及樣品采集 干酪乳桿菌、嗜酸乳桿菌與雙歧桿菌凍干菌粉有效活菌數為1×1010CFU/g， 凍干菌粉購自陜西瑞茂生物科技有限責任公司，復合菌制劑（干酪乳桿菌：嗜酸乳桿菌：雙歧桿菌=1:1:2）（Zhang 等，2022）。200 只1 日齡黃麻肉雞（公母各半，購自太谷縣田建勝養殖場）隨機分為2 組，每組5 個重復，每個重復20 只。 對照組正常飲水，益生菌組飲用水中補充1%復合菌制劑，試驗分為生長前期（1 ～21 d）、后期（22 ～42 d）兩個階段。 基礎日糧配制參照NRC（1994）和《NY/T 33-2004 雞飼養標準》（表1）。 42 日齡時進行肉雞屠宰試驗，取回腸后用生理鹽水將其沖洗干凈，分裝后液氮速凍，并于-80 ℃保存。

表1 基礎日糧組成和營養水平

1.2 RNA-Seq 測序 從對照組和益生菌組隨機抽取3 個樣本 （編號為C1、C2、C3 和P1、P2、P3）送至上海美吉生物醫藥科技有限公司，使用Illumina Novaseq 6000 測序平臺進行轉錄組測序（RNA sequencing，RNA-Seq）。

1.3 新轉錄本的發掘 使用Cufflinks 軟件對Mapped Reads 進行拼接，并與原有的基因組注釋信息進行比較，尋找原來未被注釋的轉錄區，發掘該物種的新轉錄本和新基因。

1.4 新轉錄本功能注釋 使用BLAST 軟件將發掘的新轉錄本分別與KEGG（http://www.genome.jp/kegg/）、COGs （http://www.ncbi.nlm.nih.gov/COG/）、GO （http://www.geneontology.org/）、NR（ftp://ftp.ncbi.nlm.nih.gov/blast/db/）、Swiss -Prot （ftp://ftp.uniprot.org/pub/databases/uniprot/current_release）、Pfam（http://pfam.xfam.org/）等數據庫進行比對來獲得測序所得肉雞新轉錄本的注釋信息。

1.5 可變剪接分析及差異可變剪接基因分析使用rMATS（http://rnaseq-mats.sourceforge.net/index.html）對RNA-seq 數據進行可變剪接分析。利用GO、KEGG 數據庫進行差異可變剪接基因富集分析。

1.6 SNP 分析 使用Sentieon 進行SNP 分析，得到SNP 信息。

2 結果與分析

2.1 測序數據質控結果及與參考基因組比對分析 如表2 所示，樣本過濾后的Clean reads 數占Raw reads 的比例為98.26% ～98.91%， 且Clean bases 占Raw bases 的比例為96.16% ～97.48%；測序堿基平均錯誤率均為0.026%以下；堿基識別準確率在99%以上的堿基所占百分比均在97.51%以上， 堿基識別準確率在99.9%以上的堿基所占百分比均在93.27%以上，GC 含量在樣本中所占比例為49.09% ～51.12%，沒有分離現象，說明測序產出數據質量較高，符合測序要求。

表2 序列質量和比對信息匯總統計

測序數據與參考基因組的比對率為89.04% ～91.48%， 說明不存在污染且比對基因組正確；在參考序列上有唯一位置的比對率為86.88% ～89.35%， 在參考序列上有多個位置的比對率為2.13% ～2.39%，說明基因比對效果較好，且基因覆蓋均一度較好。

2.2 新轉錄本的發掘及注釋 如圖1 所示，將測序數據與參考基因組注釋信息進行比對， 共發現1047 個新基因， 注釋到GO、KEGG、COG、NR、Swiss-Prot、Pfam 數據庫中的基因分別有276、196、323、567、250、247 個。

圖1 各數據庫注釋的新基因數

如表3 所示， 將高通量測序數據與參考基因組注釋信息進行比較，共發現20176 個新轉錄本，注釋到GO、KEGG、COG、NR、Swiss-Prot、Pfam 數據庫的新轉錄本數分別為13739、13488、17693、18764、17499、16655。

表3 各數據庫注釋的新轉錄本數

如圖2 所示， 轉錄組測序得到的新基因共276 個基因在GO 數據庫中得到歸類注釋。 分子功能主要集中在結合性、催化活性等多個分類中；在細胞組分分類中，主要集中于生物膜組分、細胞成分、細胞器等中；在生物學過程分類中，主要集中在細胞過程、代謝過程、生物調控等中。

圖2 新基因GO 分類注釋圖

2.3 肉雞可變剪接分析及差異可變剪接基因分析 如圖3 所示， 其中SE 占可變剪接類型的比例最高， 達到70.11% ～70.28%，A5SS、A3SS、MXE、RI 分別占各個階段可變剪接事件總量的6.80% ～7.75%、1.04% ～1.19%、5.80% ～6.12%、4.49% ～5.18%。

圖3 可變剪接事件統計圖

如圖4 所示，共篩選到1131 個顯著的差異剪接基因， 在SE、A5SS、A3SS、MXE、RI 事件中鑒定到的差異剪接基因數分別為588、139、176、136、92。

圖4 差異剪切基因統計圖

如圖5 所示，在SE 中外顯子排斥有327 個，外顯子包含有410 個； 在A3SS 中外顯子排斥有99 個，外顯子包含有101 個；在A5SS 中外顯子排斥有72 個，外顯子包含有87 個；在MXE 中外顯子排斥有138 個，外顯子包含有117 個；在RI 中外顯子排斥有55 個，外顯子包含有58 個。

圖5 可變剪接事件數目外顯子納入和排斥統計圖

如圖6 所示， 對獲得的差異剪接基因進行GO 富集分析，圖中展示了前20 個GO 條目，富集到了與代謝和免疫相關的GO 條目， 如大分子代謝過程的調控、 磷脂酰肌醇磷酸4-磷酸酶活性、細胞代謝過程、免疫反應的正調控、大分子代謝過程的正調控。

圖6 差異剪切基因GO 富集柱狀圖

如圖7 所示，差異可變剪接基因KEGG 富集到了以下通路， 胞吞作用、MAPK 信號通路、Rap1信號通路、GnRH 信號通路、 黏附連接、TNF 信號通路、磷脂酰肌醇信號系統、ErbB 信號通路等。

圖7 差異剪切基因KEGG 富集氣泡圖

2.4 SNP 分析 由表4 可知， 轉換SNP 所占百分比為73.53% ～73.94%； 顛換型SNP 所占百分比為26.05% ～26.47%， 在各樣品中大多數堿基取代均為轉換大于顛換。

表4 SNP 類型統計

3 討論

轉錄組學研究的是在特定發育階段或生理條件下細胞RNA 水平上基因的表達和變化情況，進而查找探究基因的功能， 它在新轉錄本預測方面有廣泛應用（Chen 等，2011），還用于小菜蛾、帝王蝶等昆蟲新測序基因組的基因注釋 （You 等，2013；Zhan 等，2011）。 Izadnia 等 （2019） 利用RNA-seq 技術鑒定羅斯708 和伊斯法罕雞兩種雞肝臟與剩余采食量（RFI）相關的關鍵基因和新的轉錄體。Liu 等（2020）在鷓鴣雞中37.6%的蛋白質編碼基因中發現了新的轉錄體。 Tang 等（2020）對來自Valgus-varus 畸形哈伯德肉雞和哈伯德肉雞的6 個雞脾臟樣本進行了轉錄組測序發現了7943 個新的轉錄本。 本研究共發現20176 個新轉錄本，1047 個新基因， 新基因中共276 個基因在GO 數據庫中得到歸類注釋。

目前， 轉錄組測序已成為分子生物學研究的常用工具， 廣泛應用于檢測新的轉錄本（Huang等，2010）、差異基因篩選（Filichkin 等，2010）、可變剪接 （Lan 等，2014） 以及SNP 篩查（Martinez等，2017）等相關研究。 真核生物中存在著大量的可變剪接（Li 等，2016），動物細胞中也存在著組織特異性（Montelli 等，2015）。Hu 等（2021）對北京鴨胸肌和腿肌進行轉錄組測序分析發現胸肌肌肉組織的可變剪接數量為32256 ～45707，腿部肌肉組織的可變剪接數量在不同時間點為35841 ～47889。Li 等（2015）發現，A3SS、A5SS、RI、SE 是20周齡青年雞和30 周齡產蛋雞肝臟中主要剪切的模式， 四類剪切事件占總數的96.5%，MEX 僅占3.5%。 Hu 等（2021）預測了漢中麻鴨不同生長階段的骨骼肌可變剪接事件，在17、21、27 d 和6 個月大時的骨骼肌中可變剪接最多 （＞10000）。Daines 等（2010）利用RNA-Seq 數據與黑腹果蠅基因組比對，識別了至少4618 個基因（31%）發生了可變剪接， 確定了7776 個潛在的新剪切事件，最常見的帶注釋的事件是外顯子跳躍（SE）。 本研究中SE 的可變剪接類型的比例最高， 共篩選到1131 個顯著的差異剪接基因，并對差異可變剪接基因進行GO 富集分析， 主要富集到代謝和免疫GO term。

SNP 是指在基因組上由單個核苷酸的變異所引起的DNA 序列多態性，其數量很多，多態性豐富， 主要發生的有2 種， 即轉換和顛換。 黃育雯（2020） 在安格斯牛下丘腦和垂體研究中發現，堿基轉換的可能性超過了顛換。 Zhang 等（2019）檢測了大橫肉雞MSTN 基因的外顯子序列， 共檢測到5 個SNPs。 張蕾等（2019）利用轉錄組測序技術對6 個高郵鴨卵巢組織進行研究發現， 轉換類型總數極顯著高于顛換類型總數。 姚娜等（2012）發現， 西藏班戈絨山羊和遼寧絨山羊山羊兩品種間有5391 個SNP 位點。本研究在6 個樣品上分別得到了可能的SNP 位點，SNP 轉換類型總數大于顛換類型總數。

4 結論

本試驗對肉雞的回腸進行了RNA-seq 測序分析，與參考基因組和參考基因比對發現了1047個新基因，20176 個新轉錄本； 共篩選到1131 個顯著的差異剪接基因， 在SE、A5SS、A3SS、MXE、RI 事件中鑒定到的差異剪接基因數分別為588、139、176、136、92；SNP 中轉換類型總數大于顛換類型總數。 本研究為發掘復合益生菌作用于肉雞回腸中的可變剪接事件提供基礎， 為進一步了解復合益生菌作用于肉雞回腸中新轉錄本的發現和SNP 位點提供了數據支撐。