bta-miR-34b/c和bta-miR-449a/b/c靶基因预测及生物信息学分析

杜鑫泽，马鑫浩，张殿琦，杜嘉伟，马 婧，谢琨成，何 杰，昝林森，2

(1.西北农林科技大学动物科技学院，杨凌 712100；2.国家肉牛改良中心，杨凌 712100；3.西安市奶牛育种中心，西安 710000)

在雄性哺乳动物中，精子发生是一个复杂且受多种途径精密调控的生理过程，主要经历4个阶段：第一阶段，A0型精原细胞有丝分裂经历A1-A4精原细胞、中间型精原细胞、B型精原细胞并最终形成初级精母细胞；第二阶段，初级精母细胞开始第一次减数分裂形成次级精母细胞；第三阶段，初级精母细胞进行第二次减数分裂形成圆形精子细胞；第四阶段，圆形精子细胞经过一系列变形分化为成熟精子[1]。在精子发生的整个过程中，雄性生殖细胞分化过程具有复杂的转录调控机制，除了编码蛋白质的mRNA外，许多非编码RNA如微小RNA(microRNA,miRNA)、piwi相互作用RNA(PIWI interacting RNA,piRNA)、环状RNA(circularRNA,circRNA)、长链非编码RNA(long noncoding RNA,lncRNA)等都发挥着重要的转录后调控作用。

miRNA是一类进化上高度保守的长度在22 nt左右的小分子非编码RNA，miRNA可以参与调控包括性成熟、睾丸发育、精子发生在内的诸多生物学过程，其在精子发生过程中以细胞特异性或阶段特异性方式表达[2]。靶向mRNA的3′-UTR以实现基因沉默是miRNA的经典作用机制[3]，miRNA常通过与靶向mRNA的3′-UTR中的互补序列结合，并通过与多蛋白RNA诱导的沉默复合物(RNA-induced silencing complex,RISC)相互作用来加速转录物衰减或翻译抑制来抑制基因表达[4]。miRNA在众多生物学途径中的作用突出了转录后基因调控的复杂性。研究表明，miRNA的表达具有较强的组织特异性，表明miRNA在不同组织、不同生物学过程中具有独特的功能。miRNA因其在基因表达的转录后调控中占据主导地位而被深入研究并广泛应用于诸多领域[5]。就畜牧生产而言，种公畜的繁殖性能被育种学家高度关注，因此，探究miRNA对家畜精子发生的影响和作用机制成为了农业科学研究的热点之一[6]。

miR-34家族由3个成员组成，分别是miR-34a、miR-34b和miR-34c。miR-34b/c在精子发生过程中具有多种功能，既可以调节减数分裂过程，也可以调节精子发生后期的变形和获能等生物学过程。研究证明，多种动物的精子中含有miR-34b/c，这些成员在精子发生[7]和早期胚胎发育[8]中发挥着重要作用。miR-449家族也由3个成员组成，分别是miR-449a、miR-449b和miR-449c，它们被定位到Cdc20b基因的第2个内含子上，形成一个miRNA簇，在不同物种间高度保守。miR-449成员在控制细胞周期和表皮分化中起主要作用[9]。报道显示，miR-34/449家族在纤毛发生[10]、细胞周期调节[11]、细胞分化[12]、癌症[13]等多种生物学过程中发挥作用。虽然miR-34和miR-449最初被指定为单独的miR家族，但miR-34与miR-449家族成员具有相同的种子序列，且miR-34和miR-449被证明可以调节一组相同的靶基因，因此可以被视为一个miRNA簇，命名为miR-34/449[14]。鉴于这2个家族成员在基因表达调控中可能一起发挥作用，且2个基因家族之间具有相互代偿机制，本研究利用生物信息学方法探究miR-34b/c、miR-449a/b/c在牛精子发生中的作用机制及调控网络，通过预测共同调控miR-34b/c、miR-449a/b/c的转录因子(transcription factor,TF)以及受miR-34b/c、miR-449a/b/c共同调控的靶基因以构建miR-34b/c、miR-449a/b/c的分子调控网络，以期为后续miR-34b/c和miR-449a/b/c调控牛精子发生的相关研究提供理论基础。

1 材料与方法

1.1 材料

牛(Bostaurus,bta)、人(Homosapiens,hsa)、小鼠(Musmusculus,mmu)、大鼠(Rattusnorvegicus,rno)、鸡(Gallusgallus,gga)和猪(Susscrofa,ssc)等物种的miR-34b/c、miR-449a/b/c及其前体的序列信息来源于 miRBase(http:∥mirbase.org/index.shtml)在线数据库。

1.2 方法

1.2.1 miR-34b/c和miR-449a/b/c的相似性分析 利用MegaX[15]软件对miR-34b/c和miR-449a/b/c前体序列分别进行多序列比对，并计算各模型组合的BIC分数(Bayesian information criterion)，选取BIC分数最低的模型组合，即最大似然估计法(maximum likelihood)、Jukes-Cantor模型和Gamma分布分别构建miR-34b/c和miR-449a/b/c的系统进化树。

1.2.2 bta-miR34/449靶基因预测 通过在线软件TargetScan 7.2(http:∥www.targetscan.org/vert_72/)[16]、miRDB(http:∥www.mirdb.org/)[17]和miRWalk(http:∥mirwalk.umm.uni-heidelberg.de/)[18]预测bta-miR-34b/c、bta-miR-449a/b/c的靶基因，取三者预测结果的交集作为bta-miR34/449的靶基因。

1.2.3 靶基因功能预测 根据靶基因预测结果，选取bta-miR-34b/c和bta-miR-449a/b/c共同作用的靶基因。分别使用在线软件DAVID 5.8(https:∥david.ncifcrf.gov/)[19]和KOBAS(http:∥kobas.cbi.pku.edu.cn/)[20]对其进行GO功能和KEGG通路富集分析(物种为牛)，收集调整后P<0.05的条目并使用R语言程序包ggplot2绘制柱状图和气泡图。

1.2.4 转录因子预测 在NCBI数据库中检索出bta-miR-34b/c和bta-miR-449a/b/c的前体序列，使用AnimalTFDB(http:∥bioinfo.life.hust.edu.cn/AnimalTFDB/#!/)[21]对其上游2 kb至上游0.1 kb的序列进行转录因子结合位点预测，获取可能调控bta-miR-34/449的转录因子。

1.2.5 构建TF-miRNA-靶基因作用网络 选取同时调控bta-miR-34b/c和bta-miR-449a/b/c的转录因子，结合靶基因预测和String在线数据库(http:∥string-db.org/)[22]构建蛋白互作网络，使用Cytoscape 3.8.0在线软件(http:∥www.cytoscape.org/)[23]构建TF-miR-34/449-mRNA作用网络图。

2 结 果

2.1 miR-34/499的相似性分析

根据miRBase数据库的查询结果，牛miR-34/499分别位于第15和20号染色体上，对比牛、人、小鼠、大鼠、鸡和猪等物种的miR-34/449序列并构建系统进化树，结合序列比对发现，miR-34/449在进化中相对保守，成熟体序列基本一致(表1、2，图1)。

表1 各物种miR-34b/c的成熟序列

表2 各物种miR-449a/b/c的成熟序列

图1 ML法构建miR-34b/c(A)、miR-449a/b/c(B)系统进化树

2.2 靶基因预测

通过在线软件TargetScan、miRDB和miRWalk对bta-miR-34b/c、bta-miR-449a/b/c的靶基因进行预测，取三者的交集分别得到了115、151、161、149和152个靶基因(图2)。

A，bta-miR-34b；B，bta-miR-34c；C，bta-miR-449a；D，bta-miR-449b；E，bta-miR-449c

2.3 靶基因功能预测

bta-miR-34b/c和bta-miR-449a/b/c共同作用的靶基因分别为73和80个，用于功能分析的靶基因为49个(图3、表3)，如SYT1、SNX15、LGI1、ELMOD1、RCAN1等。对靶基因集进行GO功能和KEGG通路富集分析，结果显示，bta-miR-34b/c和bta-miR-449a/b/c的靶基因主要富集在精子发生、钙离子依赖性胞吐的调节、胰岛素分泌的调节等GO条目，以及HIF-1信号通路、甲状腺激素信号通路、Wnt信号通路等KEGG通路(图4)。

A，bta-miR-34b∩bta-miR-34c；B，bta-miR-449a∩bta-miR-449b∩bta-miR-449c；C，bta-miR-34b/c∩bta-miR-449a/b/c。图5同

表3 bta-miR-34b/c和bta-miR-449a/b/c的共同靶基因

图4 靶基因的GO功能(A)和KEGG通路(B)富集分析

2.4 转录因子预测

预测到bta-miR-34b和bta-miR-34c共同的转录因子有77个，bta-miR-449a、bta-miR-449b和bta-miR-449c共同的转录因子有17个，以上5个miRNAs共同的转录因子有2个(图5、表4)，分别为核受体亚家族2C组成员2(nuclear receptor subfamily 2C group member 2,NR2C2)和HICZBTB转录抑制因子1(HIC ZBTB transcriptional repressor 1,HIC1)。

图5 bta-miR-34b/c和bta-miR-449a/b/c的转录因子预测

表4 bta-miR-34b/c和bta-miR-449a/b/c的转录因子

2.5 TF-miRNA-mRNA作用网络的构建

根据bta-miR-34/449预测的转录因子、靶基因以及蛋白间互作关系构建TF-miRNA-mRNA作用网络，结果显示，bta-miR-34b/c、bta-miR-449a/b/c对精子发生的调控作用通过其与上游转录因子、下游靶基因共同构成的TF-miRNA-mRNA作用网络发挥作用，bta-miR-34b/c、bta-miR-449a/b/c受到上游转录因子NR2C2和HIC1的调控，同时会抑制AXL、E2F3、DAAM1等下游靶基因的表达，以实现对精子发生的生物学调控(图6)。

图6 TF-bta-miR-34b/c、bta-miR-449a/b/c-mRNA作用网络

3 讨 论

在基因表达过程中，miRNA能够通过转录后调控影响基因的表达，miRNA是基因表达的关键调节因子，有望成为生物标志物开发的候选分子。单个miRNA可以靶向数百个mRNA并影响许多基因的表达[24]。精子发生过程是雄性动物得以产生优质精子的重要生物学过程，miRNA已被证明参与精子发生的各个时期[25]。miR-34家族可以靶向许多细胞周期调节因子，如Notch1、CDK4和MYC，且TGIF2和Notch2在精子发生过程中是miR-34c的2个重要靶基因，Notch信号能够促进生殖细胞的分化[26]。有报道称，miR-34c在精母细胞和圆形精子细胞中高表达，抑制原代精母细胞中的miR-34c可以防止生殖细胞免受睾酮剥夺诱导的细胞凋亡，在体外培养的生殖细胞中miR-34c的过度表达会触发细胞凋亡[27-28]。miR-449家族成员也在精子发生中扮演着重要角色，miR-449a、miR-449b和miR-449c在减数分裂开始期间上调，并优先在精母细胞和精子细胞中表达[29]。大量研究已经揭示了miRNA-34和miR-449家族对精子发生各个时期的调控作用及机制，但其对精子发生过程中上游转录因子和下游靶基因构成的基因表达调控网络还有待进一步挖掘。

本研究通过查询miRBase在线数据库中牛、人、小鼠、大鼠等物种的miR-34b/c、miR-449a/b/c及其前体序列信息，运用相似性分析以及miRNA转录因子、靶基因预测、富集分析等技术手段对bta-miR-34b/c、bta-miR-449a/b/c的序列保守性、上下游调控机制及功能进行分析，得到了49个共同靶基因和2个转录因子，其中，靶基因和转录因子的筛选通过在线数据库TargetScan、miRDB、miRWalk及AnimalTFDB实现，以上数据库均包含牛的相关信息或支持通过序列直接预测，尽管miRNA在物种间具有较强的序列保守性，但是使用其他物种miRNA的相关信息预测牛miRNA的靶基因或转录因子仍然不是最优选择，因此本试验通过对预测工具的合理选择在一定程度上提高了研究结果的准确性。

本研究预测到的2个转录因子(NR2C2和HIC1)在精子发生中发挥着重要作用。NR2C2也被称为睾丸核受体4(testicular orphan receptor 4，TR4)，是核受体家族的成员[30]。在小鼠睾丸中，NR2C2主要表达于初级精母细胞和圆形精子细胞中，并在精子发生的减数分裂前期和连续减数分裂前期起重要作用。NR2C2基因敲除小鼠的生殖潜能显著降低，各个阶段的精子数量明显降低，初级精母细胞降解，生精小管坏死扭曲，精子发生减数分裂前期和连续减数分裂期延长和中断，导致细胞分裂中期延长和生精细胞异常[31]。HIC1是一种锌指转录抑制因子，可通过与沉默交配型信息调节2同源物1(sirtuin 1，SIRT1)启动子结合，控制细胞生长和死亡，以响应p53依赖性凋亡DNA损伤[32]。Jenal等[33]研究发现，HIC1是细胞周期和凋亡调节因子E2F转录因子1(E2F transcription factor 1,E2F1)的新转录靶点，E2F1通过包含HIC1 P0启动子上的TATA盒在内的2个E2F DNA结合位点诱导HIC1的表达。HIC1与包含HIC1结合共有位点的E2F1启动子基因结合，通过调控E2F1转录抑制细胞生长。Uchida等[34]研究证明，HIC1缺失促进青春期前睾丸平滑肌细胞的增殖，导致成年睾丸中生精小管周围的PMC数量增加，并最终影响精子发生。

本研究筛选到的部分靶基因已有报道显示参与调控了精子发生过程。受体酪氨酸激酶AXL是TAM家族的成员，在细胞增殖、分化、凋亡和新陈代谢等生物学过程中发挥着重要作用[35]。AXL是牛睾丸支持细胞中miR-34c的直接靶基因，miRNA-34c通过靶向AXL基因调节牛支持细胞增殖和凋亡[36]。E2F3是E2F家族中的一员，是唯一被发现在精原细胞、细线期精母细胞和支持细胞中表达的E2F成员。在支持细胞特异性视网膜母细胞瘤蛋白(RB)敲除小鼠睾丸内注射shRNA-E2F3，结果显示，支持细胞中因RB缺失引起的生精功能障碍通过E2F3的体内沉默得到部分恢复[37]，暗示RB/E2F3途径在精子发生过程中发挥重要调节作用。形态发生相关激活因子1(dishevelled associated activator of morphogenesis 1,DAAM1)是一种formin家族蛋白，通过Wnt信号通路参与非支链肌动蛋白丝的成核、形态发生和器官发生等重要生物学过程[38]。Pariante等[39]评估了DAAM1与大鼠睾丸形态发生的可能关联，结果显示，在有丝分裂阶段，DAAM1与支持细胞、生殖细胞和精原细胞中的肌动蛋白共享其定位，在减数分裂期间，2种蛋白都存在于精母细胞中，而在形成的血睾丸屏障中只能检测到肌动蛋白，DAAM1在整个精子发生过程中跟随顶体系统的发展，最终保留在成熟精子的细胞质内。乳酸脱氢酶A(lactate dehydrogenase A,LDHA)在精子能量代谢过程中起重要作用，热应激能够增加猪睾丸中葡萄糖转运蛋白3(glucose transporter type 3,GLUT3)和LDHA的表达，继而增加支持细胞中乳酸产生[40]。小鼠中乳酸脱氢酶C(lactate dehydrogenase C,LDHC)基因的缺失导致精子无法获能并在精卵结合中不能穿过透明带[41]，而LDHA作为转基因引入LDHC缺失小鼠中，转基因小鼠获得了正常的精子功能和生育能力[42]。Klinefelter综合征(Klinefelter syndrome,KS)是男性性腺机能减退和性染色体异常的遗传疾病，患有KS的男性睾丸中LDHA mRNA的表达量减少与睾丸内乳酸含量的下降相关[40]。以上证据不同程度地印证了本研究预测的靶基因在精子发生中扮演着重要角色。此外，本研究预测得到的还没有被广泛报道的其他基因也具有进一步研究的价值。

靶基因GO功能富集分析表明，bta-miR-34b/c和bta-miR-449a/b/c的靶基因主要富集在精子发生、钙离子依赖性胞吐的调节、胰岛素分泌的调节等GO条目，推测bta-miR-34b/c和bta-miR-449a/b/c的靶基因在精子发生和精子获能中起作用。相关报道显示，miR-34c参与小鼠睾丸精子发生和雄性生殖细胞(male germ-line stem cells,mGSCs)凋亡的调节[7]。miR-34b/c和miR-449a/b/c对于正常精子发生和男性生育不可或缺[43]。miR-34b和miR-34c存在于与透明带结合的精子中，表明其可能与精子受精过程有关[8]。KEGG通路富集分析，bta-miR-34b/c和bta-miR-449a/b/c的靶基因显著富集在HIF-1信号传导途径、甲状腺激素信号通路和Wnt信号通路等。缺氧诱导因子(hypoxia-inducible factor,HIF)是细胞对缺氧反应的关键转录调节因子，HIF由α亚单位(HIF-1α、HIF-2α或HIF-3α)和β亚单位(HIF-β)组成，HIF-1α和HIF-2α密切相关，它们与HIF-β亚基的相互作用分别产生HIF-1和HIF-2，HIF-1和HIF-2 2种活性转录因子均与缺氧反应元件(hypoxia-response element,HRE)结合，刺激与血管生成、能量代谢、细胞生长和细胞周期进展等相关的靶基因调节生物学过程[44]。HIF-1α和HIF-2α在睾丸支持细胞中的表达已经得到证实[45]，且HIF-2参与促卵泡素(follicle-stimulating hormone,FSH)对支持细胞增殖的调节[46]。Gruber等[47]对出生2 d的小鼠进行HIF-2α消融手术，导致小鼠支持细胞无法建立完整的血睾屏障，并伴有睾丸大小和重量减少，且经消融手术处理的小鼠最终不育。本研究中靶基因富集到的另一个重要通路是Wnt信号通路，Wnt信号是一种高度保守的细胞间通信机制。典型的Wnt信号传导也称为Wnt/β-连环蛋白通路，通常被认为是一种干细胞自我更新机制[48]。原始生殖细胞的特异性分化和雄性胎儿生殖道的正常发育依赖于Wnt/β-catenin信号传导[49]。Takase等[50]研究报道，睾丸中Wnt信号通路特别有助于精原干细胞(spermatogonial stem cells,SSC)的增殖，证实了本研究结果的可靠性。

4 结 论

bta-miR-34b/c、bta-miR-449a/b/c通过上游转录因子和下游靶基因组成的分子调控网络共同调节精子发生过程。