弱精子症相关microRNAs和基因的生物信息学研究

       男性不育症是人类生殖健康面临的一个严重问题［1］。其中，弱精子症被认为是男性不育的常见原因。前向运动精子百分率低于32%即为弱精子症［2］，约40%的不育男性患有弱精子症，其主要特点是精子活力低下［3］。弱精子症由多种因素引起，然而导致弱精子症的具体分子机制尚未完全阐明。据报道，导致男性不育的遗传因素越来越多。有研究发现，CABYR和ROPN1在弱精子症患者样本中显示下调，并且两个基因的表达呈正相关，表明两个基因的共表达可能是正常精子活力的先决条件［4］。最近国内研究发现，睾丸优势表达的腺苷酸激酶家族亚组Ⅲ的AK8和AK9调控精子鞭毛的能量转运和运动功能及机制［5］。同时，研究发现精液样本中也存在大量microRNAs(miRNAs)，其表达水平的改变与精液参数相关，其中包括精子活力低下［6］。本研究通过分析已发表的研究和GEO数据库的数据，探讨了弱精子症患者的差异表达 miRNAs(DEM)和差异表达基因(DEGs)。然后，通过TargetScan预测DEMs的靶基因，再对上述的DEGs和DEMs的靶基因取交集获得最终的DEGs。笔者使用DAVID对DEGs进行了GO和KEGG分析。通过STRING构建蛋白质-蛋白质相互作用（PPI）网络，对弱精子症中DEM和DEG的相关生物学特征和分子途径提出了新的见解。miRNA-mRNA 网络进一步展示了与弱精子症相关的潜在机制。

1 资料与方法

1.1 数据检索 

1.1.1 miRNA相关数据集和文献检索     首先在GEO数据库中检索已发表的精子miRNA相关文献，未检索到可利用的数据集。然后在PubMed和Web of science上进行了系统搜索，结果发现4项符合要求的研究［7-10］。该4项研究中的样本均为精液样本，且均以生育力正常的育龄男性为对照。本研究中所用数据库及检索筛选流程见图1。

1.1.2 mRNAs相关数据集的检索     利用GEO数据库中的关键词“弱精子症（asthenospermia）”和“智人（Homo sapiens）”搜索mRNA表达数据集。通过检索，下载了数据集GSE92578，该数据集由Zhang等提供，GPL20301平台, Illumina HiSeqTM 4000测序，包括弱精子症组3个样本和正常对照组4个样本，共7个样本。

1.2 差异表达miRNAs的筛选

     对上述4项研究进行了汇总分析，然后通过联川生物云平台的韦恩图绘制工具绘制韦恩图。最后经过综合分析确定研究中差异表达的miRNAs作为DEMs。

1.3 差异表达基因筛选

      ①通过TargetScan数据库预测DEMs的靶基因。②通过GEO2R工具，对GSE92578数据集进行分析，设定差异表达基因的筛选条件为FDR<0.05和|log2FC|>1。③通过联川生物云平台的韦恩图绘制工具绘制韦恩图并确定最终的DEGs。

1.4 DEGs的GO和KEGG富集分析

      在DAVID数据库进行GO和KEGG富集分析确定DEGs功能。其中，物种选择为智人（Homo sapiens）,P<0.05表示有统计学意义。

1.5 DEGs的互作网络和枢纽基因的鉴定

      通过基因互作检索工具STRING 构建DEGs的蛋白质互作网络,同时设定PPI评分≥0.6。利用Cytoscape软件对源自STRING的PPI可视化，并利用CytoHubba插件对PPI中的节点基因进行排序，其中10位基因被认为是枢纽基因,P<0.05表示有统计学意义。

2 结果

2.1 弱精子症相关的差异表达miRNAs

     对上述4项研究进行了汇总分析（表1）并绘制韦恩图（图2，2215页），确定7个差异表达的 miRNAs即miR-27a-3p、miR-374b-5p、miR-193a-3p、miR-1-3p、miR-34b-3p、miR-26b-5p、miR-143-3p。

2.2 弱精子症相关的差异表达mRNA

      对上述4项研究确定7个差异表达的 miRNAs经TargetScan数据库进行靶基因预测后共获得9847个靶基因。在GSE92578数据集共筛选出330个差异表达基因。对上述两方面获取的基因经韦恩图取交集共获取84个DEGs（图3，2215页）。

2.3 差异表达基因的GO和KEGG分析结果

      对上述84个差异表达基因行GO和KEGG富集分析（图4，2215页）。

2.4 差异表达基因的PPI网络和枢纽基因

     利用STRING数据库构建差异表达基因的PPI网络（图5，2216页）。利用Cytoscape中的cytoHubba插件筛选出排名前10的枢纽基因（图6，2216页）。

3 讨论

   弱精子症以精子活力低下为特征，然而，影响弱精子症的遗传特征尚不清楚，需要进一步研究。有研究报道，精子中RNA包括多种mRNA、miRNA和lncRNA，是调节弱精子症的关键角色［11］。本研究通过生物信息学方法研究弱精子症中差异表达的miRNA和基因，探索弱精子症的潜在调控网络。

     本研究通过分析筛选以往的研究，确定了7个差异表达的miRNAs，包括 miR-27a-3p、miR-374b-5p、miR-193a-3p、miR-1-3p、miR-34b-3p、miR-26b-5p、miR-143-3p。富含半胱氨酸的分泌蛋白2（CRISP2）是精子中的一种重要蛋白质，在调节精子鞭毛运动中发挥作用，缺乏CRISP2的精子会表现出精子活力低下。有研究发现，miR-27a-3p通过与CRISP2的3' 非翻译区 (3'-UTR) 结合特异性靶向 CRISP2，从而抑制 CRISP2 转录后表达［12］。有研究报道miR-374b-5p在睾丸组织中表达并发挥重要作用，其表达减少可以作为精子DNA碎片指数增加的指标,表明miR-374b-5p有可能成为特发性不育男性的诊断生物标志物［13-14］。对不育三倍体鱼的miRNA的研究发现miR-193a-3p表达上调，可能通过靶向功能基因参与精子活动和睾丸发育［15］。有研究发现miR-34b家族在精子发生过程中发挥着重要作用［16］。此外，不育男性中miR-34b启动子区的甲基化频率高于生育男性，且在弱精子症患者中甲基化频率最高［17］。在关于睾丸miRNA在热应激诱导的精子发生障碍的潜在功能研究中，研究者通过功能富集分析、共表达调控网络、相关性分析和体外实验，发现miR-143-3p可能是热应激下影响精子发生的代表性潜在关键调控因子［18］。虽然未见miR-1-3p、miR-26b-5p与弱精子症相关的研究报道，但为下一步的研究提供了方向。

     本研究共获得弱精子症相关的10个hub基因，分别为AKT1、MAPK3、BRD4、DNMT3A、FURIN、LMNB2、COL5A2、COL5A3、COL11A1、COL27A1。AKT1是与各种细胞过程相关的丝氨酸/苏氨酸蛋白激酶基因AKT的3种不同亚型之一。在一项关于AKT1和AKT2对生殖影响的研究中发现，AKT1表达异常的小鼠精子表现出葡萄糖转运蛋白表达减少、细胞凋亡增加［19］。MAPK3基因编码的蛋白质是MAP激酶家族的成员。其在信号级联中发挥作用，调节各种细胞过程，例如响应各种细胞外信号的增殖、分化和细胞周期进程［20］。在一项关于大鼠支持细胞增殖和凋亡调节的研究中，MAPK3/1在支持细胞和生殖细胞之间的黏附中发挥作用［20］。不同DNMT基因的变异与特发性男性不育症有不同的关系，相关研究发现在DNMT3A基因敲除的小鼠，其胚胎致死性增加、雄性和雌性生殖细胞印记丧失、精原阶段生精停滞或无精子症［21］。有研究者得出结论， COL11A1与纤毛运动有关，并参与精子运动［22］。此外，BRD4、FURIN、LMNB2和COL5A2等相关基因均未见与弱精子症的相关研究。

      在对靶基因取交集后获取的84个基因进行GO和KEGG分析。在GO分析中，细胞组份主要富集在细胞质、细胞膜、溶酶体膜和顶体等。生物学过程的富集分析主要体现在对RNA聚合酶Ⅱ启动子的转录进行正向调节，蛋白质转移，细胞外基质组织和磷酸化作用的正向调节等。分子功能富集主要在同源蛋白绑定，蛋白激酶绑定和蛋白质丝氨酸/苏氨酸激酶活动等。在KEGG分析中，主要富集在PI3K-Akt信号通路, 人类乳头瘤病毒感染通路, 血小板激活和化学致癌－受体激活等通路。已有研究证实PI3K-Akt信号通路在男性生殖过程中发挥重要作用［23-24］。

     本研究也存在一些局限性。首先，本研究的数据来源相对单一，数据量较小。其次，分析获得的miRNAs、基因和通路有待经实验的验证。同时，miRNA 或基因中是否存在共表达现象并对于精子活力产生下调还有待进一步研究。本研究的结果为今后弱精子症的机制研究提供理论基础，以此研发新的治疗改善精子活力的药物，将有助于推动男性不育治疗的发展。