不同泌乳期牦牛乳腺组织的转录组学分析

董宝霞，杨玉莹，俞红贤，张勤文，高小龙，李莉，魏青，荆海霞*

（1.青海大学农牧学院，青海西宁 810016；2.青海大学生态环境工程学院，青海西宁 810016）

牦牛是一种特殊牛种，主要分布于平均海拔3 000 m左右的青藏高原及其邻近地区。相对荷斯坦牛、山羊、黑犏牛、马和骆驼等产乳动物而言，牦牛乳含有比其他种类的乳更高的干物质、蛋白质、脂肪、乳糖和矿物质等，且维生素含量很高，营养成分均衡，具有易消化、易吸收的特点，能满足婴幼儿日常营养需求，是一种理想的婴幼儿奶粉的原料。虽然牦牛乳具有上述诸多优点，但其产量低的缺点导致牦牛乳价格较高，限制了牦牛乳的大规模推广，如何提高牦牛的产奶量以及牦牛乳的品质是当前亟需解决的问题。

乳汁是由乳腺分泌的，是幼儿和幼畜的主要食物，其分泌主要是在催乳素、雌激素、肾上腺糖皮质激素和胰岛素等激素的综合作用下进行的，而乳腺是雌性动物抚育幼崽过程中必不可少的一部分。目前，为了提高乳畜的产奶量，相关学者利用多种技术对乳腺组织泌乳调控机制进行了大量研究。近年来，随着高通量测序技术的发展，转录组学的研究取得了很大进展，并在各个领域有了广泛应用。RNA 高通量测序技术（RNA Sequencing,RNA-Seq），是指将RNA进行逆转录以及PCR 扩增形成cDNA，然后对其进行测序，从而可以精确地定量几乎所有的转录表达。转录组技术目前已广泛应用于奶牛乳腺组织的研究中。Dai 等、Gao、Lin 等曾利用转录组测序技术对奶牛乳腺组织进行了转录组分析；Fan 等对3 头干奶期牦牛和2 头泌乳期牦牛的乳腺组织样本进行了转录组学分析，发现了一些牦牛在不同生理期的差异表达基因（Differentially Expressed Genes,DEGs）以及代谢通路等。虽然转录组技术在牦牛身上也有相关研究，但有关牦牛乳腺组织泌乳机制方面的研究相对较少。因此，为系统了解牦牛乳腺组织泌乳调控的机制，本研究利用RNA-Seq 技术，将处于不同泌乳期的牦牛乳腺组织进行转录组学测序，并进行对比分析，以期了解不同泌乳阶段的牦牛乳腺组织的DEGs 及所涉及的泌乳相关通路，为进一步研究牦牛乳腺组织中的泌乳分子机制和提高牦牛泌乳量及乳品质提供一定的理论基础。

1 材料与方法

1.1 实验乳腺样品采集 2020 年7 月于青海省祁连县某养殖场选取处于泌乳早期（Early Lactation Period,ELP）、泌乳中期（Mid-Lactation Period,MLP）、干奶期（Dry Lactation Period,DLP）的健康牦牛各3 头，其中ELP牦牛处于产后10 d 左右，MLP 牦牛处于产后70 d 左右，DLP 牦牛处于产后300 d 左右。将每头牦牛的2个后乳区表皮进行剥离，分别采集乳腺组织后迅速切成1 cm×1 cm×1 cm 的小块，放至无RNA 酶的离心管中并立即投入液氮中保存备用。

1.2 转录组文库构建及测序 将收集到的牦牛乳腺组织样品用干冰送至北京奥维森基因科技有限公司进行总RNA 提取并使用带有Oligo(dT) 的磁珠富集mRNA 以进行cDNA 文库的构建，文库构建流程如图1。然后，通过Illumina-Hiseq（HiSeqTM2500/4000）平台进行转录组测序。

图1 转录组文库构建示意图

1.3 质量控制 将Illumina 测序得到的raw reads（FASTQ Format）进行数据过滤后得到clean reads，并计算Q20（Phred 数值大于20 的碱基占总体碱基的百分比）、Q30（Phred 数值大于30 的碱基占总体碱基的百分比）及clean reads 的GC 含量。

1.4 参考基因组序列比对 比对软件选用TopHat2，将转录组测序reads 比对到基因组上，然后通过分析比对结果识别外显子之间的剪接点（Splicing Junction）。TopHat2 分析流程如图2 所示。

图2 TopHat2 分析流程图

1.5 定量基因表达水平 参考文献，采用HTSeq软件分析每个样品的基因表达水平，模型为union。计算各基因的RPKM 值。

1.6 转录组数据分析

1.6.1 差异表达分析 用DESeq R 软件包（1.10.1）分析转录组在条件不同时的差异表达，计算log2（Fold Change）及其统计检验的值。本研究以|log2（FoldChange）|>2 且qvalue<0.005 作为筛选标准来确定基因表达的显著性。

1.6.2 差异表达基因的GO 和KEGG 富集分析 GO 富集分析用GOseq 软件。同时，应用超几何检验对通路进行显著性富集分析，以KEGG pathway 为单位找出DEGs 显著富集的通路。

2 结果与分析

2.1 转录组测序数据分析

2.1.1 碱基质量分布检查 检测原始序列错误率分布情况，用Phred 数值（Phred Score，Qphred）来体现碱基质量。结果表明，不同泌乳期的牦牛乳腺组织样品的Q20 均在97.12% 以上，说明在测序时正确识别率大于99%的碱基总数占所测定的碱基总数的97.12%；不同泌乳期各样品的Q30 均在92.82% 以上，说明在测序时正确识别率大于99.9% 的碱基数占到所检测碱基总数的92.82%。测序结果显示，每个乳腺组织样品的clean bases 总数均大于9.66 G，GC 含量在46.96%～50.42%，表明测序质量较高。

2.1.2 测序数据比对分析 将每个乳腺组织样品的clean reads 比对到参考基因组上，能定位到基因组上的测序序列的比例大于99.98%。每个乳腺组织样品测序的clean reads 与参考基因组的比对率均达到95.77%以上。说明数据使用率可靠，所选择的参考基因组适合于后续分析。

2.2 基因差异表达分析及筛选 差异表达分析的目的是找出不同样本间存在的DEGs，结果见图3，每个时间段内所得到的DEGs 也是不同的，将DLP 与ELP、MLP 进行 对比，总DEGs 分别为3 137、394个；将MLP 与ELP 进行对比，总DEGs 为281个。不同泌乳阶段乳腺组织前10个DEGs 如表1、2、3 所示。

表1 DLP vs MLP 组乳腺组织前10个显著DEGs

图3 差异表达基因的火山图

从不同组别的DEGs 中可以发现，很多基因在ELP和MLP 表达量高，但是在DLP 表达量明显降低甚至不表达，如跨膜通道样蛋白5（）、乳过氧化物酶（）、糖蛋白2（）等；反之，在ELP 和MLP 表达量低或者不表达的基因在DLP 的表达量升高，如细胞色素P450 家族2 亚家族F 成员1（）、突触结合 蛋白4（）、短链脱氢酶/还原酶家族9C 成员7（）等。

2.3 GO 聚类分析和 KEGG 通路分析

2.3.1 各时间段相关上调、下调基因的GO 富集分析结果 用Gene Ontology 通过GO 富集功能分析各时间段上、下调相关基因，FDR<0.05 有意义，将DLP vs MLP、DLP vs ELP、MLP vs ELP 对照组各泌乳阶段的差异表达上、下调基因做GO 功能分析，结果显示，3组牦牛乳腺组织的DEGs 包括分子功能、生物学过程、细胞组分3 大类。在DLP vs MLP 对照组中生物学过程类别内富集的前5 条GO 类别如表4 所示。各组相关上调、下调DEGs 分别被注释到GO term，如图4 所示。结果显示，3 组牦牛乳腺组织中，最为富集的是多细胞生物过程（Multicell Organismal Process）；其次是生物过程的正调控（Positive Regulation of Biological Process）和发育过程（Development Process）。其中与细胞相关的GO 类别占大多数。分子功能类别中，富集的前10 条GO 类别全都与酶活性相关，其中最为富集的是氧化还原酶活性（Oxidoreductase Activity）、单加氧化酶活性（Monooxygenase Activity）。在细胞组分类别中，最为富集的类别是胞外区（Extracellular Region）、细胞外围（Cell Periphery）、质膜（Plasma Membrane）。

图4 GO 富集柱状图

表4 生物学过程富集前5 条GO 类别

2.3.2 各泌乳阶段相关上调、下调基因的KEGG 富集分析结果 为了进一步了解这些DEGs 所参与的生物化学代谢途径和信号传导途径，本研究利用KEGG 公共数据库进行分析（图5）。可知，在DLP vs ELP、DLP vs MLP、MLP vs ELP 中的DEGs 分别涉及29、15、22条通路，DLP vs ELP 组富集程度最高的前10 条通路中，最为富集的通路是味觉传导（Taste Transduction），其次是ECM 受体相互作用（ECM-receptor Interaction）和轴突导向（Axon Guidance）。DLP vs MLP 组富集程度最高的前10 条通路中，最为富集的通路是味觉传导（Taste Transduction），其次是药物代谢细胞色素P450（Drug Metabolism -cytochrome P450）和甘氨酸、丝氨酸和苏氨酸代谢（Glycine,Serine and Threonine Metabolism）。MLP vs ELP 组富集程度最高的前10 条通路中，最为富集的通路是卵巢类固醇生成（Ovarian Steroidogenesis），其次是甲状腺激素合成（Thyroid Hormone Synthesis）和味觉传导（Taste Transduction）。

图5 差异基因KEGG 富集散点图

表2 DLP vs ELP 组乳腺组织前10个显著DEGs

表3 MLP vs ELP 组乳腺组织前10个显著DEGs

3 讨 论

本研究采用了RNA-Seq 技术对处于ELP、MLP 和DLP 的牦牛乳腺组织样品进行了转录组学分析，结果共检测到3 812个DEGs，其中2 907个表达上调，905个表达下调。与ELP 相比，DLP 有2 643个基因表达上调，494个基因表达下调；与MLP 相比，DLP 表达上调的基因有142个，表达下调的基因有252个；而与ELP 相比，MLP 有122个基因表达上调，159个基因表达下调。GO 分析显示富集最丰富的条目与物质代谢、化学反应等相关，KEGG 分析表明DEGs 主要富集在味觉传导、ECM 受体相互作用、卵巢类固醇生成。Fan等也曾研究了牦牛干奶期和哺乳期乳腺组织的转录组，共发现360个DEGs，其GO 分析表明最丰富的条目与蛋白质结合、免疫系统等功能相关，KEGG 分析表明DEGs 主要集中于类固醇生物合成、HIPPO 信号转导、胰岛素信号转导。因此，本研究也同该研究一样发现了牦牛在不同泌乳期的DEGs，并对其进行了GO 和KEGG 分析，为今后提高牦牛的乳产量以及提升乳品质有着一定的指导意义。

3.1 物质代谢和转运与乳糖、乳蛋白和乳脂合成 将本实验所得到的不同泌乳期的DEGs 进行GO 富集分析和KEGG 通路富集分析发现，泌乳启动过程中的许多高表达的基因与物质代谢和转运活动有着密不可分的关系，与前人在牛和小鼠乳腺中的研究结果一致。乳汁主要是通过新物质合成以及从血液中吸收2 部分进行，并且这个过程中需要ATP 和相应酶的参与。血液中含有的维生素、无机盐、球蛋白以及激素等可以被乳腺获取后直接转变为乳汁成分并最终合成乳糖、乳脂和乳蛋白等物质。伴随着泌乳过程的启动，乳腺对氨基酸、葡萄糖及其他营养物质的需求量大大增加。

乳糖是哺乳动物乳汁中特有的、由葡萄糖和半乳糖结合形成的一种化合物，它的形成对乳汁的生成来说有着至关重要的作用。有关于山羊的研究表明，在分娩后1 周的乳腺比分娩前9 天所吸收的葡萄糖的量有非常显著的增加。本研究所得到的数据显示在泌乳早期-乳白蛋白的表达上调。血液中的游离氨基酸可被乳腺用来合成-乳球蛋白、酪蛋白和-乳清蛋白。本研究数据显示，与DLP 相比，ELP 乳蛋白基因表达量增高显著。牛的前胃内有大量的微生物，其发酵之后会生成很多产物，包括乙酸、丁酸等，最后合成C4～C16 脂肪酸，而且血液中的乳糜微粒和脂蛋白分解后也会产生脂肪酸，这些脂肪酸被乳腺转化生成乳脂，其几乎都是甘油三酯。脂蛋白脂肪酶（）可以将甘油三酯进行水解，所生成的产物是脂肪酸，这是机体能量的一个重要来源。泌乳初期脂蛋白脂肪酶的表达显著增加。其他与脂类合成相关的基因也在泌乳初期发生了显著变化，包括等。

3.2 泌乳性状相关候选基因的分析 二酰甘油酰基转移酶2（）是一种生物体内极其重要的酶。作为甘油酰基转移酶的作用是催化二酰甘油与脂肪酸酰基反应产生三酰甘油。基因在动物脂肪代谢旺盛的组织中表达水平很高，说明其与脂肪代谢及脂类沉积关系密切。鄢胜飞等在摩拉水牛基因上发现了9个SNPs 位点，并表明这些位点可能对摩拉水牛的生长性状和乳脂性状有着重要影响。郭将等通过PCRRFLP 技术检测草原红牛基因的遗传多态性，发现AA 型个体乳中所含有的脂肪和干物质显著高于AB型和BB 型。赵国丽等对中国荷斯坦奶牛基因第4 外显子进行了研究，发现81 bp（G/A）、474 bp（C/G）、621 bp（C/T）位点对305 d 的产奶量均有显著影响。目前现有报道中，多数是作为普通牛泌乳性状的基因，因此可以在牦牛上做进一步研究。

脂肪酸合成酶（）是一种酶，对脂肪酸合成有关键意义，在脂肪酸的形成中发挥重要作用。王小龙研究发现FASN 第34 外显子的2 种基因型对乳脂率有着显著的影响。贺一博等应用RFLP-PCR 技术对中国荷斯坦奶牛进行多态性检测发现基因存在一个突变位点g.51403203C>T（7 833 bp），而且该处存在不同基因型，参与调控奶牛乳脂的合成和代谢。马彦男等利用PCR-SSCP 技术对中国荷斯坦奶牛基因外显子34 进行了分析，结果发现2个等位基因，其中，AA 型个体305 d 时的乳脂含量显著高于AB 型个体，所以该基因在提高乳脂率方面具有重要意义。本研究结果中发现，基因在泌乳过程中表达，尤其是在ELP，其表达量对比DLP 上调显著。

诱导细胞调亡FF45 样效应因子C（），于1992 年在小鼠体内首次发现。激素敏感性脂肪酶（）是一种脂肪分解的限速酶，而直接作用于脂滴，影响HSL 在脂滴表面的所在位点，使脂质的分解速率显著降低甚至抑制其分解，最终使得脂质大量积累。有研究表明，如果在泌乳期分泌的乳汁品质高，其乳腺组织中的表达量较其他阶段是最高的，反之亦然。而且跟干奶期乳腺组织相比，泌乳初期乳腺组织中的表达量升高极显著。在本研究中经差异表达分析发现，DLP的表达量较ELP 下调了6.6 倍，所以基因的差异表达可能与牦牛乳汁中乳脂的合成有关，其可能在牦牛乳腺上皮细胞合成乳脂的过程中有着不容忽视的作用。这些结果也进一步说明本研究筛选出的差异表达基因与牦牛泌乳相关，也为筛选新的与牦牛泌乳相关的基因提供了参考。

3.3 PI3K-AKT 信号通路在泌乳过程中的作用 PI3KAKT 的下游靶点是mTOR（即哺乳动物雷帕霉素靶蛋白），即PI3K-AKT 是mTOR 的上游信号因子，所以mTOR 在调节乳蛋白的合成，特别是在翻译过程中起着关键作用。多项研究表明，AKT-mTOR 信号通路是乳蛋白的重要信号通路，从基因表达以及蛋白质合成2 方面对乳汁的合成与分泌进行调控。赵萌曾报道PI3K-AKT 信号通路与JAK-STAT 信号通路、MAPK 信号通路相互作用，这些信号途径在调控乳成分的合成和分泌等生物过程发挥重要作用。本研究结果表明有700多个与泌乳曲线变化一致的基因显著富集在PI3K-AKT信号通路，其中，87个泌乳上升阶段的差异基因和7个泌乳下降过程的DEGs 富集在此信号通路中，因此，推测PI3K-AKT 信号通路可能在牦牛泌乳过程中有着不可忽视的作用。

4 结 论

本研究应用RNA-Seq 技术对不同泌乳期牦牛乳腺组织进行了转录组学分析，共检测出3 812个DEGs，其中2 907个表达上调，905个表达下调，并对这些DEGs进行了相关功能分析。其中，GO 分析显示富集最丰富的条目与物质代谢、化学反应等相关；KEGG 分析表明DEGs 主要富集在味觉传导、ECM 受体相互作用、卵巢类固醇生成。本研究结果为进一步了解不同泌乳期牦牛乳腺组织的泌乳机制及牦牛泌乳的特殊性提供分子基础。