天台地区黄牛常年以山区放牧为主，为适应山区的环境，天台黄牛体型逐渐小型化，经过长期的自然和人工选择最终形成了微小体型的黄牛[1]。始丰微型牛是浙江省地方品种天台黄牛的一个分支，也是我国唯一的微型牛类群。动物体型与基因、营养和环境等因素相关。Nicholls等[2]通过体外饲养控制蜜蜂幼虫的饮食质量，确定了蛋白质和碳水化合物对蜜蜂成年存活时间、体重和大小存在显著影响。方冰娃等[3]利用SLAF-Seq技术研究与塔里木兔体型大小相关的受选择基因，发现磷酸酶与张力蛋白同源物基因（PTEN）、诱导信号通道蛋白1（WNT1）、原癌基因wnt7a蛋白抗体（WNT7A）、成纤维细胞生长因子13（FGF13）、核糖体蛋白S6激酶多肽1（RPS6KB1）等候选基因的功能与动物的体型大小相关，并且其中部分基因富集在调控体型相关的通路中。转录组测序技术（RNA-Seq）是转录组水平上分析差异基因表达mRNA可变剪切的重要手段[4]。贺花[5]利用高通量转录组测序技术筛选出6 800个在牛肌肉发育过程中差异表达的基因，通过GO功能分析发现，其中大部分条目涉及生长发育过程。Wickramasinghe[6]使用RNA-Seq技术筛选出奶牛两个泌乳期分别有19 094个和18 070个差异表达的基因。李建平等[7]对育肥猪进行了转录组测序并筛选出938个差异表达的基因，差异表达最显著的基因集中在萜类骨架生物合成途径和酮体形成降解途径。因此，RNA-Seq技术揭示动物自身的遗传机制，发掘关键基因及探索其调控机制[8]。本研究应用转录组测序技术分析普通天台黄牛和始丰微型牛在生长发育相关基因上的差异，探讨始丰微型牛微小体型的遗传机制。1材料与方法1.1样品采集从天台黄牛保种场选取12月龄健康天台黄牛15头。母牛组中始丰微型母牛6头，普通天台黄牛母牛4头；公牛组普通天台黄牛2头，始丰微型公牛3头。采集黄牛耳组织样放入RNA保护剂管，-80 ℃保存。1.2总RNA提取、浓度测定和质量检测提取总RNA前，对全部相关器械进行去RNA酶处理。提取牛耳组织的总RNA时，根据Trizol法说明进行提取。应用Nanodrop2000软件检测RNA的纯度与浓度。1.3上机文库制备判定RNA样品符合规格后，通过Oligo（dT）磁珠进行样本mRNA的富集，采用相同方法处理totalRNA；DNA探针与rRNA进行杂交，并结合核糖核酸酶（RNaseH）选择性消解DNA/RNA杂交链，之后用脱氧核糖核酸酶I（DNaseI）去除DNA探针，进行提纯处理后即可获得所需RNA。通过buffer对所需的RNA进行片段化处理，随机对N6引物反转录，使cDNA二链构成双链DNA。此外，补平双链DNA末端，并在5'端磷酸化，3'端构成突触的“A”末端，使其连接3'端带“T”的接头。设计特异引物对连接产物进行PCR扩增。热变作用下PCR产物会转变为单链DNA，利用段桥式引物环化单链DNA后可得到单链环状的DNA文库，此时可上机测序。1.4生物信息分析测序的原始数据使用SOAPnuke（v1.4.0）软件进行过滤，过滤原则遵循以下3点：接头污染的reads去除；低质量reads去除，即质量低于10碱基占该reads总碱基数的比例超过20%的reads；未知碱基N含量超过5%的reads去除。质控合格后，使用Hisat（v2.1.0）工具对RNA-seq reads进行参考基因组比对，通过统计reads在参考序列上的分布情况和比对率等，对比结果进行二次质控，再应用Bowtie2（v2.2.5）工具对参考基因序列与clean reads进行对比，进而得到全部基因表达量。使用DEGseq方法对两组表达的基因进行差异显著性分析，使用FDR法矫正P值为Q值，判断基因是否为差异显著基因时，遵循原则是log2FC≥2且Q值≤0.001。对差异表达基因进行生物通路分类时，参照KEGG注释结果应用R软件的phyper函数对差异表达基因做富集分析，先对P值进行求解，再通过FDR校正P值，若FDR≤0.01表示显著富集；最后对蛋白互作网络进行分析。2结果与分析2.1质控后Clean reads的统计结果（见表1）由表1可知，过滤后的reads中质量值20碱基数占比超过97%，过滤后的reads中质量值30碱基数占比超过89%，过滤后的reads占比超过99%，表明测序质量高。10.13557/j.cnki.issn1002-2813.2023.20.018.T001表1质控后Clean reads的统计结果组别过滤后的碱基总数/Gb过滤后的reads数/M过滤前的reads数/M碱基数占比/%过滤后的reads占比/%过滤后的reads中质量值30过滤后的reads中质量值20普通体型母牛comf11.0821.7521.7790.0497.9099.87comf21.0821.8621.8890.9198.3499.88comf31.0821.8621.8891.1398.4099.86comf41.0821.7221.7990.2998.0999.68始丰微型牛母牛minif11.0821.7221.7889.9198.0199.70minif21.0821.8321.8989.9698.0299.72minif31.0821.8121.8490.8098.1499.87minif41.0821.8621.8991.0298.3599.88minif51.0821.8621.8990.8498.3399.84minif61.0821.8421.8891.2798.4599.81普通体型公牛comm11.0821.8521.8891.0298.3999.84comm21.0821.8721.8991.2198.4399.86始丰微型牛公牛minim11.0821.7221.8290.0898.0299.66minim21.0821.8621.8991.1398.4499.87minim31.0821.8421.8991.0898.3799.782.2基因定量分析与测序质量验证（见图1）由图1可知，随机性分析5'端低、3'端高，主要由mRNA的降解引起，表明测序结果符合要求。10.13557/j.cnki.issn1002-2813.2023.20.018.F001图1reads在转录本上的分布2.3DEGseq差异表达基因分析（见表2）母牛组的下调基因与上调基因分别为165、298个，共计463个差异表达基因；公牛组的下调基因与上调基因分别为355、262个，共计617个差异表达基因。公牛组与母牛组差异倍数前10位的基因见表2。10.13557/j.cnki.issn1002-2813.2023.20.018.T002表2公牛组与母牛组差异倍数前10位的基因组别基因ID基因名称基因符号Log2（minif/comif）P值表达状态母牛组（comf-minif）510798长非编码RNA LOC510798LOC5107981.007 6892.989 330上调518050膜联蛋白A3ANXA31.013 0854.393 490上调784304胞苷单磷酸激酶2CMPK21.014 6019.033 010上调783907血管生成素2ANG21.020 9101.051 418上调514113尿酸氧化酶UOX1.020 9244.012 548上调516929半胱氨酸-丝氨酸蛋白CSRNP2-8.072 0543.631 863下调281357利钠肽受体2NPR2-7.072 8502.223 715下调282259胰岛素样生长因子结合蛋白1IGFBP1-6.650 0442.503 275下调505809叶酰多聚谷氨酸合成酶FPGS-6.467 8483.207 218下调521556体节发生蛋白MSGN1-5.615 4862.716 079下调公牛组（comm-minim）518186植物同源域指蛋白6PHF64.979 6620.000 004上调538673谷氨酸转运体基因SLC1A61.629 0120.000 822上调280828趋化因子配体8CXCL81.741 4860.000 019上调281308基质金属蛋白酶1MMP11.103 6200.000 021上调503553滋养层结构域蛋白2TKDP23.174 4460.000 020上调518761螺旋型肌动蛋白成核因子2SPIRE2-8.967 0794.219 514下调281348神经丝轻链NEFL-4.750 3660.000 614下调505809叶酰多聚谷氨酸合成酶FPGS-6.664 9122.692 280下调531420糖蛋白2GP2-6.467 8483.207 218下调517420胰高血糖素样肽1受体GLP1R-1.580 4410.000 313下调2.4公牛、母牛基因KEGG通路分类（见图2、图3）由图2、图3中KEGG通路可知，公牛、母牛差异显著基因主要集中在6类代谢通路中，分别是细胞加工、环境信息处理和遗传信息处理等。10.13557/j.cnki.issn1002-2813.2023.20.018.F002图2母牛组KEGG通路分类10.13557/j.cnki.issn1002-2813.2023.20.018.F003图3公牛组KEGG通路分类2.5公牛、母牛差异表达基因KEGG通路富集分析及候选基因筛选（见表3）通过R软件对差异表达基因进行KEGG通路的富集分析，发现天台牛样本的差异显著表达基因（P0.05）共富集到42条通路中，其中与生长发育相关的通路共10条。针对转录组差异表达基因，通过对其KEGG富集分析并结合相关研究文献内容，在KEGG通路中选出影响生长发育的基因。研究发现，Ras信号通路、TGF-β信号通路、MAPK信号通路、Rap1信号通路等KEGG通路中富集了较多可影响生长发育的基因，并从中筛选出多个较为重要的差异表达基因，包括TGFB2、IGFBP3、DLK1、BMP15、Rap1、IGF-1、AMPK、STAT3、CCNA1、EGF、FGF、cAMP等。10.13557/j.cnki.issn1002-2813.2023.20.018.T003表310条与动物发育相关的KEGG通路信息通路ID名称Ko04014Ras信号通路（Ras signaling Pathway）Ko04152AMPK信号通路（AMPK signaling Pathway）Ko04024CAMP信号通路（CAMP signaling Pathway）Ko04015Rap1信号通路（Rap1 signaling Pathway）Ko04022cGMP-PKG信号通路（cGMP-PKG signaling Pathway）Ko04013MAPK信号通路（MAPK signaling Pathway）Ko04310Wnt信号通路（Wnt signaling Pathway）Ko04350TGF-β信号通路（TGF-β signaling Pathway）Ko04115P53信号通路（P53 signaling Pathway）Ko04380破骨细胞分化（Osteoclast differentiation）2.6差异表达基因的蛋白网络（见图4）由图4可知，可知差异基因所表达的蛋白质之间存在相互关联。10.13557/j.cnki.issn1002-2813.2023.20.018.F004图4差异基因蛋白网络3讨论近年来，RNA-Seq技术成为在转录组水平上分析mRNA可变剪切和差异基因表达的不可缺少的工具。本研究主要利用RNA-seq技术对天台黄牛和始丰微型牛之间的差异表达基因进行了检测，并对差异表达基因进行KEGG通路富集分析，筛选出一些与生长发育相关的重要差异表达基因，如DLK1、Rap1、IGF-1、AMPK、CCNA1、TGFB2、IGFBP3等，可作为影响始丰微型牛体型大小的关键候选基因。DLK1蛋白基因在细胞分化和组织发育中起重要作用[9]。Chen等[10]发现，DLK1在胚胎发育时期的软骨细胞与骨骼细胞中高度表达。贾焕珍等[11]采用免疫组织化学染色检测肿瘤标本中DLK1、GH表达情况，结果表明，DLK1可增加血清GH水平，抑制GH3细胞增殖，对肌肉生长起到促进作用。Rap1是小分子G蛋白Ras家族成员之一，在细胞分化、增殖和附着方面起到重要的调节作用，对神经突触生长也具有促进作用[12-13]。Liu等[14]研究显示，Rap鸟嘌呤核苷酸交换因子3（RAPGEF3）是GTP酶Rap1和Rap2的鸟嘌呤核苷酸互换因子，可增强Rap1的活性，促进细胞附着和轴突再生，并通过与神经生长因子相互作用促进轴突延伸。AMPK被称为“细胞能量调节器”，是调节细胞能量代谢的开关，可调控ATP的分解与合成，对能量变化高度敏感[15-16]。AMPK与沉默信息调节因子1（Sirt1）是调节糖脂代谢的重要因子，可抑制脂质积累并提高肌肉组织的线粒体功能[17]。细胞周期调节蛋白A1（Cyclin A1，Ccna1）可促进细胞进入分裂期并且能够调节DNA合成，对细胞增殖和骨骼肌生长具有重要的调节作用[18-19]。Pakula等[20]研究表明，Ccna1表达改变可能导致面肩肱型肌营养不良症（FSHD）。胰岛素样生长因子（IGFs）可调节糖、蛋白质及脂肪的合成，还能够控制调节细胞的分化与增殖情况，故IGFs在动物生长过程具有重要作用[21]。刘子嶷等[22]研究IGF-1对体外培养猪卵母细胞成熟率的影响，结果表明，IGF-1可显著影响卵泡卵丘扩展率和第一极体排出率，改善卵母细胞的氧化还原状态。丁艳等[23]对撒坝猪IGF-1基因进行研究，并指出该基因与70~120 d日增重、180 d体重存在显著相关性，与胸围及体长等也存在显著相关，与管围不存在相关性。TGFB2属于的细胞膜受体，在大多数细胞和组织中均有表达[24]，并具有促进细胞外基质合成以及细胞增殖等作用[25]。在骨骼发育方面，TGFB2能够使骨膜间充质细胞加速分化与繁殖，使骨细胞加速繁殖，促使骨桥蛋白与Ⅰ型胶原组合。Anna等[26]研究指出，小鼠四肢中敲除该基因后，TGF-β传导信号中断，使四肢的指骨关节出现问题。通过深入研究证实，TGFB2基因会对Noggin等的基因表达进行调节，从而造成软骨细胞出现异常现象。Serra等[27]研究表明，在将小鼠的TGFB2基因敲除以后，小鼠头盖骨会出现发育不全的情况，造成胚胎死亡。IGFBP家族有7个成员，即IGFBP1~IGFBP7。IGFBP3作为IGFs的载体在循环和调节IGF活性方面起重要作用，具有抗肿瘤作用，并且与神经内分泌的分化、迁移和侵袭有关[28-29]。王玥[30]通过过表达和敲除IGFBP3后发现，参与猪脂肪细胞合成的一些基因如C/EBPα、C/EBPβ、PPARα、PPARγ、FASN等会出现表达量上升或下降，表明IGFBP3基因可间接通过调控一系列基因转录调控参与猪脂肪细胞分化过程，在脂肪生成过程中起到重要作用。4结论本研究发现，在始丰微型牛中，P53信号通路中富集的IGFBP3和TGFB2处于下调状态，抑制了胰岛素生长因子和转化生长因子的表达，进一步抑制了甲状腺激素与生长激素，诱使始丰微型牛发育缓慢，体形偏小。