绵羊睾丸差异基因及蛋白质互作网络关系研究

2022-02-15刘在霞孙燕勇付绍印何小龙张文广刘永斌

中国畜牧兽医 2022年1期

刘在霞，段仕，孙燕勇，吕琦，付绍印，何小龙，张文广，刘永斌

(1.内蒙古农业大学动物科学学院，内蒙古自治区动物遗传育种与繁殖重点实验室，呼和浩特 010018；2.内蒙古自治区农业基因组大数据工程研究中心，呼和浩特 010018；3.内蒙古自治区农牧业科学院，呼和浩特 010031；4.和田师范专科学校，和田 848000)

绵羊在全球畜牧业中占有很大比例，而繁殖效率是畜牧业收入水平的关键指标，大多数绵羊品种是单胎，季节性繁殖，繁殖率低[1]，迫切需要采用两年三产的繁殖模式。在提高绵羊繁殖效率的实践中，对母羊的研究较多，最常见的方法是利用生殖辅助技术来改善发情时间[2]、增加母羊的排卵次数[3]和进行人工授精[4]。但绵羊的成功繁殖取决于绵羊两性，获能、受精和胚胎发生过程中的失败可能是公羊精子的原因[5]。研究发现，其他物种精子中有超过14 000个转录本在受精时被转移到卵子中，这些转录影响早期胚胎发育和受孕是否成功，Hodge等[6]收集了3个绵羊品种的精液，在不同品种和射精质量对比中共鉴定出754个差异表达基因(DEGs)，这些基因在维持精子发生、受精、受孕、胚胎发育和后代生产性能中发挥着重要作用。睾丸是产生精子和分泌雄激素的重要器官，睾丸发育和精子发生是一个高度复杂和精确的过程，在发育过程中睾丸的内部结构和功能发生了一系列的变化[7-8]。随着睾丸大小的改变，产生精子的效率也会改变。营养不足影响性成熟雄性睾丸质量，降低生精效率和精子运动速度，并增加性成熟雄性绵羊的精子DNA损伤[9]。此外，睾丸发育和精子发生主要受mRNA调控，这是动态并具有阶段性的[10]。蛋白质是一切生物的物质基础，蛋白质-蛋白质互作(protein-protein interaction，PPI)在任何生长、繁殖和新陈代谢中都起着极其重要的作用，即使在单个细胞也是如此[11]，PPI通过连接相关的蛋白质，在每个单细胞中启动各种功能或结构蛋白的作用[12]。Jansen等[13]使用贝叶斯网络方法预测全基因组酵母中的PPI，用串联亲和纯化试验(TAP)验证了预测，提供了酵母PPI的全面视角。为全面探索绵羊睾丸在不同品种、不同营养水平和不同发育阶段差异基因及PPI网络，本研究使用课题组前期利用机器学习算法预测的绵羊全基因组820 072对PPI关系[14]，构建睾丸的PPI网络，识别可能调控绵羊睾丸发育和精子发生的关键基因，以期进一步提高绵羊的繁殖效率，为今后的绵羊繁殖研究提供理论依据。

1 材料与方法

1.1 绵羊睾丸转录组数据获取

绵羊睾丸转录组数据来源于NCBI的高通量测序数据存储库SRA(sequence read archive)，获取码：PRJEB19199、PRJNA282344和PRJNA421437，共获得74个转录组数据，包括50只特克塞尔羊与苏格兰黑头羊杂交后代(TB)、16只澳洲美利奴绵羊(8只在65 d内增重10%的高饲粮绵羊(M+)和8只在65 d内减重10%的低饲粮绵羊(M-))、9只不同发育阶段的小尾寒羊(ST，2、6和12月龄各3只)。

1.2 睾丸中基因表达量的计算

利用SRAtoolkit软件的fast-dump参数将高通量测序数据转变为fastq文件，利用FastQC软件对原始测序数据进行质量评估，使用Trimmomatic软件过滤接头、低质量和未知碱基数目过多的reads，获得高质量的clean data；利用Hisat2软件将clean data与绵羊参考基因组(Oar_v4.0)比对，利用Stringtie软件对转录本进行组装，之后采用R软件的Ballgown包计算基因在每个样本中的转录本表达量；使用FPKM(Fragment Per Kilobase of per Million mapped reads)值来衡量基因的表达水平。

1.3 PPI网络构建及网络中关键基因的选择与分析

PPI数据集来自于课题组前期采用随机森林(random forest,RF)和十折交叉验证法构建的预测模型，对28 592个绵羊蛋白质最终预测到820 072对蛋白质具有互作关系[14]。将820 072对互作关系中的睾丸表达基因和睾丸DEGs使用Cytoscape软件(version 3.7.2)构建睾丸PPI网络。使用Cytoscape软件中的MCODE(molecular complex detection)插件[15]对PPI网络中的节点进行密度聚类，计算网络中最重要的节点和集合(cluster)，PPI网络中选择集合标准如下：Degree Cutoff≥2;Node Score Cutoff=0.2;K-core≥2;Max Depth=100。

1.4 睾丸DEGs分析

绵羊睾丸DEGs数据集来自3类：①不同绵羊品种(8只M+和8只TB绵羊共16个RNA-Seq)，使用R语言limma软件包对2个品种(M+vsTB)进行DEGs分析，在差异基因检测过程中，将log2|FoldChange|≥1且P<0.05作为筛选标准；②通过文献挖掘16只不同饲喂营养水平的澳洲美利奴绵羊(M+vsM-)[16]；③9只不同发育阶段(2、6、12月龄)ST，分别代表睾丸发育过程中性成熟前、性成熟、性成熟后3个阶段[17]。

1.5 构建加权绵羊睾丸差异基因共表达网络

构建加权基因共表达网络(weight gene co-expression network analysis,WGCNA)[18]，为减少因测序批次不同带来的误差，选取同一批次的样本进行WGCNA分析，共表达网络的2 058个DEGs表达量来自8个M+样本。采用分层聚类和动态树切割算法对模块进行识别，最小模块数设置为30个基因，并将合并高度相似的模块的最小高度设置为0.25。

1.6 功能富集分析

为了研究模块及其基因潜在功能，使用在线生物信息学数据库Metascape (https:∥metascape.org/)对所选模块中的基因进行富集分析，并使其服从GO功能和KEGG通路富集分析，P<0.01作为筛选标准(Min Overlap=3，Min Enrichment=1.5)，利用Cytoscape插件AutoAnnotate自动识别基因集簇,并对基因集进行注释。

2 结果

2.1 绵羊睾丸转录组数据分析

2.1.1 睾丸表达基因分析通过对74个样本进行转录组分析，共鉴定到18 650个睾丸表达基因。为降低基因表达的假阳性，筛选3%的样本中FPKM≠0，且至少有一个样本FPKM≥3，最终获得11 884个基因用于后续分析，其中425个基因位于X染色体，11 150个基因位于常染色体，309个为“unplaced”。

2.1.2 睾丸PPI网络构建对睾丸中最终获得的11 884个基因在820 072对互作关系中提取PPI，仅有8 291个基因获得12 621个蛋白质(节点，Node)，产生237 366对PPI(边，Edge)，占绵羊总PPI的28.9%，并构建237 366对的PPI网络(图1A)。对PPI网络进行MCODE分析，共确定366个集合，其中集合1的分值最大为142.69，且关系对最多，是由200个Nodes和18 535个Edges构成；集合7中Nodes最多为411个，对排名前100的集合绘制柱状图(图1B)。

2.2 绵羊睾丸DEGs分析

2.2.1 绵羊睾丸DEGs的PPI网络构建不同绵羊品种(M+vsTB)共获得735个DEGs，其中上调基因324个，下调基因411个；不同营养饲喂水平(M+vsM-)共获得2 243个DEGs；不同发育阶段共获得13个DEGs。综上，去除重复的基因后共获得2 058个绵羊睾丸DEGs。对2 058个DEGs在237 366对PPI中提取互作关系，构建差异基因PPI网络，发现仅有1 644个DGEs对应2 020个蛋白质，产生46 169对PPI，对PPI网络进行MCODE分析，确定了87个集合(图2)。选择前4个得分高的集合，集合1得分最高为25.203，有13个基因；集合2得分为8.759，有26个基因；集合3得分为8.605，有45个基因；集合4得分为8.125，有12个基因(图3)。

A,睾丸DEGs的PPI网络；B,PPI蛋白集合边、点统计A,PPI network of testis DEGs;B,Nodes and Edges counts of PPI图2 睾丸DEGs的PPI网络分析Fig.2 PPI network analysis of testis DEGs

A～D，集合1～4代表PPI网络的核心基因A-D，Cluster 1-4 represented the core gene of PPI network图3 PPI网络MCODE分析Fig.3 PPI network MCODE analysis

2.2.2 构建加权绵羊睾丸DEGs共表达网络对2 058个DEGs进行WGCNA分析，以确定数据集中的共表达基因，参数设置为：beta=18，R2>0.8。根据基因表达量进行相关度的聚类，聚类度较高的基因被分配到一个模块中，不同颜色代表不同的模块，而灰色模块内基因代表其相关性较低，无法被分配到其他模块，不进行后续分析，本研究共获得7个共表达模块(图4A)，模块中的基因数目为51～929个(图4B)。

A，基因平均连锁层次聚类树状图；B，模块中的基因数量，模块由直方图颜色指定A,Average linkage hierarchical clustering dendogram of the genes；B,The number of genes in the module，modules were designated by color histogram图4 加权基因共表达网络分析Fig.4 Weighted gene coexpression network analysis

2.2.3 DEGs功能富集分析利用Metascape数据库探索模块的功能和途径，Blue模块共富集到106个功能术语，主要参与雄性生殖系统相关的功能，包括雄性配子产生(gamete generation)、繁殖(fertilization)、精子发生(spermatogenesis)、鞭毛运动(cilium movement)、细胞周期的过渡(cell cycle transition)、AMPK信号通路(AMPK signaling pathway)等。使用Cytoscape插件AutoAnnotated可视化网络，把相似度>0.3的功能术语用边连接，共产生14个功能集合术语，富集网络图显示了功能集合内和集合间功能富集结果的相似性，网络中的节点代表一个功能术语，边代表功能之间的连接关系，同一颜色节点代表一个功能集合富集的过程，发现它们高度连接并聚成一个完整的网络(图5)，其中配子产生(gamete generation)包括：生殖细胞发育(germ cell development)、精子发生(spermatogenesis)、精子细胞分化(spermatid differentiation)、精子发育(spermatid development)、雄性配子产生(male gamete generation)和多细胞生物体内参与生殖的细胞过程(cellular process involved in reproduction in multicellular organism)。表明Blue模块内基因是参与精子发生和精子细胞发育过程的重要的基因集合。

图5 睾丸差异基因富集结果全局网络图Fig.5 Global network diagram of testicular DEGs enrichment result

2.2.4 关键基因的筛选对2 058个DEGs通过PPI网络的MCODE和WGCNA联合分析，考虑到精子发生和精子细胞的发育过程相关基因主要分布在Blue模块，因此利用Blue模块内基因和PPI网络关键的4个集合，确定了25个基因(表1、图6)，发现仅来自不同营养水平的基因有19个：基质重塑相关8(matrix remodeling associated 8,MXRA8)、富含亮氨酸重复序列和免疫球蛋白样结构域3(leucine rich repeats and immunoglobulin like domains 3,LRIG3)、神经细胞黏附分子2(neural cell adhesion molecule 2,NCAM2)、丝氨酸/苏氨酸激酶17b(serine/threonine kinase 17b,STK17B)、V-集和跨膜结构域(V-set and transmembrane domain containing,VSTM4)、恩贝酸(embigin,EMB)蛋白、免疫球蛋白超家族成员3(immunoglobulin superfamily member 3,IGSF3)、CDC样激酶(CDC like kinase 4,CLK4)、VRK丝氨酸/苏氨酸激酶2(VRK serine/threonine kinase 2,VRK2)、淀粉样前体蛋白(amyloid beta precursor protein,APP)、磷脂酰肌醇-4-磷酸3-激酶催化亚单位2β(phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 beta,PIK3C2B)、尼沙林(nischarin,NISCH)、分拣微管连接蛋白2(sorting nexin 2,SNX2)、红细胞膜带4.2蛋白(erythrocyte membrane protein band 4.1 like 2,EPB41L2)、SHC适配器蛋白(SHC adaptor protein 4,SHC4)、小核糖核蛋白多肽B2(small nuclear ribonucleoprotein polypeptide B2,SNRPB2)、SRA茎环相互作用的RNA结合蛋白(SRA stem-loop interacting RNA binding protein,SLIRP)、RNA结合基序单链相互作用蛋白2(RNA binding motif single stranded interacting protein 2,RBMS2)和RNA结合基序蛋白46(RNA binding motif protein 46,RBM46)；仅来自不同品种的基因有4个：生长因子受体结合蛋白14(growth factor receptor bound protein 14,GRB14)、酪氨酸激酶2(tyrosine kinase 2,TYK2)、异种核糖核蛋白A1样 (heterogeneous nuclear ribonucleoprotein A1-like,LOC101109899)和ELAV类RNA结合蛋白1(ELAV like RNA binding protein 1,ELAVL1)；来自不同营养水平和品种的基因有2个：受体酪氨酸激酶3(receptor tyrosine kinases3,TYRO3)和肝细胞黏附因子(hepatic and glial cell adhesion molecule,HEPACAM)。25个基因分布于13条染色体上，且各基因间存在互作关系，如SLIRP与LOC101109899基因、RBM46、RBMS2和ELAVL1基因存在互作关系。

表1 Blue模块25个关键基因

最外圈表示绵羊全基因组27条染色体的编号与长度，格越长表示染色体越长，同时外圈展示了25个关键基因在染色体上的位置；圈里的热图代表关键基因的表达量，表达量越高，颜色越深；最内圈的连线表示关键基因间互作关系The outermost ring indicates the number and length of 27 chromosomes in the whole genome of sheep，the longer the grid is,the longer the chromosome length is,meanwhile,the outermost ring shows the location of 25 key genes on the chromosome；The heat map in the circle represents the expression of key genes，the higher the expression,the darker the color；The innermost ring lines represent interactions between key genes图6 关键基因Circos图Fig.6 The key genes Circos map

SLIRP基因表达量最高为100.3，RBM46基因表达量为42.47，ELAVL1、RBM46和SLIRP基因显著富集到AMPK信号通路和芳香族化合物分解代谢过程；TYK2基因表达量为39.18，主要参与白细胞介素-12、-23和-27介导的信号通路；TYRO3基因表达量为21.54，主要参与配子产生、精子发生、蛋白磷酸化。

3 讨论

PPI在任何生命的生长、繁殖和新陈代谢中都起着极其重要的作用，睾丸发育和精子发生是影响雄性绵羊生育力的重要因素。虽然已经有大量关于睾丸发育和精子发生的研究，但对睾丸PPI网络的研究较少，本研究对睾丸表达的11 884个基因和2 058个DEGs分别构建绵羊PPI网络，发现睾丸的PPI网络可较为紧密地形成一个大的网络，只有少部分关系稀疏，对2 058个DEGs提取互作关系，仅有1 644个DGEs对应2 020个蛋白质，产生46 169对PPI，通过MCODE、WGCNA和Metascape分析，发现Blue模块是与精子发生相关的基因集。 Silva等[19]构建了人类睾丸PPT网络，共鉴定出1 778个蛋白质及它们之间的32 187对相互作用，并有少量孤立的簇；富集分析显示，最重要的生物学过程类别都与生殖相关，如有性生殖(P=3.71E-156)、精子发生(P=5.36E-147)和受精(P=4.33E-44)，与本研究结果相近。

睾丸是转录最复杂的器官，其表达的基因数量是所有哺乳动物器官中最多的，据报道，睾丸在人和其他物种中广泛的转录超过所有蛋白质编码基因的80%，长链非编码RNA(lncRNA)基因和基因组的基因间序列，包括假基因和转座因子等，向睾丸的倾斜表达更为明显，对人类和小鼠精子发生的单细胞转录组研究发现，广泛的转录是通过转录扫描机制来纠正DNA损伤，从而维持雄性生殖系中的DNA序列完整性。在精子发生过程中表达的基因在转录链上表现出较低的突变率，且在群体中具有较低的多样性[20]。Soumillon等[21]对人、恒河猴、小鼠、负鼠和鸡的不同组织器官进行转录组测序，发现睾丸分别表达20 322、20 078、21 819、18 198、14 574个蛋白质编码基因，且常染色体蛋白质编码基因在睾丸中的转录频率高于其他器官。Clark等[22]对3头成年公羊和3头成年母羊的主要组织和5类细胞进行转录组测序，发现有19 921个基因至少在一个组织中表达。本研究通过严格的表达量筛选最终获得11 884个睾丸表达基因，由于与参考基因组比对、过滤或表达量筛选时参数和条件与前人不同而产生差异，由此可见睾丸表达基因较为丰富。

哺乳动物睾丸是先天免疫的特殊器官，分泌多种免疫调节因子，对外来和自身抗原耐受[23]，Toll样受体(Toll-like receptors，TLR)在支持细胞和生殖细胞中表达，在被配体激活后引发睾丸的先天免疫，Tyro3基因负调控支持细胞中的TLR3信号转导[24]，降低炎症因子表达水平，识别与结合凋亡细胞，控制睾丸对病原体的先天免疫应答[25]。本研究通过MCODE分析发现，Tyro3基因位于分值最高的集合1，表明其与周边基因的密集程度很高，处于关键地位，且其在不同营养水平和不同品种组都是差异基因，表明其在精子发生中可能起关键作用。通过功能富集分析发现,Tyro3基因主要参与配子产生、精子发生、蛋白磷酸化等生物学过程。唐红梅[26]研究发现，Tyro3基因在哺乳动物的神经、免疫、生殖和血管等组织中广泛表达，并在神经发育、免疫调节、精子发生、NK细胞分化、血小板功能及血管发生中发挥着重要的作用。

ELAVL1基因位于1号染色体上，又名HuR或Hua，与mRNA的3′-非翻译区结合，调节转录的稳定性。Ghosh等[27]研究发现，小鼠出生后ELAVL1基因的缺失可导致造血器官萎缩、肠绒毛丢失、阻塞性小肠炎，并在10 d内死亡。Chi等[28]对小鼠敲除和过表达HuR基因发现，HuR基因在原始生殖细胞中的失活与正常的减数分裂后细胞形成和精子细胞成熟是矛盾的，此外，在圆形精子细胞中过表达HuR基因延缓了精子细胞的分化，敲除ELAVL1基因使生殖细胞分化缺陷、精子完全失活，导致雄性不育[29]。本研究中，ELAVL1基因处于集合4，主要参与AMPK信号通路。研究发现，在支持细胞中AMPK信号通路调节能量代谢，连接复合体的稳定性和增殖，一旦平衡被打破，睾丸的微环境和精子的质量都会受到影响，如在α1 AMPK整体敲除小鼠中，精子表现为头部异常、鞘弯曲和活动能力受损[30]。以上均表明ELAVL1基因与精子发生有关。

RBM46基因位于7号染色体上，是一种新的生殖细胞特异性因子，可调节减数分裂和细胞凋亡，在性腺组织的生殖细胞中特异表达。研究发现，RBM46基因突变的斑马鱼中，RBM46基因的特异性破坏导致突变体的雄性化和不育，虽然突变体的精原细胞基本正常，但精子发生受损，减数分裂未发生，表明RBM46基因在斑马鱼生殖细胞第一次减数分裂过程中起着至关重要的作用[31]。Zhai等[32]研究发现，RBM46基因分布在小鼠胚胎干细胞的细胞质中，通过靶向β-Catenin mRNA降解调节胚胎干细胞分化，β-Catenin是Wnt信号通路中的一个关键效应子，RBM46基因转录后调控在胚胎干细胞分化中起重要作用，RBM46基因过表达导致胚胎干细胞向滋养外胚层分化，而敲除RBM46基因导致胚胎干细胞向内胚层分化。本研究中RBM46基因表达量为42.47，处于集合4，显著富集到AMPK 信号通路和芳香族化合物分解代谢过程，并与SLIRP、LOC101109899、RBMS2和ELAVL1基因存在互作关系，但其调控分子机制尚不清楚。SLIRP基因是核受体(NRs)抑制因子，其存在于核内与NR靶基因相关的复合物中，在骨骼肌、心脏、肝脏和睾丸中表达量最高[33]，本研究中SLIRP基因表达量为100.3，与此结论相一致。 Colley等[34]研究发现，当敲除SLIRP基因的纯合雄性小鼠与野生型雌性杂交时，平均产仔数比野生型雄性和雌性小鼠杂交产仔数减少了1/3，且敲除小鼠的精子运动能力明显低于野生小鼠，电镜观察发现敲除小鼠精子中段/环连接处断裂，线粒体形态发生改变，表明SLIRP基因的缺失导致与精子结构和线粒体形态受损相关的弱精子症。Shan等[35]研究表明，SLIRP基因调节雄性繁殖能力且与精子运动呈正相关，与正常精子相比，弱精子症患者精子中SLIRP基因的mRNA和蛋白质表达量降低，氧化损伤增加，能量代谢降低。

TYK2基因位于5号染色体上，是JAK激酶蛋白家族中第一个鉴定的成员，JAK家族还有3个成员：JAK1、JAK2和JAK3[36]。TYK2基因在细胞因子转导过程中起着关键作用，广泛参与免疫、细胞增殖、生长、分化、迁移和细胞凋亡等过程[37]，本研究中TYK2基因显著富集到参与白细胞介素-12、-23和-27介导的信号通路，与前人研究结果相似。D’Cruz等[38]通过免疫印迹和激光扫描共聚焦显微镜检测到JAK家族的4个成员中仅TYK2基因在人精子尾部中表达，表明其可能调节精子活力。本试验中TYK2基因表达量为39.18，处于集合2，与RBM46、SLIRP等基因存在互作关系，但未对JAK1、JAK2和JAK3基因在绵羊精子尾部的表达情况进行验证，仍有待进一步研究。