不同阻滞阶段非梗阻性无精子症的生物信息学分析

王万伦，刘桐佳，张婷，李硕，边艳超，肖瑞

据世界卫生组织估计，目前全球约有10%～15%的夫妇患有不孕症，其中男性因素约占所有不孕症病因的一半[1-2]。无精子症是男性不育的原因之一，可分为两大类：输精管道机械性梗阻（梗阻性无精子症）和原发性生精衰竭[非梗阻性无精子症（nonobstructive azoospermia，NOA）]。约60%的无精子症患者最后被诊断为NOA[3]。NOA是男性不育中最严重的一种，表现为精子生成过程被破坏[4-5]；精子发生过程不活跃，生精细胞阻滞在不同发育阶段，如减数分裂前期阻滞阶段（premeiotic arrest stage，PRE）、减数分裂期阻滞阶段（meiotic arrest stage，MEI）、减数分裂后期阻滞阶段（postmeiotic arrest stage，POST）[6]。NOA可继发于多种疾病，包括隐睾、睾丸易位、染色体异常、Y染色体微缺失和其他导致睾丸功能衰竭的遗传疾病[7-9]，但其分子机制还未被阐明[5，10]。虽然显微解剖睾丸精子提取（micro-TESE）是NOA的标准治疗方法，但约50%的NOA患者的精子提取不成功[11]。因此，亟需阐明精子发生过程的明确分子机制，并发现有效的NOA诊断标志物或治疗靶点。二代测序技术的不断成熟有助于从基因层面进一步了解疾病的发生、发展过程。本研究通过对比20例不同阻滞阶段的NOA患者睾丸组织及4例正常生精男性睾丸组织的表达谱数据，筛选出差异基因（differentially expressed genes）及枢纽基因（Hub基因），并且预测其关键通路，以便为临床治疗男性不育寻找潜在的生物标志物和治疗靶点。

1 资料与方法

1.1 数据来源本研究分析的基因表达数据集来自GEO（Gene Expression Omnibus）数据库（https://www.ncbi.nlm.nih.gov/geo/）。从该数据库中共检索到22篇有关人类无精子症的系列数据。经过仔细审查，选择了1个基因表达谱（GSE45885）。此表达谱基于GPL6244平台[HuGene-1_0-st]Affymetrix Human Gene 1.0 ST Array[transcript(gene)version]，所有的数据都可以在网上免费获得，且该研究不涉及本文作者在人类或动物身上进行的任何实验。

1.2 差异基因数据处理GEO2R在线分析工具（https://www.ncbi.nlm.nih.gov/geo/geo2r/）检测NOA与正常样本之间的差异基因，计算调整后的P值和。符合临界值标准（P＜0.05、≥2）的基因被认为是差异基因。使用miRDB（http://mirdb.org）预测微小RNA（microRNA，miRNA）的靶基因，对每个数据集进行统计分析，并使用韦恩图线上工具（bioinformatics.psb.ugent.be/webtools/Venn/）识别交叉部分。

1.3 差异基因的GO（Gene Ontology）及KEGG（Kyoto Encyclopedia of Genes and Genomes）信号通路富集分析GO分析是大规模功能富集研究的常用方法。基因功能可分为生物过程（biological process）、分子功能（molecular function）和细胞成分（cellular components）。KEGG是一个广泛使用的数据库，其存储了大量关于基因组、生物途径、疾病、化学物质和药物的数据。本研究使用DAVID（Database for Annotation，Visualization and Integrated Discovery）工 具（https://david.ncifcrf.gov/）对差异基因进行GO富集分析和KEGG信号通路富集分析。设定P＜0.01、count＞5认为数据有统计学意义。

1.4 蛋白质相互作用（protein-protein interaction，PPI）网络图构建及关键基因鉴定本研究使用STRING数据库（http://string-db.org/）分析PPI网络。为评估潜在的PPI关系，此前确定的差异基因被映射到STRING数据库。提取PPI的最低要求交互评分为0.9。随后用Cytoscape软件（www.cytoscape.org/）对PPI网络进行可视化。连接度越高的节点对维持整个网络的稳定性越重要。CytoHubba是Cytoscape中的一个插件，用来计算每个蛋白节点的程度。在本研究中，前10名的基因被确定为枢纽基因，并列出PPI评分。

2 结果

2.1 确定差异基因GSE45885基因表达谱中包含4例正常生精样本，20例NOA患者样本，其中包括2例PRE样本，7例MEI样本，11例POST样本。数据标准化后（见图1），以矫正后P＜0.05、≥2为标准确定差异基因。基于该标准，PRE组筛选出463个上调基因，12个下调基因；MEI组只有2个下调基因，没有上调基因；POST组筛选出3个下调基因，没有上调基因（见图2A）。值得注意的是，本研究对3个阻滞阶段的NOA下调基因取交集，发现2个相同的下调基因：MIR15A和MIR509-3（见图2B）。用在线工具miRDB预测这2个miRNA的靶基因，发现MIR15A有1 415个靶基因，MIR509-3有282个靶基因，两者有58个共同靶基因（见图2C）。

图1 NOA数据归一化矫正后

图2 NOA数据集可视化处理

2.2 差异基因的功能富集分析使用DAVID对PRE中的上调基因以及MIR15A和MIR509-3的共同靶基因（见图3）进行GO功能富集分析和KEGG信号通路富集分析。GO富集分析显示PRE组中的上调基因主要富集在生物过程中，包括精子发生（spermatogenesis）、精子细胞发育（spermatid development）和精子的运动（sperm motility）等；而在细胞成分中，PRE组的上调基因主要富集在顶体囊泡组（acrosomal vesicle）、纤毛运动（motile cilium）、微管骨架（microtubule cytoskeleton）和微管（microtubule）形成等功能中（见图4）。MIR15A和MIR509-3共同的靶基因主要富集在生物过程中，包括通过质膜黏附分子黏附同质细胞（homophilic cell adhesion via plasma membrane adhesion molecules）、神经系统发育（nervous system development）、细胞黏附（cell adhesion）；在细胞成分和分子功能中也有一定程度富集，分别为质膜的组成部分（integral component of plasma membrane）和钙离子结合（calciumionbinding）。本研究进一步对这些共同靶基因进行KEGG信号通路富集分析发现，这些基因主要富集在磷酸肌醇代谢（inositolphosphatemetabolism）以及心肌细胞肾上腺素能信号传导（adrenergic signaling in cardiomyocytes）通路中（见表1）。

表1 共同靶基因的GO及KEGG信号通路富集分析

图3 MIR15A和MIR509-3的共同靶基因

图4 PRE组上调基因GO富集分析

2.3 PPI网络构建及关键基因筛选使用STRING数据库预测PRE组上调基因的PPI网络，然后使用Cytoscape软件v3.7.1（https:/cytoscape.org/）可视化来自STRING数据库的PPI网络（见图5）。利用Cytoscape中的CytoHubba插件，根据程度（Degree）计算方法对PPI网络中的节点进行排序，其中前10个基因被认为是枢纽基因，分别为DNAI1（10分）、DNAI2（9分）、PGK2（7分）、ROPN1L（7分）、ARMC4（7分）、DRC1（7分）、DNAAF3（7分）、CABYR（6分）、ZMYND10（6分）和CCDC65（6分）。

图5 PRE组上调基因的PPI网络

3 讨论

3.1 精子发生阻滞的时期可能与差异基因的数量有关本研究从GEO数据库中筛选GSE45885数据集，并利用生物信息学方法鉴定该数据集中3个阻滞阶段NOA的差异基因，共筛选出475个差异基因，PRE阶段的差异基因数量最多，MEI和POST阶段差异基因的数量明显减少，提示精子发生阻滞的阶段可能与差异基因的数量有关联，差异基因的数量越多，精子发生阻滞的阶段越靠前。本研究对PRE中上调的差异基因做GO富集分析发现，在PRE中上调基因主要参与生精过程中极其重要的多个环节，这些基因功能异常可能会导致精子在减数分裂前期发生阻滞；PRE表达上调的差异基因中筛选出10个枢纽基因，检索文献发现DRC1是纤毛和鞭毛中微管连接蛋白与动力蛋白调控复合体（nexin-dyneinregulatorycomplex）组装的重要调控因子。DRC1缺失将导致纤毛缩短和运动障碍[12]。也有病例报道ARMC4缺失会导致精子鞭毛异常，轴突组织紊乱，外动力蛋白臂（outer dynein arm）缺失[13]，最终导致男性不育。可见以上2个基因缺失会引起生殖异常，而本研究中两者表达均上调，推测DRC1和ARMC4表达异常增加可能也会导致男性不育；可将DRC1和ARMC4作为NOA研究的潜在靶基因，进行后续的功能研究。

3.2 MIR15A和MIR509-3表达下降可能是精子发生阻滞的启动因素本研究发现3个阻滞阶段NOA的下调基因都为miRNA，并且MIR15A和MIR509-3为3个阻滞阶段NOA的共同下调基因。miRNA是由18～25个核苷酸组成的小单链RNA，miRNA缺乏典型的开放阅读框，不能翻译成蛋白质，却在细胞的生长发育、蛋白质分泌和基因调控中发挥关键作用；在生殖系统中，miRNA也发挥着极其重要的作用[14-15]。某些miRNA下调可能会导致男性不育。有研究表明MIR15A参与了多种卵巢功能的调控，如卵巢颗粒细胞的增殖与凋亡、孕激素、雄激素和雌激素的释放[16]，并且可以抑制前列腺癌的发生[17]，但是在男性不育中暂未见文献报道。有研究表明NOA和少精子症患者的MIR509表达水平比可生育对照组低[18]。本研究发现，在3个阻滞阶段的NOA中，MIR15A和MIR509-3表达均下调，提示MIR15A和MIR509-3可能对精子发生的进程具有极其重要的调节作用。本研究还发现，MIR15A和MIR509-3有58个共同的靶基因，这些共同靶基因主要富集在神经系统发育、细胞黏附和钙离子结合等功能中。已有研究证明，神经系统与生殖系统之间有着密切的关联[19]，某些细胞黏附分子的缺乏也会导致雄性不育，如细胞黏附分子1[20-21]。提示MIR15A和MIR509-3可能通过这些共同靶基因协同调控精子发生的进程，两者表达下调可能启动精子发生阻滞，这有待于进一步的实验研究验证。

总之，本研究对3个阻滞阶段NOA患者与正常男性睾丸组织之间的差异基因进行了生物信息学整合分析，发现3个阻滞阶段的NOA中差异基因的数量有明显差异，但是3个阻滞阶段的NOA又具有共同的下调基因，这些基因可作为后续NOA发病机制研究的潜在靶基因。