雌性阿什丹牦牛背最长肌三个出生后年龄段的比较全长转录组学与非靶向代谢组学分析

《Biomolecules》：Comparative Full-Length Transcriptomic and Untargeted Metabolomic Profiling of Longissimus Dorsi Muscle in Female Ashdan Yaks at Three Postnatal Ages

　　

为了表征雌性阿什丹牦牛背最长肌中与年龄相关的分子差异，本研究结合牛津纳米孔技术（ONT）全长转录组测序和非靶向液相色谱-质谱联用（LC-MS）代谢组学，对27只独立动物进行了分析，这些动物分为三组：10日龄（M0）、6月龄（M6）和54月龄（M54），每组9只。在BH q < 0.05且|log2FC| > 1的条件下，三组成对比较分别鉴定出1999、3480和2350个差异表达基因，其中M0与M54之间的差异最大。由这些基因集重建的网络识别出高连通性候选基因，包括涉及M0的比较中的ENSBGRG00000003669，以及M6与M54比较中的MRTO4和BYSL。在1842个保留的有注释的LC-MS峰特征中，分别有1487、1474和1285个在BH q < 0.05水平下具有差异性；其中676、589和452个还显示出至少两倍的差异。单独的基因过代表性和定量代谢物成员图谱检验支持了涉及M0的两个比较中的氧化磷酸化，以及M0与M6比较中的嘌呤代谢。丙氨酸、天冬氨酸和谷氨酸代谢显示出重复出现的代谢物图谱差异，但没有显著的基因过代表性。在1860个匹配的基因-特征对中，977个合并关联通过了BH q < 0.05的检验。转录本资源包括4642个新候选转录本和14个独特的差异剪接候选事件。这些发现描述了与年龄相关的转录和代谢变异，并为后续研究优先提供了通路和计算候选对象。

研究背景与目的

牦牛是青藏高原最具代表性的家畜，适应低氧寒冷环境，其肉以高蛋白含量和独特风味著称，在高原牧区经济中扮演重要角色。骨骼肌约占体重的40%，其发育不仅决定产肉量，还通过肌纤维类型和肌内脂肪沉积影响肉质。已有研究利用全转录组分析揭示了牦牛肉与脂肪组织中与肌内脂肪沉积相关的ceRNA调控网络，也有研究通过转录组和DNA甲基化分析鉴定了无角牦牛不同发育阶段背最长肌中与肌肉生长相关的差异表达基因，如TMEM8C、IGF2和MUSTN1。然而，这些研究主要关注单一组学层面，尚缺乏将全长转录组与代谢组联合分析以全面揭示年龄相关分子变化的报道。此外，现有研究多采用第二代测序技术，难以准确捕获全长转录本信息和可变剪接事件。因此，本研究旨在通过整合ONT全长转录组测序和非靶向LC-MS代谢组学，系统描绘雌性阿什丹牦牛三个关键出生后年龄点（10日龄、6月龄、54月龄）背最长肌的分子图谱，鉴定与年龄相关的差异表达基因、代谢物及潜在的关键通路和调控网络，为牦牛肉品质改良提供分子基础。

研究设计与主要结论

研究人员选取了27只雌性阿什丹牦牛，分别在10日龄（M0）、6月龄（M6）和54月龄（M54）三个时间点各取9只独立个体，采集其背最长肌组织。同时进行ONT全长转录组测序和非靶向LC-MS代谢组学分析。研究发现，M0与M54之间的差异表达基因数量最多（3480个）。蛋白质互作网络分析识别出高连通性候选基因，如ENSBGRG00000003669（在涉及M0的比较中）和MRTO4、BYSL（在M6 vs M54中）。代谢组分析显示大量差异丰度峰特征。跨组学通路分析表明，氧化磷酸化在涉及M0的两个比较中均获得基因和代谢物层面的共同支持，嘌呤代谢在M0 vs M6中获得双重支持。此外，研究还发现了4642个新转录本候选和14个独特的差异剪接候选事件。这些结果为理解牦牛肌肉发育的年龄相关分子机制提供了重要资源，并筛选出可用于后续验证的通路和候选分子。

主要关键技术方法

1. 1.
   样本来源与分组：所有27只雌性阿什丹牦牛均来自青海省大通牦牛场，在2026年4月于自由放牧管理条件下采样，分为10日龄（M0）、6月龄（M6）和54月龄（M54）三组，每组9只独立动物。
2. 2.
   ONT全长转录组测序与分析：提取背最长肌总RNA，利用Oxford Nanopore Technologies平台进行全长转录组测序。原始信号经Dorado转换为FASTQ序列，经Pychopper和NanoFilt处理后，使用minimap2比对至牦牛参考基因组（LU_Bosgru_v3.0）。采用Pinfish和StringTie进行转录本组装，gffcompare进行新转录本注释，TransDecoder预测编码区。利用SUPPA进行差异剪接分析。
3. 3.
   非靶向LC-MS代谢组学分析：采用超高效液相色谱-四极杆静电场轨道阱质谱联用仪进行数据采集，使用Waters ACQUITY UPLC BEH Amide色谱柱分离。原始数据经ProteoWizard转换为mzXML格式，并通过内部R流程进行峰检测和代谢物注释。
4. 4.
   差异表达与功能富集分析：使用PyDESeq2（基于DESeq2流程）进行基因和转录本的差异表达分析，阈值设为BH q < 0.05且|log2FC| > 1。采用超几何检验进行KEGG和GO功能过代表性分析。利用STRING v12.0构建差异表达基因编码蛋白的互作网络。
5. 5.
   整合分析与统计验证：通过配对Spearman相关分析评估基因与代谢物的关联。采用正交偏最小二乘判别分析和置换检验验证代谢组模型的可靠性。通过定量代谢物成员图谱检验评估代谢通路的整体变化。

研究结果

3.1 全长转录组图谱与转录本资源

• •
  3.1.1 样本结构与新转录本支持：主成分分析显示三个年龄组的基因表达谱明显分离，PC1和PC2分别解释了50.81%和21.27%的方差。共获得4642个新转录本候选，其中3573个具有功能注释，3875个在至少3个动物中估计计数≥10。
• •
  3.1.2 可变剪接候选与SSR计数单位：SUPPA分析共生成23733个结构剪接事件，经多重检验校正后，在M0 vs M6、M0 vs M54和M6 vs M54中分别鉴定出3、8和7个差异剪接候选，对应14个独特事件。
• •
  3.1.3 转录因子注释与lncRNA预测：提供了631个共有lncRNA预测和543个转录因子注释条目，作为后续研究的计算资源。
• •
  3.1.4 差异基因表达与下游分析：在M0 vs M6、M0 vs M54和M6 vs M54中分别鉴定出1999、3480和2350个差异表达基因。功能过代表性分析分别识别出32、31和18个KEGG通路。蛋白质互作网络分析显示，ENSBGRG00000003669在涉及M0的比较中连通性最高，而MRTO4和BYSL在M6 vs M54中共享最高连通性。
• •
  3.1.5 差异转录本表达与基因映射富集：分别鉴定出2499、4261和2868个差异表达转录本，映射到基因后的富集分析分别得到36、27和9个KEGG通路。

3.2 LC-MS代谢组学图谱与年龄相关丰度差异

• •
  3.2.1 保留峰特征与总体样本图谱：经质控后保留37901个非标准峰特征，其中1842个具有注释。主成分分析显示三个年龄组代谢物谱分离，PC1和PC2解释了58.48%的方差。
• •
  3.2.2 内部交叉验证与标签置换结果：OPLS-DA模型的三折交叉验证平均Q2值分别为0.9992、0.9981和0.9975，标签置换检验均显著（p = 0.001）。
• •
  3.2.3 大效应峰与描述性丰度聚类：在具有两倍以上差异的大效应子集中，M0 vs M6、M0 vs M54和M6 vs M54分别包含676、589和452个峰特征。K-means聚类将1842个注释峰分为3个簇，展示了不同的年龄相关丰度模式。
• •
  3.2.4 跨转录组和代谢组层面的通路证据：氧化磷酸化在M0 vs M6和M0 vs M54中同时获得基因过代表性和代谢物成员图谱检验的支持。嘌呤代谢仅在M0 vs M6中获得双重支持。丙氨酸、天冬氨酸和谷氨酸代谢在所有比较中均显示代谢物图谱关联，但未达到基因过代表性显著性。

3.3 匹配的基因-代谢物关联

在氧化磷酸化、嘌呤代谢和丙氨酸、天冬氨酸和谷氨酸代谢三个焦点通路中，共定义1860个独特的基因-特征对。合并Spearman分析鉴定出977个通过多重检验校正的关联，其中232个|rho| ≥ 0.7。但在调整年龄组后进行年龄内置换检验，无任何一对通过校正（最小校正q值为0.11904），表明合并关联主要反映了年龄组间的差异。

讨论与结论

讨论总结：本研究系统描述了雌性阿什丹牦牛三个出生后年龄点的背最长肌分子图谱。最大的差异表达基因和转录本集出现在M0与M54之间，但由于年龄间隔不等且为独立个体，这并不证明转录变化随年龄递增。氧化磷酸化和肌肉细胞骨架在涉及M0的比较中过代表，而PPAR信号在M6 vs M54中得到支持，这与不同年龄阶段的肌肉结构和代谢需求差异一致。研究提供了631个共有lncRNA预测和543个转录因子注释条目作为潜在调控背景，但并非每个条目都参与年龄依赖的肌肉网络。ONT测序获得了4642个新转录本候选，但读段支持、完整边界和生物学功能是不同的概念。差异剪接分析保留了少量候选事件，表明转录多样性远大于有统计学支持的差异剪接子集。蛋白质互作网络分析鉴定的高连通性候选基因，如RPS2、RPS9、MRTO4和BYSL，可能涉及蛋白质合成机器，但这只是关于翻译能力的假设而非证据。代谢组学分析显示大量差异峰特征，但氨基酸和核苷酸代谢提供了生物学背景。氧化磷酸化获得双重支持，但其5个测量化合物谱仅覆盖该通路的一部分。PPAR基因过代表性与年龄相关的脂质过程相容，但富集不代表通路激活或脂肪沉积量。匹配的基因-代谢物关联分析显示，合并关联主要源于年龄组间差异，而非直接的调控联系。

结论翻译：全长转录组学和非靶向代谢组学分析鉴定了雌性阿什丹牦牛在10日龄、6月龄和54月龄时背最长肌中与年龄相关的分子差异。最大的差异表达基因和转录本集出现在M0与M54之间。由差异表达基因衍生的蛋白质关联网络优先识别出高连通性候选基因，包括涉及M0的比较中的ENSBGRG00000003669，以及M6 vs M54中的MRTO4和BYSL。独立的基因和代谢物图谱分析支持了涉及M0的两个比较中的氧化磷酸化以及M0 vs M6中的嘌呤代谢。连同4642个新转录本候选和14个独特的差异剪接候选，这些结果为比较牦牛肌肉生物学提供了资源。它们确定了用于后续研究的通路和计算候选对象，但并未建立连续的发育机制、因果调控或对肉质的影响。本横断面研究比较了雌性阿什丹牦牛三个相距甚远的年龄点，未包含中间时间点或相应的肌肉生长、肌内脂肪、肌纤维特性或肉质测量指标。饮食、摄入量以及繁殖或其他生理状态无法完全与年龄区分。测序批次和LC-MS进样顺序记录缺失，无法评估或调整这些技术因素。转录本边界、优先基因和推断的分子关系未经qRT-PCR、蛋白质组学或功能实验独立验证。额外的审稿后通路和关联分析属于探索性研究，这些发现应被解释为在所采样条件下的年龄相关图谱。