2023年1月6日，Cell子刊《iScience》在线发表了广西壮族自治区水牛研究所邓廷贤副研究员、河南农业大学动物科技学院反刍动物生产科技创新团队付彤教授和刘深贺博士联合研究的重大科研成果。该研究采用Direct RNA测序（DRS）技术解析了公牛睾丸发育的转录组图谱，鉴定到不同发育阶段睾丸发育和精子发生的hub基因，强调了RNA加工和表达的复杂性以及参与睾丸发育和精子发生的转录调控。论文第一作者为刘深贺博士，第一通讯作者为邓廷贤博士，并列通讯作者为付彤教授。龙8总区基因转录组事业部李亚丽、杨柳等人担任共同作者。

西门塔尔牛以其生长速度快、肉质好和环境适应性强等特征闻名世界。在过去的二十年中，西门塔尔牛与一些中国地方肉牛品种进行杂交，以提高其生产效率。公牛睾丸大小的线性增加伴随着细胞区室的成比例增长对建立精子发生，早期成熟和潜在优质育种公牛等选择至关重要。

本研究旨在使用DRS技术探讨西门塔尔牛在不同发育阶段的睾丸转录组动态变化和复杂性。阐明了m5C介导的选择性剪切事件调节西门塔尔公牛睾丸发育和精子发生的潜在分子机制。该研究结果可为今后制定优良的公牛育种计划提供策略。

材料：9头临床健康的西门塔尔公牛。分为三组：青春期前（TY 0，n = 3，1周龄小牛）、青春期（TY 1，n = 3，约1岁公牛）和青春期后（TY 2，n = 3，约2岁公牛）。

主要研究方法：DRS

测序平台：Nanopore PromethION

1、 不同阶段睾丸发育形态学研究

对3个发育阶段的牛睾丸形态学特征（睾丸重量，TM；睾丸长度，TL；睾丸宽度，TW）进行检测，结果表明不同发育阶段的睾丸表型参数存在显著差异。公牛睾丸组织形态学研究表明，青春期前后睾丸组织形态差异显著：生精小管的直径随睾丸发育逐渐增大（图1）。

图1睾丸组织形态学和精液质量参数检测

2、 DRS数据的统计

DRS测序总产生了5077万条长reads。测序Reads平均长度分别为TY0、TY1和TY2的1663、1381.5和1382.83 bp（图2B）。和基因组比对后，所有样本中比对上的Reads长度分布为：TY2具有最高的中位长度(899个核苷酸)，其次是TY1(886个核苷酸)和TY0(717个核苷酸)(图2C)。

图2 DRS数据统计

3、 DRS揭示牛睾丸转录组特征

DRS测序reads共鉴定到56,560个牛的编码基因，包含97,802个转录本（28,999个已知转录本和68,803个新转录本)(图3A）。鉴定的新基因数目（45,072个）是已知基因（21,427）的2倍以上。已知转录本和新转录本的中位长度分别为2,304和836 bp（图3D）。对比NCBI中注释（Btau_5.0.1）的参考转录本，优化后转录本的中位长度为3,237 bp（图3G），表明牛睾丸转录本编码区（CDS）的实际长度长于Btau_5.0.1中的注释长度。

图3 牛睾丸转录组的复杂性特征

4、 睾丸发育期间DEG和DEI之间存在较大差异

PCA显示，两份样品（TY1A和TY2A）的生物重复性较差，剔除后进行后续分析（图4A）。任意两组之间比较共鉴定了15,636个DEG和23,932个DEI，其中9,641个基因在DEG和DEI之间重叠（图4B）。TY 0-TY 1和TY 0-TY 2之间共同表达的DEG和DEI占比为79.43%和78.87%。TY 0属于性成熟前的阶段，而TY 1和TY 2两组中的动物可以产生精子。由此可以推断，大多数DEG和DEI参与睾丸发育和精子发生。

图4睾丸发育过程中，DEGs和DEIs存在较大差异

5、 RNA加工事件对睾丸重量相关hub基因的影响

WGCNA分析结果显示DEG和DEI数据集分别检测到3个和2个共表达模块。所有模块均与睾丸形态和精液表型性状显著相关。进一步鉴定与TM相关的模块中的枢纽基因，DEG和DEI之间重叠的5,043个基因被最终定义为hub基因（图5B）。Poly(A)长度与这些hub基因的TPM呈显著负相关（图5D）。hub基因共鉴定了12,792个AS事件和2,795个DAS事件，最丰富的AS事件是SE（23.63%），最多的DAS事件是A3（图5E）。不同的DAS事件在TY0和TY1/TY2组之间的表达水平具有显著差异，但是TY1和TY2之间没有差异（图5F）。以上结果表明，枢纽基因可以通过不同的AS事件调节表达水平，从而影响睾丸重量。

图5 WGCNA分析及睾丸重量相关的hub基因的Poly(A)尾、选择性剪接分析

Hub基因中平均每个转录本的m6A和m5C修饰位点数分别为4.12和25.58，它们在终止密码子和3’UTR区域附近富集（图6G-H）。具有差异m6A和m5C位点的基因，占总hub基因的41.46%和28.71%，在不同比较组中具有不同的表达模式（图6J）。这些hub基因在DRS和MeRIP-Seq验证数据中观察到具有相似的表达趋势。

图6 睾丸重量相关的hub基因的甲基化修饰分析

6、 睾丸重量相关hub基因的功能分析

hub基因在有性生殖中富集最多，其次是雄性配子生成、精子发生、多生物生殖过程、鞭毛精子活力等（图7A）。对富集到Top10 GO terms的基因进行PPI分析，共有87个hub基因分布在3个PPI网络中（图7 B-C）。根据综合评分，筛选出与睾丸发育和精子发生相关的前20个hub基因进行qPCR验证，显示出相似的表达水平趋势（图7D）。

图7 Top20 hub基因的功能注释分析

7、 MAEL基因通过转录后修饰和剪接参与睾丸发育

Top20 hub基因的Poly(A)平均长度随睾丸发育而不断增加（图8A）。此外，m6A和m5C均有助于性成熟前后top20枢纽基因的表达水平（图8B）。然而，只有3个枢纽基因（MAEL、ROPN1L和SRATA6）在所有发育阶段显示发生DAS事件，MAEL显示最高发生率（图8C）。

MAEL的转录本在不同发育阶段的表达水平不同，其中rna3047的表达水平高于gene2544.t2（图8D）。共3种AS事件（包括A3、A5和SE）同时发生在2544.t2转录本，且在睾丸发育的不同阶段中，这些事件都有显著差异（图8E）。此外，研究者还发现这些DAS事件的剪接位点上存在丰富的m5C位点，且无论何种DAS发生，TY1和TY2组中的m5C位点数量均显著高于TY0组（图8E）。以上结果指出m5C位点可能促进DAS事件的发生，进而产生不同的转录本来增加MAEL基因的表达。

该研究首次公布了公牛睾丸转录组图谱，探究了m5C介导的选择性剪切事件调节西门塔尔公牛睾丸发育和精子发生的潜在分子机制，解码了睾丸发育和精子发生的枢纽基因。该研究结果可为今后制定优良的公牛育种计划提供策略。