基于转录组测序挖掘广西麻鸡生长发育相关基因

陈 婷，卢晟盛，杨丽丽，李添宝，张 敏，江明生

(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个生长发育关键基因。