基于轉錄組測序挖掘廣西麻雞生長發育相關基因

陳 婷，盧晟盛，楊麗麗，李添寶，張 敏，江明生

(1.廣西大學動物科學技術學院，南寧 530004；2.廣西大學農牧產業發展研究院，南寧 530004)

廣西麻雞具備優質肉雞的良好產肉性能、風味佳、營養高，近年來在市場占有份額越來越大，母雞的營養優于公雞[1]，胸肌和腿肌具有豐富的氨基酸，且腿肌比胸肌風味更佳[2]。腿肌是雞主要的產肉部位之一，其在肌纖維生長發育[3]、屠宰性能和肉品質[4]、風味[5]等方面有相關報道。廣西大學動物科技學院動物遺傳育種實驗室在前期生產試驗中發現，廣西麻雞母雞生長速度慢、產蛋較少、母雞就巢性強、繁殖性能差等問題，其生長高峰期在60日齡左右。因此，加強對廣西麻雞的遺傳育種研究具有重要的意義。

雞的生長發育是一個重要的經濟性狀，在家禽生產中，隨著家禽選育程度的不斷提高，對雞生長發育的調控機制應該得到重視。生長發育受遺傳和環境等多種因素的影響，在一系列細胞因子和轉錄因子的調控下按照嚴格的進程完成，出生后肌肉的形成對生長起決定作用。通過轉錄組學研究方法已有學者對雞的肉質[6]、脂肪沉積[7]、肌內脂肪代謝[8]、精子活力[9]、長鏈脂肪酸變化[10]、繁殖性狀[11]、肝臟抗氧化酶活性[12]、多不飽和脂肪酸[13]等相關差異表達基因進行研究以加速品種的選育和改良。在雞生長發育的研究中，腿肌轉錄組測序的研究相對較少。本研究以廣西麻雞母雞出生時(1日齡)和生長高峰期(60日齡)腿肌為研究對象，對2個不同生長階段的腿肌組織進行轉錄組測序，從mRNA水平上對差異表達基因進行研究，對肌肉生長發育過程中差異表達基因進行篩選，并進行功能注釋以挖掘影響生長發育的關鍵基因，以期為闡明廣西麻雞生長發育的分子機制提供參考。

1 材料與方法

1.1 樣品采集

1和60日齡健康廣西麻雞母雞各3只，由廣西南寧良慶區農利來種禽有限公司提供。屠宰后，采集腿肌組織(1日齡腿肌組織分別記為1LM-1、1LM-2、1LM-3；60日齡腿肌組織分別記為60LM-1、60LM-2、60LM-3)，迅速放入―80 ℃保存備用。

1.2 RNA提取

用TRIzol法分別提取6只雞的腿肌組織總RNA，利用Bioanalyzer 2100和RNA 6000 Nano LabChip Kit分析總RNA質量和純度，選取D260 nm/D280 nm>1.8、D260 nm/D230 nm>1.0、RNA濃度>50 ng/μL的樣品用于后續試驗。

1.3 轉錄組測序分析

樣品檢驗合格后，以mRNA為模板合成雙鏈cDNA，雙鏈cDNA經過末端修復，加A尾修飾并連接測序接頭，再利用Oligo磁珠對其片段大小進行篩選和純化，用UDG酶消化雙鏈。PCR反應程序：95 ℃預變性3 min；98 ℃變性15 s，60 ℃退火15 s，72 ℃延伸30 s，共8個循環；72 ℃終延伸5 min，形成300 bp±50 bp的cDNA文庫。使用Illumina NovaseqTM6000進行轉錄組測序，利用Cutadap軟件對原始數據進行質量控制，利用HISAT2軟件將 clean reads比對到基因組上，使用FPKM 值作為基因表達量。

1.4 差異表達基因分析

利用DESeq2軟件進行基因的差異表達分析，設定閾值為log2|FoldChange|>0，校正P值(Padj)<0.05，當FoldChange≥1.5或≤0.05且P值<0.05時被定義為差異表達基因(DEGs)。

1.5 差異表達基因GO功能和KEGG通路富集分析

利用GO數據庫對篩選出的差異表達基因進行功能注釋，利用KEGG數據庫對差異表達基因的信號通路進行富集分析，篩選出與生長發育相關的功能基因。

1.6 實時熒光定量PCR驗證

隨機選取視黃醇結合蛋白7(RBP7)、肌鈣蛋白T2(TNNT2)、谷氨酸脫羧酶1(GADL1)、細胞表明糖蛋白(CD44)、肌肉生長抑制素(GDF8)、磷酸化酶激酶調節亞基β(PHKB)、熱休克蛋白H1(HSPH1)、RCSD結構域包含1(RCSD1)共8個差異表達基因進行實時熒光定量PCR驗證。根據GenBank中各基因序列，以GAPDH作為內參基因，利用Oligo 7.0軟件設計引物，引物信息見表1。引物均由南寧捷尼斯生物科技有限公司合成。使用實時熒光定量PCR驗證轉錄組測序結果。PCR反應體系10 μL：SYBR Green Ⅰ 5 μL，上、下游引物各0.4 μL，cDNA模板1 μL，ddH2O 3.2 μL。PCR反應程序：95 ℃預變性30 s；95 ℃變性5 s，60 ℃退火30 s，65 ℃延伸30 s，共45個循環。采用2-ΔΔCt法計算目的基因相對表達量。用SPSS 25.0軟件進行t檢驗，用GraphPad Prism 8.0軟件進行作圖。

表1 引物信息Table 1 Primer information

2 結 果

2.1 轉錄組測序數據分析

通過文庫構建和Illumina NovaseqTM6000平臺測序，原始數據經過質量控制后過濾后，6個樣本分別獲得34 296 198～41 647 642條clean reads，有效比對率在88.33%以上，Q20條含量在99.90%以上，Q30含量在97.51%以上，GC含量在50.00%～51.50%，說明測序數據可靠，可用于后續研究。

表2 轉錄組測序質控結果Table 2 Quality control results of transcriptome sequencing

2.2 差異表達基因分析

通過DESeq2軟件對1和60日齡廣西麻雞腿肌轉錄組測序數據進行分析，結果發現，與1日齡相比，60日齡共獲得差異表達基因2 304個，其中998個基因上調，1 306個基因下調(圖1)。

圖1 差異表達基因火山圖Fig.1 Volcano map of differentially expressed genes

2.3 差異表達基因GO功能和KEGG通路富集分析

2.3.1 GO功能富集分析 GO功能富集分析發現，差異表達基因GO富集的條目分為生物過程(biological process，BP)、細胞組分(cellular component，CC)和分子功能(molecular function，MF)3類，顯著富集的GO條目有572條，其中生物過程381條、細胞組分70條、分子功能121條，生物過程占大多數。由表3可知，生物過程中差異表達基因主要參與肌動蛋白成核的正調控、成肌細胞分化的負調控、成纖維細胞生長因子刺激的細胞反應、橫紋肌組織發育，差異表達基因的數目分別為5、7、5、3和3個；分子功能中差異表達基因主要參與肌動蛋白結合和膠原結合，差異表達基因數目分別為45和12個；細胞組分中差異表達基因主要參與細胞質基質、肌鈣蛋白復合物和肌球蛋白絲，表達差異基因數目分別為28、4和3個(表3)。

表3 部分生長發育相關GO條目Table 3 Partial growth and development related GO terms

續表

2.3.2 KEGG通路富集分析 KEGG通路富集分析發現，共獲得34個顯著富集的通路。由圖2可知，主要顯著富集的通路包括：氧化磷酸化、心肌收縮、Apelin信號通路、心肌細胞中的腎上腺素能信號轉導、MAPK信號通路、卡波西肉瘤相關皰疹病毒感染、脂肪酸降解、煙酸和煙酰胺代謝、晝夜夾帶、GnRH信號通路、多巴胺能突觸、胃酸分泌、糖尿病并發癥中的AGE-RAGE信號通路、血管平滑肌收縮、鞘脂代謝、糖酵解/糖異生、亨廷頓病、cAMP信號通路等。

圖2 差異表達基因KEGG通路富集分析Fig.2 KEGG pathway enrichment analysis of differentially expressed genes

2.3.3 生長發育相關基因篩選 通過GO功能和KEGG通路富集分析，獲得5個生長發育相關基因，分別為球蛋白重鏈10(MYH10)、肌球蛋白鏈15(MYH15)、成纖維細胞生長因子10(FGF10)、成纖維細胞生長因子16(FGF16)、GDF8基因(表4)，其中MYH10、MYH15和FGF10基因為上調表達差異基因，FGF16和GDF8基因為下調表達差異基因。

表4 生長發育相關基因篩選結果Table 4 Screening results of growth and development related genes

2.4 實時熒光定量PCR驗證

為了驗證轉錄組測序結果的可靠性，隨機選取8個在轉錄組測序中的差異表達基因進行實時熒光定量PCR驗證。由圖3可知，TNNT2、RBP7和RCSD1基因表達水平上調，GDF8、CD44、GADL1、PHKB和HSPH1基因表達水平下調，與轉錄組測序結果一致，證明轉錄組測序結果可靠。

圖3 實時熒光定量PCR驗證轉錄組測序結果Fig.3 Validation of RNA-Seq results by Real-time quantitative PCR

3 討 論

本研究對廣西麻雞腿肌組織進行轉錄組測序，共獲得2 304個差異表達基因，其中998個基因上調，1 306個基因下調，結合GO功能和KEGG通路富集分析對篩選出的差異表達基因進行挖掘，獲得MYH10、MYH15、FGF10、FGF16、GDF8共5個生長發育相關基因。He等[14]在生長快、慢的16周齡邊雞腿肌組織的差異轉錄本中篩選出108個差異表達基因，獲得了17個上調基因和91個下調基因。Yin等[15]在海陽黃雞12、16 d胚胎及1日齡和10周齡4個生長時期腿肌中篩選出6 150個差異表達基因，肌肉抑制素(MSTN)、肌源性分化1(MYOD1)、肌球蛋白重鏈6(MYF6)、肌球蛋白重鏈5(MYF5)和胰島素生長因子1(IGF1)等基因是雞生長發育的關鍵基因，這與本研究篩選的基因有所不同，說明不同生長時期所受到的基因調控會有所差異，但都可以為雞生長發育研究提供一定的參考。

MYH10和MYH15是肌球蛋白超家族的成員。MYH10具有多種功能，參與調節胞質分裂、細胞運動和細胞極性，在細胞骨架重組中起主要作用，其在小鼠神經系統和心臟的形成過程中至關重要[16]，在心臟發育和胚胎存活中起重要作用[17]。MYH15參與肌動蛋白結合和細胞骨架運動活性，Zhang等[18]研究3個品種雞肌肉生長發育的機制發現，與生長較慢的品種相比，生長較快的品種MYH15基因為顯著上調基因，還有研究發現MYH15在胚胎期肌肉中不存在，并且在出生后肌肉中可檢測到[19]。本研究發現，MYH15基因為顯著上調基因，說明MYH10基因對出生后雞的生長發育起關鍵調控作用。FGF10基因參與胚胎發育、細胞增殖和細胞分化，Rivetti等[20]發現，FGF10在小鼠中對誘導乳腺基板和白色脂肪組織的形成至關重要，FGF10基因對耳朵的誘導和生長也至關重要，FGF10基因突變小鼠在內耳分化過程中表現出輕微的缺陷[21]，并且可促進雞靜聲神經節神經突生長和神經元存活[22]。FGF16基因參與胚胎發育、細胞生長、形態發生，Ahsan等[23]研究表明，FGF16基因是雞生長發育的候選基因，推測其對雞的生長調控可能有重要作用，但是也有研究報道FGF9和FGF16基因可協同作為胚胎心肌細胞的生長因子[24]。GDF8基因具有骨骼肌生長負調節功能，有研究表明，GDF8基因可抑制肌源性衛星細胞的活化、增殖和分化，有研究發現GDF8對兩種不同骨骼肌的雞成肌衛星細胞的影響有所差異，胸大肌肌肉中分離的衛星細胞比股二頭肌的細胞中肌肉抑制素的抑制作用更敏感[25]，并且GDF8基因在雞胚胎的各種組織中廣泛表達[26]。本研究發現GDF8基因為顯著下調基因，進一步表明GDF8基因的表達有負調控作用，可影響雞肌肉的生長發育。

GO功能富集分析發現，差異表達基因顯著富集到了572個GO條目，其中包含了生長發育過程的非常重要的條目，如肌動蛋白成核的正調控、成肌細胞分化的負調控、成纖維細胞生長因子刺激的細胞反應、橫紋肌組織發育、肌動蛋白結合、膠原結合、細胞質基質、肌鈣蛋白復合物和肌球蛋白絲等，其中MYH10和MYH15基因富集在細胞質基質，FGF16和FGF10基因富集在成纖維細胞增殖的正調控，GDF8基因富集在成肌細胞分化的負調控，說明雞生長發育需要復雜的基因調控。KEGG通路富集分析發現，MYH10、MYH15、FGF10、FGF1基因富集在心肌收縮、MAPK信號通路、心肌細胞中的腎上腺素能信號傳導、緊密連接、肌動蛋白細胞骨架信號通路，其中心肌收縮、MAPK信號通路、心肌細胞中的腎上腺素能信號傳導富集靠前。MAPK信號通路對細胞增殖、分化、遷移、衰老和細胞凋亡起著至關重要作用[27]，另外緊密連接是細胞連接的重要組成部分，它構成離子和分子通過細胞的屏障，調節質膜頂端和基底蛋白和脂質的運動[28]，緊密連接參與控制細胞增殖和基因表達的作用[29]。本研究也發現了MAPK信號通路和緊密連接被顯著富集，說明MAPK信號通路和緊密連接對雞生長發育起到關鍵作用。

4 結 論

本試驗通過對廣西麻雞2個生長時期腿肌進行轉錄組分析，共獲得2 304個差異表達基因，其中998個基因上調，1 306個基因下調。通過GO功能和KEGG通路富集分析共獲得MYH10、MYH15、FGF10、FGF16和GDF8共5個生長發育關鍵基因。