基于转录组测序研究调控黄羽肉鸡干毛性状关键基因和信号通路

2024-01-17姬改革陈智武单艳菊刘一帆屠云洁邹剑敏章明巨晓军束婧婷张海涛唐燕飞蒋华莲

中国农业科学 2024年1期

姬改革，陈智武，单艳菊，刘一帆，屠云洁，邹剑敏，章明，巨晓军，束婧婷✉，张海涛，唐燕飞，蒋华莲

1 江苏省家禽科学研究所/江苏省家禽遗传育种重点实验室，江苏扬州 225125；2 广西金陵农牧集团有限公司，南宁 530000；3 江苏立华牧业股份有限公司，江苏常州 213168；4 广西富凤农牧集团有限公司，南宁 530024

0 引言

【研究意义】近年来，动物疫病和新冠疫情的爆发使中国的肉鸡消费结构发生了重大变化。在政府的鼓励和引导下，消费者购买黄羽肉鸡的消费习惯正由活鸡市场向生鲜市场逐步转变[1]。黄羽肉鸡是指含有中国地方品种鸡血缘的有色羽鸡，肉质普遍优于白羽肉鸡，是中国肉鸡养殖业的一个重要组成部分。与白羽肉鸡相比，处于上市日龄的黄羽肉鸡的羽毛拔出较为费力，拔毛后会有部分色素、毛根、甚至血液等残留在皮肤表面，严重影响冷鲜上市鸡的屠体外观。这主要与鸡体表覆盖的羽毛成熟度相关。羽毛成熟度也就是通常所说的干毛程度，是消费者用来间接判定肉鸡品质的依据之一，也能够影响上市活鸡的售价[2]。“干毛”是指鸡体表羽毛发育成熟时，羽毛根部毛管收缩变细变空，不含羽髓液的一种羽毛特征[2]。处于干毛状态的羽毛被拔出后，肉鸡屠体表面毛囊细小，无内容物残留，屠体表面看起来干净美观[3]。企业更倾向于培育干毛较早且一致性较好的肉鸡品系。【前人研究进展】羽毛是由毛囊发育而来，毛囊是羽毛生长发育的控制中心。成年鸟羽毛生长发育、脱毛、换羽再生等，都与毛囊周期发育相关。羽毛毛囊的生长周期大致可分为生长期和静止期[4]。处于静止期的毛囊，羽毛仍然完整保留在毛囊中，直到脱落进入下一个生长期。羽毛的再生受到季节、昼夜节律、激素和损伤（如拔毛）等多个因素的影响[5]。羽毛的再生是通过毛囊干细胞的增殖和分化驱动的。已知多个信号通路如FGF、Wnt/β-catenin、BMP 和Shh 等在羽毛生长发育和再生过程中具有重要作用[6]。研究发现Wnt信号是维持羽毛再生所必需的，存在于真皮乳头细胞的Dkk2 和BMP/TGF-β 等信号分子通过微调Wnt 信号的表达而调控羽毛的再生[7]。LIN 等[8]研究发现LncRNA 通过靶向抑制Wnt 通路相关基因表达调节羽毛再生和轴形成。通过转录组测序技术，发现Wnt，Shh，Notch，BMP 和细胞黏附分子之间的相互作用参与再生羽毛的形态发生[9-10]。这些信号分子的表达受到时间和空间的严格控制，相关基因表达紊乱会导致羽毛生长速度下降、羽毛质量差等问题[11]。色素生成信号通路基因在鸟类羽毛颜色图案多样性方面发挥了重要作用[12-13]。【本研究切入点】目前禽类羽毛发育相关的研究，大多数都集中在羽毛再生、羽毛颜色图案分子调控机理等方面。关于黄羽肉鸡羽毛毛囊周期调控，羽毛干毛性状的调控机制等研究较少。【拟解决的关键问题】本研究拟通过转录组测序技术，以干毛和未干毛羽毛毛囊皮肤组织为材料，初步筛选与羽毛干毛性状相关的重要候选基因和信号通路，为揭示羽毛毛囊周期调控分子机制和开发干毛性状相关分子标记辅助黄羽肉鸡育种提供基础资料。

1.材料与方法

1.1 样品制备

试验于2020 年6—9 月进行，试验地位于江苏省家禽科学研究所仪征试验基地。实验鸡来自江苏立华牧业股份有限公司、江苏省家禽科学研究所等单位联合培育的配套系花山鸡。在相同的饲养条件下，采用标准化程序饲喂和免疫，自由饮水。观察鸡羽毛生长状况，发现大多数鸡的背羽在5 周龄时呈现出未干毛的特征，9 周龄时呈现干毛特征。因此，分别在5 周龄和9 周龄各选择3 只背羽发育一致性较好的母鸡作为试验材料。将鸡放血屠宰后去毛。以背中线为轴，采集左侧皮肤组织2 cm×2 cm 大小的皮肤样迅速放入液氮后转入-80 ℃冰箱保存备用。采集右侧同一部位的皮肤，用10%甲醛溶液固定，固定后的皮肤组织按照文献[14]中的方法进行脱水、石蜡包埋、苏木精和伊红（Hematoxylin and eosin，H.E）染色并拍照观察。

1.2 总RNA 提取与转录组测序

根据试剂盒的说明，使用TRIzol 试剂（Invitrogen,USA）从每个皮肤样本中提取总RNA。分别用NanoDrop 2000（Thermo，美国）和Agilent 2100 Bioanalyzer（Agilent Technologies，CA，美国）测试RNA 样品的浓度和完整性。质检合格后，每个样品取1 μg 总RNA，进行文库构建，共构建6 个文库，利用北京百迈客生物科技公司Illumina 平台进行测序。测序获得的原始数据（raw reads）经过滤，去除接头、含有ploy-N的等低质量碱基序列，获得质量较高的序列（clean reads）进行下一步分析。原始数据已经上传至NCBI SRA 数据库，存储号为PRJNA739526。将得到的clean reads 用Hisat2（v2.1.0）软件比对到鸡参考基因组（Galgal 6.0）上。根据比对到基因组上的位置信息，统计每个转录本的read count 信息，采用FPKM 值对基因表达水平进行定量[15]。

1.3 差异表达基因数据富集分析

差异表达基因（DEGs）分析采用DESeq2 (v1.30.0)软件进行。以Pvalue＜0.05 和|fold-change|＞2 为差异基因的筛选标准。获得的差异表达基因利用topGO 软件[16]进行 GO（gene ontology）富集分析。使用clusterProfiler 软件[17]来进行差异表达基因 KEGG（Kyoto Encyclopedia of Genes and Genomes）的通路显著性富集分析。

1.4 基因共表达网络分析和hub 基因筛选

采用在线STRING 软件（http://string-db.org/）对筛选到的差异基因进行分析，将互作评分大于0.9 的蛋白关系对，导入cytoscape3.3 软件绘制蛋白相互作用网络图。利用Cytoscape 中的CytoHubba 应用程序，通过4 种分析算法（MCC、DMNC、MNC 和Degree）计算网络中的hub 基因[18]。

1.5 GSEA 基因富集分析

GSEA 基因富集分析是对两组样本中的所有基因进行输入，找出具有协同差异的基因集，不需要指定阈值（P或FDR值）来筛选差异基因，兼顾了GO/KEGG富集分析中容易遗漏掉的差异表达不显著却有着重要生物学意义的基因。利用GSEA v3.0（http://www.broadinstitute.org/gsea/index.jsp）软件对所有基因进行分析。根据标准化富集得分（normalized enrichment score，NES）、错误发现率（false discovery rate，FDR）等对富集通路进行筛选分析，主要针对KEGG 代谢和信号转导途径进行分析。

1.6 差异表达基因实时荧光定量PCR 验证（RT-qPCR）

根据NCBI 查询到的基因序列，设计引物，由上海生工（生物）技术有限公司合成。采用Trizol 试剂盒提取组织RNA，反转录为cDNA，采用SYBR Green I 法，以β-actin 为内参进行荧光定量PCR 验证。反应体系20 μL：10 μL qPCR Master Mix，上下游引物各0.4 μL（10 μmol·L-1），cDNA 模板2 μL，50×ROX Reference Dye 0.4 μL，加ddH2O 补足20 μL。反应条件按照试剂盒（ChamQ SYBR Color qPCR Master Mix，诺唯赞）说明书推荐的条件设置PCR 扩增条件，并采集溶解曲线。以5 周龄皮肤组织为对照，采用2-ΔΔCt法计算目的基因表达量。引物具体信息见表1。

表1 引物序列表Table 1 The qPCR primers were used in this study

2 结果

2.1 不同发育时期羽毛毛囊组织形态学比较

不同发育时期羽毛形态变化如图1 所示，当羽毛处于未成熟状态时，羽轴根部颜色乌黑，根部羽髓液充足。随着肉鸡日龄的逐渐增加，至9 周龄时羽毛逐渐发育成熟，羽毛根部逐渐中空变细，颜色变浅，羽髓已被吸收，羽毛鞘逐渐角质化，羽毛已经从羽鞘中突破出来，也就是通常所说的干毛。

图1 鸡背部羽毛不同发育时期形态变化Fig.1 Morphological changes of chicken back feathers at different developmental periods

毛囊组织H.E 染色结果如图2 所示，在未干毛的羽毛毛囊组织中（图2-A），基底部真皮乳头（dermal papilla）、真皮髓质（dermal pulp）在毛囊中清晰可见，能够不断为羽毛生长提供营养。在干毛的羽毛毛囊组织中（图2-B），真皮乳头已经萎缩；髓质已经退化至底部，几乎看不到明显的髓质。

图2 鸡背部羽毛毛囊不同发育时期组织切片H.E 染色Fig.2 Histological H.E staining of chicken back feather follicles at different developmental stages

2.2 不同发育时期毛囊皮肤组织转录组测序数据统计

不同发育时期背部皮肤组织样本进行转录组测序。未干毛皮肤（H1 组，生长期毛囊）3 个样本共获得136 001 102 个raw reads；干毛皮肤（H2 组，静止期毛囊）3 个样本共获得132 540 460 个raw reads。进行质控过滤后两组分别获得 127 225 676 个和124 836 360 个clean reads。两组Q30 数据都在92%以上，GC 含量为49%左右。与参考基因组比对结果显示，生长期和静止期两组总比对率都高于94%，其中2%左右的序列比对到基因组多个位置，92%左右的序列比对到基因组唯一位置。说明测序获得的数据质量可靠，可信度高，可以用于后续分析。相关信息具体见表2。

表2 测序信息统计及与参考基因组比对Table 2 Summary of sequencing reads statistics and mapping to the reference genome

2.3 差异表达基因功能富集分析

以|foldchange|＞2 和P＜0.05 作为标准进行差异表达基因筛选，以生长期为对照组，共筛选到942 个DEGs，其中有558 个基因在静止期皮肤组织表达量上调，384 个基因表达量出现下调（图3）。与羽毛发育相关的基因如WNT6、FOXO3、FZD4、LRP6、LEF1、BMPR2、EDA2R等被筛选到。对差异表达基因进行GO 功能富集分析，分别显著富集到426 个生物学过程，81 个细胞组分和95 个分子功能（P＜0.05）。前30 个被显著富集的生物学过程如图4-A 所示，DEGs主要富集到细胞分裂、有丝分裂细胞周期的调节、有丝分裂染色体凝集、有丝分裂姐妹染色单体分离、DNA 几何变化和DNA 复制等有关的生物学过程中。Wnt 信号通路、胰岛素、骨骼肌发育、脂肪酸代谢和运输相关的生物过程也被显著富集（P＜0.05）。细胞组分富集到中央纺锤体复合体、浆膜、纺锤体微管等部位。分子功能富集到ATP 结合、DNA 复制原点结合、微管结合、肌动蛋白结合等方面。

图3 不同发育时期毛囊（生长期vs 静止期）皮肤组织差异表达基因火山图（|foldchange|＞2 和 P＜0.05）Fig.3 Volcano plots of DEGs between feather follicles at different developmental periods (growth and rest stage)in chicken skin samples

图4 差异表达基因前30 个显著富集GO 生物学过程（A）和KEGG 富集通路（B）Fig.4 The top 30 significantly enriched GO biological processes (A) and KEGG enrichment pathways (B) for differentially expressed genes

KEGG 富集结果如图4-B 所示，共有13 个信号通路被显著富集（P＜0.05），最显著富集的通路包括细胞周期，卵母细胞减数分裂，DNA 复制和同源重组都与细胞周期活动有关（P＜0.05）。调控细胞衰老、凋亡的p53 信号通路，与毛囊生长发育有关的MAPK、TGF-β 信号通路也被显著富集到（P＜0.05）。

2.4 差异表达基因共表达互作网络（PPI）构建与hub基因筛选

利用string 数据库进行差异表达基因共表达网络构建，并通过Cytoscape 对结果进行展示，获得202个节点和916 条边组成的网络（图5）。进一步采用CytoHubba 插件中的MCC、MNC、Degree 和EPC 算法获得得分较高的前10 个hub 基因（表3）。采用Venn对4 种算法获得的前10 个基因进行分析，发现有6个基因CDK1、MAD2L1、BUB1、CCNB2、PLK1及BUB1B在4 种算法中都能找到，说明这6 个基因在蛋白互作网络中占有重要位置，是重要的hub 基因。

图5 差异表达基因蛋白互作网络Fig.5 Differentially expressed gene protein-protein interaction network

表3 通过CytoHubba 插件4 种算法在PPI 网络中排名前10 的基因Table 3 Top 10 genes identified in the protein-protein interaction (PPI) network by CytoHubba plug in 4 algorithms

2.5 GSEA 分析

进一步采用GSEA 分析对生长期和静止期皮肤毛囊样本中的所有基因进行分析，结果显示，5 个通路包括DNA 复制、细胞周期、嘧啶代谢、同源重组、氧化磷酸化在生长期毛囊皮肤样本中显著富集（|NES|＞1,FDR＜0.25）。其中，DNA 复制和细胞周期通路被极显著富集（图6）（FDR＜0.01），这2 个通路富集的基因表达在生长期毛囊样本中上调。在DNA 复制基因组通路中，MCM2、DNA2、FEN1、PRIM2、MCM3等被确定为核心富集基因。在细胞周期通路中，PPI 网络中的枢纽基因如CDK1、BUB1和BUB1B都被确定为核心富集基因。同时富集到黏着斑、紧密连接、FOXO、MAPK、胰岛素、促性腺激素释放激素（GNRH）和TGF-β等30 个通路在静止期毛囊皮肤组织中表达量上调（|NES|＞1,FDR＜0.25）。其中MAPK 和TGF-β信号通路在差异表达基因KEGG 富集分析中也被筛选到（图7）。

图6 GSEA 富集在生长期毛囊皮肤组织中上调的DNA Replication（A）和Cell Cycle（B）通路Fig.6 GSEA enrichment of DNA replication (A) and cell cycle (B) pathways up-regulated in feather follicle skin tissue at the growth phase

图7 GSEA 富集在静止期毛囊皮肤组织中上调的MAPK（A）和TGF-β（B）通路Fig.7 GSEA enrichment of MAPK (A) and TGF-β (B) pathways up-regulated in feather follicle skin tissue at the resting phase

2.6 差异表达基因实时荧光定量PCR 验证

分别以处于生长期和静止期的毛囊皮肤组织各8 个为材料，从富集的信号通路中随机选择与细胞周期、毛囊生长发育等相关的基因进行荧光定量PCR 检测，计算log2（Fold change, FC）值，并与转录组测序结果比较，如图8 所示，以生长期为对照，BUB1、BIRC5、CDK1、NEK2、CDC6基因在静止期皮肤组织中表达量都是下调的；ANXA1、FZD4、BCL2 和EDA2R 基因在静止期皮肤组织中表达量都是上调的。8 个DEGs 的荧光定量PCR 检测结果与转录组测序结果变化趋势基本一致。

图8 差异表达基因荧光定量PCR 验证Fig.8 RT-qPCR validation of differentially expressed genes

3 讨论

3.1 羽毛干毛性状与毛囊生长周期的关系

毛囊的形态结构十分复杂，控制着羽毛生长速度、颜色等性状，同时具有独特的生长特性，能够进行周期循环。本研究对不同时期花山鸡羽毛发育情况进行观察，发现5 周龄时背部羽毛毛根粗大，毛管内充满髓液，处于未干毛状态。9 周龄时羽毛根部收缩变细，无羽髓液，符合干毛形态特征。一个完整的毛囊主要包含表皮和真皮组织。真皮髓质位于毛囊中央，富含神经、血管及色素，当羽毛进入生长期，髓质随着血管通道的扩张而增大，为羽毛生长提供营养。真皮乳头是永久性结构，位于毛囊的基底部，随着毛囊的生长周期性的增大或萎缩。本研究结果显示，未干毛和干毛羽毛毛囊组织学形态变化基本符合毛囊生长期和静止期的特征，与LIN 等在羽毛再生中的研究结果是一致的[3,19]。说明羽毛的生长发育受到毛囊生长周期调控。

3.2 细胞周期通路参与毛囊生长发育周期调控

对处于生长期和静止期毛囊皮肤组织进行转录组测序，共检测到了942 个差异表达基因。羽毛毛囊不是独立存在的，能够与周围的血管、肌肉、神经和脂肪等组织联合，共同为生长中的羽毛提供养分。因此，差异表达基因富集分析结果显示，除了常见的毛囊生长相关通路，肌肉发育、脂质代谢、血液循环等生物学过程也被富集。CHEN 等[5]报道，羽毛毛囊的真皮乳头和表皮细胞存在性激素受体，鸟类羽毛周期的转换受到性激素的调控。胰岛素和GNRH 信号通路在GSEA 结果中被富集到，说明胰岛素和性激素在毛囊周期调控中具有重要作用。WANG 等[20]研究认为，免疫细胞在人和动物毛发周期调控中具有重要作用。本研究结果显示，淋巴结形态发生通路也被富集到。TAN 等[21]研究报道，VEGF 和VEGFR-2 基因在黄羽肉鸡干毛和未干毛皮肤组织中表达差异显著。本试验差异表达基因集中未筛选到这2 个基因，可能与所用黄羽肉鸡品种及日龄等不同有关。GO 和KEGG 同时显著富集到众多与细胞增殖和凋亡相关的生物学信号通路。如p53 信号通路具有细胞凋亡、DNA 修复、细胞周期阻滞、细胞衰老等多种生物学效应[22]。

构建差异基因蛋白互作网络，采用CytoHubba 插件识别PPI 网络中的hub 基因，发现CDK1、MAD2L1、BUB1、CCNB2、PLK1 和BUB1B 在网络中占有重要地位。MAD2L1 和BUB1B 是关键的有丝分裂检查点基因，当染色体脱离有丝分裂纺锤体时，该检查点保证了染色体的及时和正确分离[22-23]。PLK1 属于PLK激酶家族成员，调节有丝分裂中心体成熟和分离，纺锤体组装形成，染色体排列和分离，以及细胞分裂和脱落，在细胞周期调节、细胞增殖以及抑制细胞凋亡等方面具有重要作用[24]。BUB1 能够招募PLK1 到动粒中以稳定动粒-微管连接，从而确保染色体准确分离[25]。CDK1 是驱动有丝分裂重塑的主要激酶[26]，通常在G2/M 期表达量达到高峰，在G1/S 期表达量显著下调，参与对细胞存活至关重要的多种生物学过程，包括G2/M 转换，检查点激活，DNA 修复和DNA 复制[27]。WU 等[28]在鸡羽毛毛囊再生中的研究证实，从生长期至静止期CDK1 基因表达量逐渐下降，同时伴随着髓质和真皮乳头细胞数量的降低。CCNB2 是一种细胞周期调节因子，可与CDK1 形成复合物，调节细胞周期的G2/M 期，在有丝分裂的启动中起重要作用[29]。CCNB2 在G1 期合成，然后在早期、中期和晚期迅速减少[30]；CCNB2 的异常表达会导致G2/M 检查点功能在细胞周期中失效，使细胞增殖减少[31-32]。CDK1、BUB1和BUB1B等是GSEA 结果中显著富集通路的核心基因。在鹅羽毛再生中的研究发现，与细胞周期相关的基因如CCNB3、CCNY、CDK3 等也被筛选到[33]。说明CDK1、BUB1和BUB1B等与细胞周期转换相关的基因在调控鸡羽毛生长发育中具有重要作用。

GSEA 结果与差异表达基因富集结果基本一致，与细胞增殖相关的DNA 复制和细胞周期通路被显著富集。这可能是因为毛囊表皮和真皮干细胞周期性的激活和休眠在毛囊生长周期的转换中发挥着重要作用。在起始阶段，毛囊干细胞在周围信号分子的刺激下恢复生长，进入生长期，真皮乳头周围的上皮层持续增厚，上皮基底层细胞增殖旺盛，分枝区细胞增殖、分化为羽枝和羽小枝；转为静止期，真皮乳头和毛囊宽度都会萎缩，髓质间充质细胞程序性死亡，并向羽毛近端退去，上皮基底层细胞增殖大大减少[10]。整个生长周期伴随着大量细胞的增殖、分化和凋亡。DNA复制和细胞周期通路的相关基因在生长期毛囊皮肤组织中表达量都是上调，与预期结果基本一致。说明细胞周期调控相关基因在毛囊周期转换的调控中发挥了重要作用。

3.3 毛囊生长发育相关信号分子在毛囊周期调控中的作用

毛囊的形态发生是胚胎期表皮和真皮相互作用的结果。成年后毛囊干细胞周期性的激活和休眠是由表皮和真皮细胞间的微环境所决定的[34]。在小鼠中的研究认为毛囊进入生长期主要受Wnt/β-catenin 信号传导的调控，而静息期主要受TGF-β/BMP/Smad 活性调控[35]。采用转录组测序技术比较羊不同生长周期的皮肤组织[36]，发现BMP、TGF-β、MAPK 信号通路在调节羊毛囊周期生长的过程中具有重要作用。本研究结果发现，差异表达基因被显著富集到Wnt、MAPK和TGF-β 信号通路。其中，MAPK 和TGF-β 信号通路在GSEA 结果中也被显著富集到。说明这3 个信号通路在鸡羽毛周期发育过程中起重要作用。Wnt家族成员众多，Wnt 配体及其下游信号分子在毛囊的定位、形态发生，周期再生及后期羽毛分枝形成过程中都发挥了重要作用[6,37-38]。研究表明[5]，羽毛再生过程中，毛囊真皮干细胞的周期性激活，是由Wnt信号介导的。DKK 是经典的Wnt 信号通路抑制剂。在鸡羽毛毛囊中过表达DKK2，能够使毛囊中细胞增殖减少，细胞凋亡增加，最终导致羽毛再生延迟[7]。SONG 等[39]在鹅胚中的研究认为，TGF-β1 作为一种功能性细胞因子，能够抑制真皮成纤维细胞增殖，促进真皮成纤维细胞凋亡。BMPs 是分泌型信号分子转化生长因子（TGF）-β 超家族的成员，通过与其受体结合发挥生物学作用。研究发现，BMP2 在牛静止期毛囊皮肤中高表达[40]，能够抑制毛囊进入生长期。BMPR2 是BMPs 受体，推测BMPR2 基因在鸡静止期毛囊皮肤组织中表达升高能够抑制羽毛毛囊向生长期转变。

FANG 等[41]研究发现，鸡不同长速羽毛皮肤组织差异表达基因主要富集在Wnt、TGF-β 和MAPK 等信号通路。比较不同发育时期鸭颈部皮肤组织差异表达miRNA[42]，发现高表达差异miRNA 大多与细胞周期调控有关，预测的靶基因富集在黏着斑、MAPK、Wnt、Notch 和TGF-β 等信号通路。目前关于TGF-β 和MAPK信号通路在家禽羽毛毛囊周期转换中的作用机制研究还相对较少。在鼠上的研究发现[43]，MAPK 信号通路的激活能够调控毛囊干细胞进入静止期。SU等[44]研究认为，MAPK 信号通路能够通过调控细胞的增殖和分化来调节哺乳动物毛囊生长周期转换。冯云奎等[45]研究证实，miRNA 通过靶向MAPK 信号通路相关基因表达调控羊毛囊干细胞的增殖和凋亡。MAPK 被认为是BMPs 的下游信号通路之一[46]。提示TGF-β/BMP/MAPK 信号通路基因可能通过调控细胞周期相关基因表达，进而影响鸡羽毛毛囊生长周期的转换。