崂山奶山羊不同泌乳期乳腺组织可变剪接的鉴定与分析

关键词： 奶山羊； 乳腺； 可变剪接； RNA测序

中图法分类号： S827 文献标识码： A 文章编号： 1000-2324（2024）01-0108-08

可变剪接（Alternative splicing，AS）是指转录后形成的前体RNA经过剪接体和剪接因子的相互作用，最终形成多种成熟的mRNA 的过程[1]，是一种广泛提高转录组和蛋白质组复杂性的分子机制，控制发育程序以及高等真核生物对环境的反应[2-5]。随着RNA测序技术和分析工具的最新发展，RNA-seq 数据提供了丰富的研究资源，可用于揭示不同生物体在许多生物过程中的新颖可变剪接事件和可变剪接调控机制[6-8]。目前，山羊可变剪接机制的研究报道较少，有关可变剪接事件在奶山羊乳腺发育与泌乳生理活性中的机理尚不清楚。因此，研究可变剪接调控奶山羊乳腺生理的机制，可为更好地理解乳腺发育和泌乳调控提供必要的理论基础[9]。

当前，在畜禽体内可变剪接机制的功能探索研究表明，可变剪接在动物体的生长发育过程中起到重要的调控作用。例如：对4个不同发育时期猪睾丸组织进行可变剪接鉴定和分析，发现第1 个外显子可变剪接和最后1 个外显子可变剪接事件类型与睾丸素合成和分泌密切相关[10]。通过对不同发育时期北京鸭胸肌和皮脂RNA-seq数据进行可变剪接分析，发现胸肌特有可变剪接基因显著富集到骨骼肌发育相关的GO条目中，而皮脂特有的可变剪接基因显著富集到油脂合成的GO条目中[11]。研究人员对山羊五种组织类型（心脏，肾脏，腿部肌肉，肝脏和脾脏）不同发育阶段差异可变剪接基因进行分析，结果显示鉴定得到的差异可变剪接基因主要参与器官的功能和发育[12]。

本实验室以崂山奶山羊乳腺组织为研究模型，基于RNA-Seq 数据分析了可变剪接事件对不同泌乳时期乳腺生理的功能，揭示了其在乳腺生理学和发育过程中基因表达的复杂调控网络[13-18]，探索可变剪接调控奶山羊泌乳生理的分子机制，为奶山羊的乳腺发育和泌乳机理研究提供科学理论基础，同时也为奶山羊育种改良工作提供理论支持。

1 材料与方法

1.1 试验样本及转录组数据

本研究RNA-Seq 数据为本实验室先前课题研究数据，已上传NCBI SRA 数据库（SRA 号：SRA218662）。试验样本均来自于青岛奥特种羊场，四岁龄、第三胎次的5 只崂山奶山羊，分别在泌乳初期（产后20 d），泌乳高峰期（产后90 d），和泌乳末期（产后210 d），采集不同泌乳时期乳腺组织。试验羊只自由获取食物，无任何疾病，并在相同条件下饲养和管理。经全身麻醉后通过外科手术分别采集三个泌乳时期乳腺组织，并迅速放入无RNase 冷冻管，然后立即浸入液氮罐中保存。RNA提取、混池测序质量检测和序列比对质量较好[14]。

1.2 可变剪接特征分析

本研究中使用ASprofile 软件根据转录本的剪接位点对三个文库分别进行可变剪接事件和基因的鉴定[11]。分析结果存在7 种可变剪接事件类型（图1）：可变5′或3′端剪接、单内含子滞留、多内含子滞留、单外显子跳跃、多外显子跳跃、第1 个外显子可变剪接、最后1 个外显子可变剪接。

1.3 差异可变剪接基因检测及注释分析

使用rMATS 软件检测三个文库之间差异可变剪接基因[20]。FDR≤0.05 为差异显著性判断标准，小于该阈值的差异可变剪接基因定义为显著性差异基因。

通过GO 数据库和KEGG 数据库进行生物学功能分析[21-22]，筛选满足P≤0.05 的GO条目与KEGG 通路进行下游分析。

1.4 可变剪接调节蛋白和目的基因调控网络的构建

基于GO和KEGG注释结果，筛选与乳腺发育及泌乳相关的可变剪接基因，在STRING数据库[23]中进行蛋白互作分析。然后使用Cytoscape软件[24]构建基因间的调控网络。

2 结果与分析

2.1 测序数据质量评价与分析

对测序数据经统计与质控分析后，泌乳初期、高峰期和末期文库中获得原始序列分别为36 336 892 条、35759380条和36469596条。与山羊参考基因组比对结果表明，泌乳初期比对到基因组上的序列数为34864794条，高峰期33888773条，末期34963722条，比对率分别为96.79%、95.60%和96.64%（表1）。

2.2 可变剪接的鉴定

泌乳初期、高峰期和末期三个文库共鉴定到151 503 个可变剪接事件，对应17 627 个可变剪接基因，占山羊总注释基因的64.92%。其中，泌乳高峰期鉴定到基因数目最多，共计16 943个，泌乳初期和末期分别为16 085、11 509 个（图2-A）；可变剪接频率分析表明：泌乳高峰期鉴定到可变剪接频率为3.24，低于泌乳末期的4.04，泌乳初期最低，为3.11，三个泌乳期平均可变剪接频率为3.46。可变剪接类型分析表明：三个泌乳期均以第1 个外显子可变剪接、最后1 个外显子可变剪接、多外显子跳跃和可变5′或3′端剪接为主要可变剪接事件，四种事件类型分别占到泌乳期可变剪接事件总数目的94.16%、93.79% 和94.52%；多内含子滞留事件数目在三个泌乳期中最低，为稀有事件类型，分别占比0.04%、0.02%和0.03%（图2-B）。三个泌乳期可变剪接基因在1-29 号染色体中均在18 号染色体分布最多，在27号染色体分布最少（图2-C）。

2.3 差异可变剪接基因分析

三个泌乳期共有可变剪接基因14783个，占剪接基因总数目83.87%。泌乳初期与高峰期共有可变剪接基因数目最多，泌乳初期与末期数目最少，分别为816 和168。特有可变剪接基因数目在泌乳高峰期最多、泌乳初期次之、泌乳末期最少，分别为983、317 和198（图3-A）。初期与高峰期、初期与末期和高峰期与末期三组分别鉴定到5 848、4870 和5365 个差异可变剪接基因，不同泌乳期间共鉴定到差异可变剪接基因6323个，其中共有差异可变剪接基因4 018个，占总数目63.55%（图3-B）。在初期与高峰期间检测到1173个显著性差异可变剪接基因，初期与末期间检测到884个，高峰期与末期检测到1 107 个，外显子跳跃类型基因数目在三组中均占比最高，其次为3′端可变剪接位点和5′端可变剪接位点，三种类型基因数目分别占每组显著性差异可变剪接基因总数的91.56%、86.06%和90.87%（图3-C）。

2.4 可变剪接基因GO注释分析

为了进一步探究差异可变剪接基因的功能，对初期与高峰期、初期与末期和高峰期与末期三组5 848、4870和5365个差异可变剪接基因进行分析。在细胞组分中，主要有438 个基因与胞质相关，397 个基因与细胞核相关，263个基因与核质相关，92个基因与高尔基体相关。在生物学过程中，主要有130 个基因与以DNA为模板的转录调控相关，84个基因与RNA聚合酶II 启动子转录的负调控相关，78 个基因与GTPase 活性的正调控相关，61个基因与蛋白质磷酸化相关；在分子功能中，主要有188个基因与金属离子结合相关，183 个基因与蛋白质结合相关，176 个基因与ATP 结合相关，131个基因与poly（A）RNA结合相关（图4）。

此外根据注释结果，在差异可变剪接基因注释结果中筛选到了一些可变剪接调节蛋白和乳腺发育及泌乳相关的可变剪接调节蛋白48 个，筛选到与乳腺发育及泌乳相关的基因共计65个，其中23个基因注释到乳腺上皮细胞增殖，15个基因注释到乳腺导管形态发生中分支的调控，9个基因注释到乳腺腺泡发育，5个基因注释到乳腺脂肪发育，10个基因注释到乳腺导管末梢乳芽生长（表2）。

2.5 可变剪接基因KEGG通路分析

KEGG通路分析结果表明，泌乳初期与高峰期、泌乳初期与末期、泌乳高峰期与末期差异可变剪接基因富集结果具有非常高的重合性，主要涉及泛素介导的蛋白水解、胞吞作用、甘油磷脂代谢、催产素信号通路、钙信号通路、类固醇生物合成、氨基糖和核苷酸糖代谢、内质网中的蛋白质加工、半乳糖代谢等乳腺组织生理活性密切相关的通路（图5）。

2.5 可变剪接基因调控网络

根据差异可变剪接基因的GO 和KEGG 的注释结果，我们筛选出65 个和乳腺发育及泌乳相关的目的基因，同时根据GO注释结果筛选的48 个可变剪接调节蛋白，通过STRING 分析，得到46 个节点，64 对调控关系，利用Cytoscape 构建可变剪接调节蛋白和目的基因调控网络；基因PTEN 有14 个关系对，为本网络核心基因（图6）。

3讨论

众多研究证明，可变剪接在所有高等真核生物基因表达过程中广泛存在，鉴定分析可变剪接已成为真核基因调控领域的研究重点[25-26]。研究人员发现人类大约有95% 的基因发生可变剪接[27]，鸭约有46.12%的基因发生可变剪接[11]，牛约有21%的基因发生可变剪接[28]；本项研究共鉴定到151 503 个可变剪接事件，对应17 627 个可变剪接的基因，可变剪接基因数目占山羊总注释基因的64.92%，通过对比分析可以发现在不同物种之间，可变剪接基因占总注释基因数目存在巨大差异，推测可能和可变剪接事件的时序性表达、物种特异性有关[29]。可变剪接不仅在真核生物体基因表达过程中广泛存在，同时在生物体生长发育过程中存在特异性表达等特点[30]。在本研究中，泌乳高峰期鉴定到剪接事件和基因数目均明显高于泌乳初期和泌乳末期，结合已有研究报道乳腺组织在不同泌乳时期生理结构和功能会发生显著性变化[31-32]，我们推测可变剪接事件的变化可能和乳腺组织不同泌乳期的生理活性密切相关。

可变剪接基因的生物学功能分析在可变剪接研究中具有重要作用。本研究中，我们筛选到和乳腺细胞增殖、分化与凋亡等密切相关的GO条目，以及一些影响乳腺上皮细胞内物质转运和组织代谢等密切相关的GO条目，不同泌乳期特有剪接基因的功能具有明显特征，这与研究人员之前报道一致[11]，这些研究表明特有剪接基因功能具有高度的适应性。差异可变剪接基因GO注释和KEGG通路分析中，筛选出一些如：催产素信号通路、钙信号通路、糖代谢和磷脂代谢等与乳腺上皮细胞泌乳能力密切相关的通路。通过构建调控网络，筛选PTEN为参与乳腺发育与泌乳调控网络的核心基因。前期研究表明，干扰PTEN 基因，可影响奶山羊乳腺上皮细胞脂质合成相关基因的转录及脂肪酸组成的影响，显著影响乳脂合成[33]。在奶牛乳腺上皮细胞中，PTEN基因可靶向调节PI3K-AKT信号转导通路，影响乳腺上皮细胞β-酪蛋白，甘油三酸酯和乳糖的分泌，从而参与奶牛乳腺组织的发育和泌乳[34-35]。

4结论

在奶山羊乳腺不同泌乳期，可变剪接的发生具有阶段特异性，可变剪接基因与乳腺组织生理活性特点密切相关，参与了奶山羊乳腺组织发育及泌乳功能的调控作用。PTEN可参通过构建的核心调控网络参与乳腺发育与泌乳调控。