基于生物信息学方法的小鼠精子发生关键基因筛选

李 欣 ，田 佳，脱征军，邵怀峰，张伟新，李委奇，王 琨，徐 寒，温 万*

（1.中国农业科学院北京畜牧兽医研究所，北京 100193；2.宁夏畜牧工作站，宁夏银川 750105）

转录组是某个物种或特定细胞在某一功能状态下产生的所有RNA 总和，通过对不同样本的转录组差异研究，不仅可以了解特定时间和空间条件下疾病的发生发展状态，也能作为研究疾病诊断和预后标志物的重要研究手段[1]。目前，RNA-seq 已经成为转录组研究中最为常见的手段，来自不同动物的多个器官的大量RNA-seq数据集陆续出现[1-2]，这些数据保存于NCBI 数据库中。精原细胞处于雄性生殖细胞的早期发育阶段，能不断进行有丝分裂，分化为精母细胞。精母细胞经减数分裂成为单倍体的精细胞。单倍体的精细胞经历形态变化、细胞器重排以及代谢途径和方式的改变等，最终形成终末分化的成熟精子，携带着由基因和环境共同决定的遗传信息。研究显示，Sycp3 与精子减数分裂染色体重组过程密切相关[3]，Crem 在精母细胞到圆形精子的转换过程起到开关的作用[4]。Zuo 等[5]对小鼠精子变形过程中的圆形精细胞、长形精细胞及附睾尾精子进行了RNAseq 和生物学分析，获得了大量选择性剪接事件、差异表达基因和新转录本信息。虽然对精子发生已有大量研究，但是从精原细胞依次到成熟精子的整体研究很少，且整个过程的分子机制仍不清楚。模式动物小鼠容易饲养，繁殖率高，遗传上有较高的纯合度，代谢类型、生理疾病等与人类接近，也具有完整的精子变形过程，样本取材方便。因此，本研究通过整合公共数据库中的转录组数据，将小鼠精原细胞和变形期精细胞和成熟精子进行转录本比较，以期全面揭示小鼠精子发生过程中基因表达的动态性与时空特异性，为后续全面研究精子发生调控机制提供理论基础。

1 材料和方法

1.1 数据来源 测序数据来自于GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/），登录号为GSE43717。小鼠转录组原始数据信息如表1 所示。

1.2 实验方法

1.2.1 数据预处理及差异表达基因筛选 将测序得到的原始序列（Raw reads）进行过滤去杂，去除测序接头序列、重复冗余序列和低质量序列，获得高质量序列数据（clean reads）。然后用TopHat 2 将各样品质控后得到的clean reads 与小鼠参考基因组（GRCm38.p6）进行比对，比对过程中允许有2 个碱基的错配。利用FPKM 值反映对应基因的表达量，使用R 软件的Ballgown 包进行差异分析。差异表达基因的筛选首先用取对数后的倍数变化（Fold change,FC）和统计检验P值初选，然后用FDR（False Discovery Ratio）校验方法对P值进行假阳性检验。当FDR≤0.05 且|log2(FC)|＞1 时，则为差异表达基因。

表1 小鼠转录组原始数据信息

1.2.2 差异基因的GO 分析和KEGG 富集分析 通过WebGestalt（WEB-based GEne SeT AnaLysis Toolkit，http://www.webgestalt.org/）分析工具对筛选出的差异表达基因进行GO 和KEGG 富集分析。

1.2.3 关键基因筛选 筛选出关键通路，将所包含的基因进行FPKM 值排序，取前18 名的基因进行详细的功能分析。

2 结果

2.1 差异表达基因筛选结果 对雄性小鼠不同阶段生殖细胞的转录组数据进行预处理，并通过阈值筛选得到2个不同阶段的差异基因，从精原细胞到圆形精细胞共获得10 785 个差异表达基因，其中上调基因6 176 个，下调基因4 609 个（图1-A）；从圆形精细胞到成熟精子共获得8 545 个基因，其中上调基因3 770 个，下调基因4 775 个（图1-B）。

2.2 差异基因通路富集分析 由图2 可知，精原细胞到圆形精细胞过程的差异基因主要参与了生物合成过程（生物过程）、质膜和囊泡部分（细胞组分）及阴离子配位（细胞功能）等过程。由图3 可知，圆形精细胞到成熟精子的差异基因主要参与了细胞器组织和芳香族化合物生物合成（生物过程）、水解酶活性（细胞功能）等过程，同时，可以明显看出圆形精细胞变为成熟精子的基因多集中于生物过程。

图1 差异基因火山图

精原细胞到圆形精细胞的差异基因，富集的KEGG信号通路主要是肌动蛋白细胞骨架的调节和丝裂原活化蛋白激酶（Mitogen-activated Protein Kinase,MAPK）信号通路，而圆形精细胞到成熟精子上调差异基因主要是紧密连接和MAPK 信号通路。

2.3 关键基因的筛选 鉴于GO 和KEGG 富集分析获得的条目较多，本研究根据精子发生过程的关键点筛选出相关性较强的8 个通路，并对这些基因进行了查重去冗余，获得628 个差异基因；再将638 个基因在3 个阶段的FPKM 值进行排序，取FPKM 值最大的前18 个基因进行研究，该18 个基因所属的GO 条目详细信息见表2 和图4。表2 所示的7 个通路是精子发生过程的关键条目。该条目中，Reproduction（GO:0000003）和Reproductive process（GO:0022414）条目中的基因可能与繁殖相关；Motile cilium（GO:0031514）、Sperm flagellum（GO:0036126）和Cilium（GO:0005929）条目中的基因可能与精子尾部形成相关；Golgi apparatus（GO:0005794）条目中的基因可能与精子顶体的形成相关；Spermatogenesis（GO:0007283）条目中的基因可能与精子发生相关。

图2 精原细胞到圆形精细胞过程的上下调基因GO/KEGG 富集分析

图3 圆形精细胞到成熟精子过程的上下调基因GO/KEGG 富集分析

表2 精子发生过程的关键条目

3 讨 论

本研究通过分析GSE43717 转录组测序数据，从精原细胞到圆形精细胞筛选到的上调基因比下调基因多1 576 个，而从圆形精细胞到成熟精子筛选到的下调基因比上调基因多1 005 个，这说明在早期单倍体细胞、圆形精子细胞中合成了许多新的mRNA，但在晚期单倍体细胞、长形精细胞中，染色质固缩进而导致结构的改变使转录失活[7]。通过GO 分析发现差异基因主要参与精子质膜和囊泡形成、生物合成、水解酶活性过程。差异表达基因主要富集的KEGG 信号通路主要是肌动蛋白细胞骨架的调节、紧密连接和MAPK 信号传导通路。其中MAPK 信号传导通路是贯穿精原细胞到精子整个过程的共有部分，该通路不但参与了精子生成与凋亡、运动、获能和顶体反应等过程，还是附睾中精子成熟的重要通路。这些重要生物学过程以及生物通路的发现，是根据精子发生过程关键点的转录组测序、差异基因和FPKM 值较大基因筛选而来的，可见在研究生物发育过程中，选择合适的时间点、差异基因筛选标准以及基因表达量对机理的研究至关重要。

图4 关键基因的热图分析

3.1 参与精子发生相关基因Prm2、Selenof、Tnp1、Tnp2、Ropn1l、Pafah1b2、Spink2、Prm1参与精子发生过程。在精子发生过程中，mRNA 的定位和翻译被认为是以一种阶段性的方式调节的。有报道显示，Prm2mRNA 位于圆形精细胞的拟染色体，在长形精细胞的胞浆进行翻译，说明Prm2在整个精子发生过程中具有时空调节的特性[8]。而且Prm1、Prm2、Tnp1 和Tnp2均是在染色质重塑过程中取代组蛋白的基本染色体蛋白质，它们在精子变形前期被转录并以惰性状态存在，并在长形精子细胞中被激活翻译[9-10]。丝氨酸蛋白酶抑制剂Kazal 2 型（Spink2）基因属于SPINK 家族。研究表明，Spink2在小鼠睾丸中表达量较高，粗线期精母细胞阶段首次显著表达，而在附睾中表达量显著降低[11]。Ropn1l是构成精子鞭毛纤维鞘的关键基因之一，直接影响到精子活力和生育力[12]，这也是精子发生过程的重要一环。另外，Selenof可能参与二硫键形成相关的氧化还原反应，质控经内质网折叠所产生的蛋白质[13-14]，进而影响精子发生过程。

3.2 与精子尾部形成相关基因 精子尾部结构与精子运动功能密切相关，具有正常尾部结构的精子才能游动，尾部每一个部分的发育均受到相关基因的精确调控，如果发生基因缺失或突变会导致相应编码蛋白的改变，从而导致精子尾部结构异常或产生功能障碍[15]。Ldhc、Hsp90aa1、Gstm5、Csnk2b、Spa17均与精子尾部形成相关。Hsp90aa1 属于分子伴侣，参与ATP 酶活性相关的功能循环，从而诱导目的蛋白的构象变化，导致其活化[16]。Hsp90aa1 位于人精子颈部、尾部中段等区域，在精子获能时表达量增加，在胞内钙稳态以及蛋白酪氨酸磷酸化方面发挥重要作用；在氧化应激条件下，可通过激活一氧化氮合酶来保持精子活力[17]。Gstm5 属于谷胱甘肽（GST）家族，在纤维鞘中表达，促进生殖细胞的增殖和分化，同时，在氧化应激过程中可能发生改变，以对抗自由基对睾丸和生精细胞的攻击[18]。Spa17是一种在睾丸组织中特异性表达的蛋白，定位于初级精母细胞、次级精母细胞和成熟精子中，而精原细胞、支持细胞和间质细胞中无表达[19]。成熟精子中Spa17 主要表达于精子尾部[20]，推测其与精子运动有关。另外，Ldhc和Csnk2b未见在精子中的相关研究。

3.3 与繁殖相关的基因 与繁殖相关的基因包括Txnrd3、Tcp1、Gpx4、Dkkl1和Dbil5。动物硫氧还蛋白还原酶Txnrd 是还原型辅酶Ⅱ（NADPH）依赖性蛋白质，属于吡啶核苷酸二硫化物氧化还原酶家族，含有黄素腺嘌呤二核苷酸（FAD），已知的有Txnrd1、Txnrd2 和Txnrd3。其中，由Txnrd3基因编码的同工酶主要存在于睾丸和发育早期的精子细胞中，在结构上有一个谷胱甘肽还原型结构域，推测它参与精子成熟过程中的氧化还原反应[21]。Tcp1 即t 复合多肽1，在精子发生过程中表达量较高，定位于精子顶体囊泡中，影响精细胞形成过程中顶体的形成[22]。Gpx4 即谷胱甘肽过氧化物酶4，是一种膜相关的谷胱甘肽过氧化物酶，在哺乳动物睾丸中表达量较高，其亚型均参与精子发育[23]，线粒体Gpx4（mGpx4）敲除小鼠精子质量下降，进而导致不育，该基因的所有亚型全部敲除后引起胚胎致死[24]。Dkkl1的mRNA 主要产生于精母细胞，其蛋白定位于顶体，可促进精子穿透透明带[25]。而且，Dkkl1在弱精子症患者的精子中表达明显降低[26]，推测其在精子运动中起重要作用。另外，研究显示，Dbil5 是一种睾丸特异性蛋白，在精子发生的后期表达[27]，大鼠Dbil5 的mRNAs 翻译在精细胞早期受到强烈抑制，在精细胞后期被激活[28]，推测其可能通过影响精子后期的变形，进而影响繁殖。

4 结 论

本研究通过对小鼠精子形成期3 个阶段的差异表达基因进行分析，发现FPKM 值靠前的差异基因主要集中在精子发生、精子尾部形成以及繁殖等条目。尤其与精子发生相关，包含鱼精蛋白在内的8 个基因是学者研究比较全面的，体现了这些基因的重要性。此外，对繁殖相关的基因研究也较多，主要关系到精子受精能力引起的不育现象。总之，精子形成期间发生了大量基因的变动，关注期间关键基因可为精子变形机制的阐明提供理论基础。