麦洼牦牛不同生长阶段肝脏差异表达基因分析

傅芳 王利* 字向东

（1 西南民族大学，青藏高原动物遗传资源保护与利用教育部和四川省重点实验室，成都 610041）

（2 西南民族大学，动物科学国家民委重点实验室，成都 610041）

肝脏是哺乳动物机体内最大的消化腺及代谢器官，利用内分泌因子不断调节新陈代谢平衡，在各种营养或有害物质代谢、蛋白质合成和消化液合成等过程中发挥重要功能(Schlegelet al.,2012; Shappellet al., 2019; Jalaliet al., 2020)。动物生长发育过程中，外源性生物活性物质和营养物质对肝脏的成熟和代谢进行调节(徐函等，2018)。已有研究探讨了Hippo 信号通路参与调控新生仔猪肝脏的生长代谢(尹聪，2013)，并利用RNA-Seq 技术分别解析了不同发育阶段猪(贺伸，2017)、鸡(Liet al.,2015)等畜禽肝脏组织的生长发育机制，以理解其生长发育机制和规律，但不同发育阶段牦牛肝脏组织的生长发育机制及其基因的表达至今未见研究报道。因此，研究高原动物肝脏的生长发育机制有助于了解高原哺乳动物肝脏生长发育进化进程及遗传基础，为牦牛肝脏生长发育调控提供理论依据。

牦牛(Bos grunniens)分布于青藏高原，能适应高海拔高寒低氧的气候环境条件，可为人类提供优质的乳、肉、皮、绒等多种畜产品而获得“高原之舟”之美誉(Miaoet al., 2015)。麦洼牦牛主要分布在四川省阿坝藏族羌族自治州红原县、若尔盖县、阿坝县的草原牧区，是我国动物遗传资源中珍稀的畜种资源和基因库(李世林等，2014)。已有研究探索了牦牛生长发育过程中的肺脏(何俊峰等，2009)、肾脏(王婷等，2017) 和皮肤(于川等，2017)等形态学和组织学变化，而在分子水平上的研究仅限于肺脏(Xinet al.,2019)和肌肉(Maet al.,2019; Xinet al., 2019) 的生长发育机制，未见牦牛肝脏的生长发育相关研究。因此，本研究对不同发育阶段的麦洼牦牛肝脏进行转录组测序，旨在深入理解牦牛肝脏组织生长发育过程的变化，为高原牦牛肝脏的生长发育机制研究提供数据和借鉴。

1 研究方法

1.1 实验动物与样品采集

实验动物分为1 日龄组(LD)、15 月龄组(LM)和5 岁龄组(LY)的健康麦洼牦牛，每组3 个生物学重复，体重分别为(12.35 ± 1.85) kg、(98.88 ± 2.50) kg 和(268.55 ± 27.82) kg，2019 年7 月采自四川省阿坝藏族羌族自治州红原县。空腹24 h 后快速放血处死，取LD 组、LM 组和LY 组麦洼牦牛肝脏组织放入液氮中冷冻，并于-80℃保存备用。

1.2 RNA提取、文库构建及高通量测序

用RNAiso Plus(TaKaRa，日本)提取总肝脏组织RNA，用1%琼脂糖凝胶电泳及Nano-Photome‐ter®分光光度计(Implen，美国)分别检测RNA 的完整性和提取质量，其OD260/280值为1.8~2.2，则质检合格，可进行下一步分析。采用去除核糖体RNA 的方法构建本次实验所需牦牛肝脏组织链特异性文库(Parkhomchuket al., 2009)，随后分别使用Qubit2.0 Fluorometer (Invitrogen，美国) 和Agi‐lent 2100 bioanalyzer (Agilent Technologies，美国)对文库进行初步定量及质量检测，最后采用qRTPCR (Bio-Rad，美国)对文库的有效浓度进行准确定量以保证文库质量。由北京诺禾致源有限公司使用二代高通量测序平台Illumina (HiSeqTM2500)(Illumina，美国)进行测序。

1.3 转录组测序数据处理及生物信息学分析

用HISAT2 v2.0.5 软件(https://ccb. jhu. edu/software/hisat2/manual. shtml#usage) 将质控后的有效测序数据(Clean Reads) 与牦牛参考基因组(BosGruv2.0：ftp://ftp. ensembl. org/pub/release-97/fasta/bos_mutus/) 进行快速精确的比对。获得匹配数据(Mapping Reads) 后比对结果进行质量评估，得到各样品的比对效率及在基因组上的定位信息(Mortazaviet al., 2008)。 分 别 使 用StringTie(1.3.3b) 软件和Subread 软件中的featureCounts 进行新转录本组装(Perteaet al., 2015) 和基因表达水平定量分析(Garberet al., 2011; Liaoet al.,2014)，并用FPKM (Fragments Per Kilobase of tran‐script sequence per Millions base pairs sequenced) 对测序深度和基因长度进行了校正(Brayet al.,2016)。基于负二项分布的DESeq2 (1.16.1) 软件对Read Count 进行组间表达差异基因统计分析(Loveet al., 2014)，标准为组间样品|log2 (Fold Change)| > 0 & padj < 0.05。采用层次聚类对基因的FPKM 值进行聚类分析，对行进行均一化处理。采用Cluster Profiler (3.4.4) 软件对差异基因集进行GO 功能富集分析与KEGG 通路富集分析，预测其可能参与的生物学过程和功能，并进行相应分类及统计。

1.4 qRT-PCR差异验证

采用qRT-PCR 验证转录组差异表达基因的表达，根据KEGG 富集分析选择显著富集(P=0.000 001 130) 的通路富集基因进行验证，以β-actin作为内参，引物详见表1。qRT-PCR 扩增体系为10 μL，模板cDNA 1 μL，上、下游引物各0.8 μL，ddH2O 2.2 μL，TB GreenTMPremix Ex TaqTMⅡ5.2 μL。扩增程序为：95℃3 min；95℃10 s，53.4℃(β-actin)55℃(验证基因)20 s，72℃20 s，共39 次循环；溶解曲线65℃~ 95℃每5 s增加0.5℃。将qRT-PCR 验证结果与转录组测序结果进行关联分析，检验RNA-Seq 数据可靠性。采用2-ΔΔCT法进行分析，通过SPSS 24.0 统计软件对数据进行显著性分析。

表1 qRT-PCR引物信息Table 1 Primers information of the genes for qRT-PCR

2 结果

2.1 牦牛肝脏RNA的提取与cDNA文库质量检测

取牦牛肝脏组织RNA，用琼脂糖凝胶电泳检测，可观察到两条清晰的条带，其OD260/280值为1.92 ~ 2.18，说明本次提取的RNA 完整性较高、质量较好。文库构建完成后，测得各肝脏组织样本的RNA 浓度均在1 032 ng/μL 以上。将文库稀释至1.5 ng/μL，测得文库有效浓度高于2.0 nmol/L，结果表明cDNA 文库质量合格，构建成功。

2.2 测序数据的分析与注释

通过双端测序法测得LD 组、LM 组和LY 组牦牛肝脏转录组平均每个测序样本约16.07 G、14.73 G 和16.19 G 原始碱基数(Raw Bases) [约10.72 亿条、9.82 亿条和10.80 亿条原始数据(Raw Reads)]，质控后得到约15.48 G、14.14 G 和15.57 G 有 效 测 序 碱 基 数 (Clean Bases) ( 约10.32 亿条、9.42 亿条和10.38 亿条有效测序数据)。测序所得数据碱基错误率均小于0.03%，且Q20 大于96.44%、Q30 大于91.00%，GC 含量在51.10% ~ 53.22%之间。以上结果表明该实验得到了高质量的转录组测序结果。将有效测序数据与参考基因组(BosGruv2.0) 进行序列比对，获取其在参考基因组上的位置信息。3 个实验组牦牛肝脏样本比对到参考基因组上的匹配数据分别占有效测序数据的87.13%、84.42%和81.03%，表明所选参考基因组作为参考进行组装注释可以满足分析需求。

2.3 新转录本预测与基因表达水平分析

使用StringTie (1.3.3b) 软件进行新转录本组装，组装后对新转录本进行Pfam、SUPERFAMI‐LY、GO、KEGG 等数据库注释，结果发现共有552 个新基因(表2)。RNA-Seq 相关性检查结果表明，每组样品间基因表达水平的相关系数R2均大于0.951，即各组样品之间表达模式的相似度较高，增加了实验的可靠性，符合实验要求。

表2 部分新转录本结果Table 2 Some results of new transcript

2.4 差异表达基因分析

将LD 组、LM 组和LY 组牦牛肝脏基因表达水平分析得到的数据采用DESeq2(1.16.1)进行分析，LD 组 与LM 组 筛 选 出1 144 个DEGs，其 中634 个DEGs 在LD 组中表达较高、510 个DEGs 在LM 组中表达较高；LD 组与LY 组筛选出1 228 个DEGs，其中751 个DEGs 在LD 组中表达较高、477 个DEGs 在LY 组中表达较高；LM 组与LY 组筛选出483 个DEGs，其中271 个DEGs 在LM 组中表达较高、212 个DEGs 在LY 组中表达较高(图1)。差异基因聚类分析结果表明表达模式相近的基因(纵坐标) 或样本(横坐标) 聚集效果明显，即同一基因(纵坐标) 在不同样本(横坐标) 中的表达情况一致(图2)。

图1 不同生长阶段牦牛肝脏差异表达基因统计. LD：1日龄组；LM：15月龄组；LY：5岁龄组Fig.1 Distribution of differentially expressed genes in yak at dif‐ferent growth stages. LD: 1-day-old; LM: 15-month-old; LY: 5-year-old

图2 不同生长阶段牦牛肝脏差异表达基因聚类热图. 红色代表高表达量，蓝色代表低表达量；LD：1日龄组；LM：15月龄组；LY：5岁龄组Fig. 2 Heat map of cluster of differentially expressed genes in yak at different growth stages. Red and blue represent high and low expres‐sion,respectively;LD:1-day-old;LM:15-month-old;LY:5-year-old

2.5 GO功能与KEGG通路富集分析

基于上述对生长过程中差异表达基因在3个生长阶段的表达模式分析，对3个生长阶段的高表达基因进行差异表达基因功能分析，结果显示，LD组、LM 组和LY 组中分别有325 个、85 个和84 个表达水平高于其他两组的显著差异表达基因(图3)。

图3 三个生长阶段两两比较分析获得的高表达编码基因数韦恩图. LD：1日龄组；LM：15月龄组；LY：5岁龄组Fig.3 Venn diagram of paired comparison of genes with higher expression among stages. LD:1-day-old;LM:15-month-old;LY:5-year-old

根据GO 功能分类富集，将注释基因分为生物过程(biological process, BP)、细胞组成(cellu‐lar component, CC) 和分子功能(molecular func‐tion,MF)三大功能类。LD 组、LM 组和LY 组高表达基因注释结果中分别富集到102 个、104 个和134 个显著差异注释条目(P< 0.05)，其中分别含有59个、51个、67个BP类，8个、12个、15个CC类和35 个、41 个、52 个MF 类，均选取最显著的30 个亚类绘制柱状图进行展示(图4A、C、E)，结果显示其占比最多的GO 条目是氧化还原相关过程、发育过程和代谢过程。根据KEGG 通路富集结果显示，LD 组、LM 组和LY 组高表达基因注释结果中分别富集到19 个、13 个和19 个显著差异注释信号通路(P<0.05)，均选取前10 个显著KEGG通路绘制散点图进行展示(图4B、D、F)，结果显示其占比最多的KEGG 通路是PI3K-Akt信号通路、粘着斑和ECM-受体相互作用。

图4 不同生长阶段高表达基因GO富集分析柱状图与KEGG富集分析散点图. A：LD组高表达基因GO富集；B：LD组高表达基因KEGG富集；C：LM 组高表达基因GO 富集；D：LM 组高表达基因KEGG 富集；E：LY 组高表达基因GO 富集；F：LY 组高表达基因KEGG 富集；GO富集分析中横坐标代表GO富集条目，纵坐标-log10(P value)代表GO富集条目显著性水平；KEGG富集分析中点的大小代表注释到KEGG信号通路上的基因数，颜色从红到紫代表富集的显著性大小，横坐标Gene Ratio代表注释到KEGG 信号通路上的差异基因数与差异基因总数的比值，纵坐标代表KEGG信号通路；LD：1日龄组；LM：15月龄组；LY：5岁龄组Fig.4 Histogram of GO enrichment analysis and scatter diagram of KEGG enrichment analysis of highly expressed genes in different growth stag‐es. A: Enrichment of GO in LD group; B: Enrichment of KEGG in LD group; C: Enrichment of GO in LM group; D: Enrichment of KEGG in LM group; E: Enrichment of GO in LY group; F: Enrichment of KEGG in LY group; In enrichment analysis of GO, X-axis represents GO terms and -log10 (P value) of Y-axis represents significance level of GO terms; In enrichment analysis of KEGG, the size of the dots represents the number of genes annotated on KEGG pathways,the color from red to purple represents the significance of the enrichment,and Gene Ratio of X-axis represents the ratio of the number of DEGs annotated to KEGG pathways to the total number of DEGs and Y-axis represents KEGG pathways; LD: 1-day-old;LM:15-month-old;LY:5-year-old

2.6 qRT-PCR差异验证结果

从高通量测序结果中挑选8 个DEGs(CYP7B1、PGFS2、CYP1A1、UGT2C1-1、UGT2C1L、HSD11B1、CYP2C19、UGT2C1-2) 进行qRT-PCR 验证。根据标准曲线结果可知，其扩增效率E 值为95% ~ 100%、相关系数R2达到0.997 以上，扩增效率和相关系数线性关系满足分析要求。计算RNA-Seq 的FPKM 值与qRT-PCR 结果之间的皮尔逊相关性，其皮尔逊相关系数r为0.869，满足分析要求。结果表明，qRT-PCR 验证结果与RNA-Seq 测序结果一致，表明基于转录组分析DEGs 的表达结果可靠(图5)。

图5 RNA-Seq 差异表达基因的qRT-PCR 验证. *代表LM 组与LD 组相比差异显著P<0.05，#代表LY 组与LD 组相比差异显著P<0.05；LD：1日龄组；LM：15月龄组；LY：5岁龄组Fig.5 Verification of RNA-Seq by qRT-PCR. *indicates that the difference is significant between LD and LM(P<0.05),#indicates that the dif‐ference is significant between LD and LY(P<0.05);LD:1-day-old;LM:15-month-old;LY:5-year-old

3 讨论

在已有的牦牛肝脏转录组的相关研究中，以牦牛肝脏为研究对象探讨镉胁迫下牦牛肝脏组织转录组特征(刘祥军，2017)，但还未见不同生长阶段牦牛肝脏转录组的相关报道。采用RNA-Seq 技术对1 日龄、21 日龄、2 岁龄的猪肝脏组织进行转录组测序，分别筛选出396 个、260 个和324 个表达水平高于其他两组的DEGs (贺伸，2017)。对20 周和30 周的鸡肝脏组织进行了转录组测序，分别筛选出1 485 个和1 082 个高表达DEGs (Liet al., 2015)。本实验采用RNA-Seq 技术对1 日龄、15月龄和5岁龄的健康麦洼牦牛肝脏进行转录组测序，分别筛选出325 个、85 个和84 个表达水平高于其他两组的DEGs。与先前研究结果一致，1 日龄肝脏相比15 月龄和5 岁龄肝脏筛选出的DEGs最多，其原因可能是1日龄与15月龄和5岁龄牦牛之间发育程度相差较大导致，说明1日龄肝脏发育与15 月龄和5 岁龄相比更为关键。该研究结果表明，不同生长阶段牦牛肝脏差异表达基因的变化与生长发育的分子调控机制密切相关，为进一步探索牦牛肝脏生长发育相关新基因提供参考。

通过转录组学的方法筛选出了3个生长阶段牦牛肝脏的差异表达基因，功能分析的结果显示3个年龄阶段髙表达的基因富集至各个生长阶段非常重要的生物学过程中。对筛选出的前30项GO条目分析发现，占比最多的GO 条目是氧化还原相关过程、发育过程与代谢过程。氧化还原酶种类丰富，作为酶中数目最多且最重要的一个分支，广泛应用于生物催化级联反应中，并对营养物质的转化以及矿质营养元素的循环具有重要意义(伍开琳和陈永正，2018)。同时，机体参与氧化还原反应伴随着大量电子和铁离子的传递和流动(Junget al.,2016)。肝脏发育过程中，初生时功能尚不健全，随着生长功能逐渐完善，至体成熟阶段时功能健全，代谢来自食物的各种成分，包括有毒有害物质也主要是在肝脏中进行(仲召鑫等，2019；石李梅等，2020)。1 日龄与15 月龄的高表达基因显著富集至各种氧化还原酶活性、钙离子结合等氧化还原相关过程和单、多细胞等发育过程中(P<0.05)，表明这两个阶段的肝脏组织均处于快速生长发育的年龄阶段。其中，1 日龄时，部分基因富集至胚胎发育调控过程，表明此时的肝脏组织还处于胚胎发育调控的延伸阶段；15 月龄时，部分基因富集至T细胞受体信号通路等免疫过程中，表明此时的肝脏组织已经承担部分机体的免疫功能。代谢过程是生命所需的最基本活动之一，并参与细胞或动物体内必不可少的生理活动(Lempradlet al., 2015)。肝脏是机体物质代谢的中枢器官，是许多物质合成、分解和转化等代谢过程发生的主要场所(王绍庆，2012)。5 岁龄时的高表达基因主要显著富集至肝脏营养物质代谢过程中(P<0.05)，包括氨基酸代谢、小分子代谢等。表明各个生长阶段高表达基因显著富集过程与肝脏生长发育特性和功能相吻合，由此推测氧化还原酶活性、钙离子结合等氧化还原相关过程，单、多细胞等发育过程与氨基酸代谢、小分子代谢等代谢过程在不同发育阶段牦牛肝脏的生长发育中具有重要的调控作用。

对筛选出的生长发育相关的KEGG 通路分析发现PI3K-Akt信号通路、粘着斑和ECM-受体相互作用占比最多。PI3K-Akt 信号通路可调节DNA 甲基化在调控细胞的增殖分化、生长发育、适应环境和疾病发生发展上发挥重要作用(梅传忠，2010；邓大君，2014)。粘着斑激酶(focal adhesion kinase,FAK)主要参与整合素激活的FAK 信号转导通路，通过ECM 途径浸润邻近组织、发生粘着斑的动态改变，对肿瘤细胞发生发展进程都发挥着重要作用(Liuet al.,2018;Weiet al.,2020)。ECM与其相应受体结合后能够向细胞发出信号从而启动或调节基因的表达，该信号传导在组织器官形态学发生、胚胎发育、细胞增殖分化等过程中起重要作用(王思源等，2019；Zhanget al., 2020)。而机体发育、组织修补和血管生成等多种生理状态均需要进行胶原代谢降解胶原(包新见，2006)，这与本实验结果在PI3K-Akt 信号通路、粘着斑与ECM-受体相互作用中主要富集了胶原蛋白α 链(COLA)家族基因结果一致。其次，整合素α 亚基(ITGA) 也富集于PI3K-Akt 信号通路、粘着斑与ECM-受体相互作用中。而整合素可与胞外配体结合，影响黏着斑与细胞骨架结合形成ECM-整合素-细胞骨架蛋白的复合物，介导PI3K-Akt 等多种细胞内信号通路(杨洪娟，2017)，这与本实验结果在PI3K-Akt信号通路、粘着斑与ECM-受体相互作用中富集了整合素α 亚基(ITGA)家族基因结果一致。由此可知，PI3K-Akt 信号通路、粘着斑和ECM-受体相互作用在各生长阶段牦牛肝脏大量富集主要是与细胞受体介导的细胞内信号传导活动相关，这也更加说明了在不同生长阶段牦牛肝脏生理活动在分子层面上表现出的巨大差异性和特殊性。

本研究利用RNA-Seq 技术探讨不同生长阶段牦牛肝脏的生长发育机制，牦牛肝脏生长发育过程可能是通过调节氧化还原过程、发育过程、代谢过程，以及调节PI3K-Akt 信号通路、粘着斑和ECM-受体相互作用信号通路的基因进行，为深入研究高原动物肝脏生长发育机制奠定理论基础。