RNA-seq转录组测序分析不同季节对槟榔江水牛精液品质的影响

耿红红，张敬旸，李　莲，王根林

(南京农业大学 动物科技学院，南京 210095)



RNA-seq转录组测序分析不同季节对槟榔江水牛精液品质的影响

耿红红，张敬旸，李莲，王根林*

(南京农业大学 动物科技学院，南京 210095)

摘要：本研究旨在利用RNA-seq技术对不同季节槟榔江水牛全血mRNA进行分析，找出不同季节精液品质候选基因。选取5头生长状况相近的3～5岁槟榔江水牛，进行春、夏、秋、冬4个季节的精液品质(鲜精活力、精子畸形率、射精量、精子密度)检测。根据精液品质，选取差异极显著的春(BL)、夏(SBL)两个季节的水牛样品进行RNA-seq转录组测序分析。颈静脉采集全血，抽提mRNA并组内混合制样，利用RNA-seq技术分析不同季节槟榔江水牛全血mRNA差异表达。根据差异倍数>1.5，P<0.05的标准进行筛选，得到差异基因，再利用DAVID 6.7软件对差异基因进行GO和KEGG分析，通过RT-PCR验证基因芯片结果的准确性。结果表明：差异表达基因84个，其中上调基因56个，下调基因28个；GO分析发现，这些差异基因与应激反应、免疫反应有关；KEGG分析表明，差异基因显著富集到肿瘤坏死因子(TNF)和T细胞受体信号通路。其中，神经营养因子受体类似物死亡结构域(NRADD，又叫热休克蛋白因子1 (HSF1))、原癌基因FOS和JUN，与精液品质相关性强，可作为进一步研究的靶基因。

关键词：RNA-seq；槟榔江水牛；差异表达基因；精液品质

槟榔江水牛(Betelnut-jiang buffalo)是中国发现唯一的河流型水牛，主产于腾冲槟榔江上游，槟榔江水牛体型较小，饲料消耗较低，而且产奶量相对较高，牛奶营养丰富，是比较好的乳用品种。当前，人工授精技术及冷冻精液的日益普及，具有优良性状的种公牛越来越受到重视。槟榔江水牛作为改良产奶量低的沼泽型水牛的优良父本，探讨其遗传优势有现实意义。

种公牛的精液品质受到季节、品种、年龄、营养和饲养水平等要素的影响，其中季节变化，尤其是温度湿度的变化对精液品质的影响最为显著。研究表明，牲畜自身的生理状态和精子发生需要一个适宜的外界环境作为保障[1]。季节变化中夏季温度较高时会极大的影响公牛的精液品质和生产效益[2-3]。衡量精液品质的指标主要有射精量、精子密度、原精活力、冻后精子活力、畸形率和有效精子数等。相关数据显示，种公牛夏秋季采精量明显低于冬春季，其中采精量夏季<秋季<冬季<春季[4]。西门塔尔种公牛原精活力和精子密度随季节变化明显，具体为春季至秋季呈下降趋势，尤其在6～9月份，达到最低值，冬季后又明显上升[5]。M.Sadegh等[6]通过对印度水牛夏季和冬季繁殖性能的比较发现，冬季精液采集量、密度、精子活力和直线运动精子数均高于夏季，而精子头部畸形率则夏季高于冬季。由此可见，季节变化会引起精液品质表观参数发生规律性的变化。

随着快速发展的分子生物学技术，加之候选基因法在动物遗传育种中的广泛应用，人们不仅可以在DNA、RNA及蛋白表达水平上研究影响精液品质的候选基因或与其紧密连锁的分子标记，还可以将分子标记与多基因效应合并使用进行标记辅助选择，从动物育种角度进行提前选择，最终借助遗传的改变改善精液品质，以提高经济效益。转录组测序技术(RNA-seq)主要是用高通量测序技术将mRNA、sRNA 以及ncRNA 等的序列测定出来，以期能够全面快速地获得几乎所有的所需研究样本的转录本，并反映出它们的表达水平[7]。RNA-seq在真核生物的基因表达调控、医学和动植物遗传育种等方面的研究具有不可估量的潜力[8-9]。季节变化对公牛精液品质的研究大多数停留在表观指标的规律性变化，而对于分子层面甚少涉及，这对选育工作极为不利，利用转录组学的技术研究精液品质亦鲜见报道。故本研究对5头成年中国槟榔江水牛公牛在春、夏、秋、冬4个季节进行精液品质的跟踪检测，并采用RNA-seq技术进行差异表达基因筛选分析，探究季节对槟榔江公牛精液品质的影响，为种公牛的选育提供理论依据。

1材料与方法

1.1样本收集

自湖北省武汉市兴牧科技有限公司种公牛繁育场选取具有相同的饲养条件，年龄3～5岁，体型相近；健康，无患病史的成年槟榔江水牛种公牛5头。

1.2精液品质相关指标的测定

根据国家牛冷冻精液检测标准(GB/4143-2008)的相关规定，对精液指标包括射精量、鲜精活力、密度、畸形率进行测定。具体方法：采用台牛法采集精液，射精量即采即测，一周两次；精子密度采用密度仪测定，直接读取数据；精子活力通过电视显微装置观察，每样片观察3个视野，并观察不同液层，进行评估；畸形率的测定采用姬姆萨染色法，并用血球分类计数器进行计数，主要步骤包括制片、固定、染色、镜检，镜检时每个抹片至少观察200个精子。

1.3RNA的提取、建库及测序

采用Trizol法进行槟榔江水牛外周血总RNA的抽提；提取所得RNA通过紫外分光光度计分析、琼脂糖凝胶电泳和Agilent 2200 Bioanalyzer进行质检 (A260 nm/A280nm≥1.5，A260 nm/A230 nm≥1，RNA完整性≥7)，若电泳发现有基因组DNA污染，则用DNA酶消化处理；送广州锐博生物技术公司完成RNA测序。

1.4芯片结果数据分析与处理

根据 GeneChip Scanner3000 扫描芯片结果，用Command Console Software 3.1 读取得到的原始数据，采用GeneSpring Software 11.0归一化处理质控合格的数据(MAS 5.0)进行数据分析。本研究筛选的标准为差异倍数值(Fold change值，即FC值)>1.5且P<0.05，利用网络数据库及相关生物信息学软件对差异表达基因进行Pathway分析、GO分析以及层次聚类分析等，对数据中所涉及的分子功能、分子代谢通路以及可能的表达调控机制进行分析。

1.5差异表达基因的qRT-PCR验证分析

1.5.1反转录生成cDNA 使用TaKaRa PrimeScript®RT reagent kit with gDNA eraser 反转录试剂盒将总RNA反转录合成cDNA，20 μL PCR反应体系：mRNA模板4 μL，Premix Taq Version 2.0 4 μL，ddH2O 12 μL。PCR反应条件：37 ℃ 15 min，85 ℃ 5 s，4 ℃保存。

1.5.2荧光定量PCR 采用SYBR Green方法进行实时定量PCR反应。所需引物用Primer 5软件设计，由南京金斯瑞生物科技有限公司合成，引物信息详见表1，20 μL反应体系：cDNA模板2 μL，ROX(SYBR Green Master)10 μL，ddH2O 7.2 μL，上下游引物各0.4 μL。在7000 Sequence Detection System(Applied Biosystems)仪上进行，每个样本均重复3次。使用仪器自带软件检测各样品CT值，并取平均值。mRNA表达水平用2-ΔΔCT方法来表示。

表1　RT-PCR引物序列

2结果

2.1不同季节槟榔江水牛精液指标差异分析

槟榔江水牛不同季节精液指标差异的分析结果见表2。由表2可知，槟榔江水牛采精量在季节之间差异不显著(P>0.05)。精子活力春冬和夏秋季节差异均不显著，春冬与夏秋季节差异极显著，在4个季节的变化为春季>冬季>秋季>夏季，说明夏季的到来会对槟榔江水牛精子活力产生影响，这种影响一直会持续到秋季，直到冬季原精活力才会恢复到之前的水平。槟榔江水牛精子密度春冬季节差异显著，并与其他季节差异极显著，且春季>冬季>秋季>夏季，表明夏季槟榔江水牛精子密度明显降低，夏季结束后，槟榔江水牛精子密度逐渐恢复到原有水平。夏季槟榔江水牛冻后精子畸形率夏季与其他季节相比差异极显著，春秋和冬季则差异不显著，且4个季节冬季<春季<秋季<夏季。

通过对槟榔江公牛不同季节精液相关指标的分析发现夏季是槟榔江公牛精液品质最差的季节，秋季次之，冬春季节差异不显著，选择具有代表性的夏季和春季作为试验组和对照组进行转录组测序。

表2　槟榔江水牛不同季节精液指标比较(n=5)

同行肩标不同小写字母表示不同季节差异显著(P<0.05);不同大写字母表示不同季节差异极显著(P<0.01)

Lowercase superscripts in the same row indicate significant differences (P<0.05)；upper-case superscripts indicate extremely significant differences (P<0.01)

2.2差异表达基因的筛选

利用RNA-seq技术检测试验组和对照组(SBL∶BL)间差异表达基因，结果发现，试验组比对照组有显著差异表达基因84个，其中56个基因表达上调，28个基因表达下调(Fold change>1.5，P<0.05)，部分结果见表3。

2.3候选基因的GO功能分类

利用DAVID Bioinformatics Resources 6.7中的Functional Annotation Chart模块将差异表达基因按照3大功能板块，包括生物学过程(BP)、细胞组分(CC)和分子功能(MF)进行限定和描述。结果(图1)发现，共涉及到215条GO Terms(P<0.05)。差异表达基因主要与应激反应、免疫反应、脂多糖代谢、RNA加工拼接及蛋白磷酸化相关。脂多糖代谢与炎症反应和凋亡相关。FOS、JUN、NRADD主要参与应激反应、免疫反应和蛋白磷酸化过程。DDX5主要参与RNA加工拼接过程。在细胞组成分中，差异表达基因主要参与膜细胞器的形成，在分子功能中，差异表达基因主要与结合功能相关。

2.4候选基因的KEGG通路富集分析

在生物体内，不同基因相互协调行使其生物学功能，根据KEGG分析有助于进一步分析差异表达基因中存在的显著性富集的代谢通路注释。本研究使用 DAVID Bioinformatics Resources 6.7软件中的“KEGG Functional Annotation Chart”模块分别对差异显著基因进行通路富集分析，结果共得到90条通路，其中显著富集的通路有18条，主要参与到剪接体、炎症性肠病、雌激素及TNF(肿瘤坏死因子)等通路中(图2)。

2.5差异表达基因的RT-PCR初步验证

选出8个差异显著基因进行荧光定量PCR检测，结果见图3。NRADD、FOS、JUN、DDX5、DUSP1、UBE2B为上调差异表达基因，KLK1和MMP9为下调差异表达基因。且定量结果与RNA-seq结果趋势一致，证明转录组在筛选差异表达基因方面的可靠性。

3讨论

表3　转录组测序前10个上下调基因

图1　GO分析结果Fig.1　The GO analysis result

图2　KEGG分析Fig.2　The KEGG analysis result

种公牛精液品质会随着季节变化而变化。本研究结果显示，与其他3个季节相比，夏季精子活力较差，畸形率较高，密度较低。与其他季节相比，春冬季节精子活力较高，畸形率最低，密度较高，春冬季节相比，春季精液各指标略优于冬季。夏季是槟榔江公牛精液品质最差的季节，秋季次之，冬春季节差异不显著。张汝等[10]在对槟榔江水牛精液品质的研究中发现，槟榔江水牛冷冻精液解冻活力和精子密度春、秋、冬季极显著高于夏季，采精量夏季低于秋季，夏季与冬季的差异不显著，这与本研究结果相一致。而夏季对公牛精液品质的影响可能与其温度湿度较高引起机体一系列的变化有关系。桑润滋等[11]研究表明，夏季温度较高时，容易引起公牛热应激，可造成荷斯坦种公牛精液品质显著下降，使原精活力、冻精活力、精子密度、活精子百分率和顶体完整率分别比春季下降3.13%、11.1%、34.8%、15.0%和17.8%，精子畸形率比春季上升25.5%。在精子生成过程中，热诱导能引起生精细胞凋亡，睾丸内精母细胞和圆形精子细胞对温度升高最敏感，是凋亡的主要细胞[12]，夏季高温刺激，有可能引起生精细胞的逐渐调亡，最终造成精子密度降低，活力较差。

1.NRADD;2.JUN;3.DUSPI;4.KLK1;5.FOS;6.DDX5;7.UPF2B;8.MMP9A.RT-PCR结果；B.RNA-seq结果A.RT-PCR result；B.RNA-seq result图3　RT-PCR检测Fig.3　Confirmation of RT-PCR

本试验选取具有代表性的夏春季节槟榔江水牛作为试验组和对照组，通过转录组测序，并在测序评估的基础上，共筛选到84个差异表达基因，包含56个上调基因和28个表达下调基因。对差异表达基因进行筛选、功能注释和富集分析后发现，参与应激反应、免疫反应相关基因HSF1、FOS和JUN在夏季试验组中表达量显著高于春季对照组。其中NRADD，又叫HSF1，属于热休克转录因子[13]，是最早被发现、最具代表性、最重要的热休克转录因子[14]。HSF1在热应激反应中的主要功能是在热休克基因的表达过程中与相应启动子结合，启动基因的转录过程，最终促进热休克蛋白(Heat shock proteins，HSP)的表达并使其表达增加[15]。本研究中夏季槟榔江水牛外周血中HSF1表达量较春季显著上调，且槟榔江水牛在夏季精子数量减少，活力下降，可能与夏季高温高湿气候引起槟榔江水牛出现热应激反应，并激活了热应激相关生物学过程和通路，导致HSF1表达量上升，以期对抗温度升高的不良环境对自身造成的损害。由此，笔者推断HSF1与精液品质相关。HSF1转基因小鼠生精障碍的机制为HSF1激活热休克蛋白(HSP60)和热休克蛋白(HSP105)，后者诱导着生殖细胞的大量凋亡，致使雄性小鼠无成熟精子，且研究表明由热诱导造成的雄性不育多数与HSF1的活化相关[16]。且HSF1能够启动雄性生殖细胞中的细胞程序化死亡，将雄性小鼠生精细胞阻滞在减数分裂期[17]。此外，J.Korfanty等[18]在小鼠模型中发现HSF1的过度表达引起精母细胞的大量凋亡并导致雄性不育。

FOS和JUN均参与TNF信号通路和雌激素信号通路。TNF信号通路主要与信号传导，细胞凋亡等生物学过程相关。包括热应激、趋化因子和炎症因子等都能够激活该通路。本研究中槟榔江水牛夏季精液品质下降可能是与TNF信号通路的激活相关。该通路的激活引起FOS和JUN基因的上调，而这些基因的上调能够对抗细胞的凋亡，降低对精液品质的影响。由此推断FOS和JUN基因上调可能与精液品质相关。K.M.Robertson等[19]研究表明，雌激素具有防止未变态精子细胞凋亡，维持成熟精子功能的作用。FOS(也叫C-FOS)和JUN(也叫C-JUN)是激活蛋白-1 (Activator protein-1，AP-1)类转录因子，属于原癌基因，其蛋白产物作为核转录因子在细胞增殖和分化中发挥重要作用。原癌基因对精液品质的影响，一方面是通过原癌基因转录产物的特异性或不同形式，另一方面是原癌基因在精子发生过程中不同水平的表达[20]。研究发现，多种应激可影响C-FOS和C-JUN的表达，从而诱导C-JUN氨基末端激酶(C-JUN N-terminal kinase，JNK)和p38丝裂原活化蛋白酶(p38 mitogen-acdvated protein kinase，p38)的表达，最终调控细胞的凋亡[21]。雌激素、促性腺激素及多种生长因子促进间质细胞的生长发育和睾酮的分泌是通过FOS和JUN家族的癌基因表达增强而实现的。人绒毛膜促性腺激素(hCG)和多种生长因子能刺激睾丸间质细胞的生长发育和睾酮分泌中均伴有原癌基因(包括FOS)的表达增强，且已证实C-FOS在睾丸组织中有广泛的表达，且在生精周期中有周期特异性[22]。本研究中与春季对照组相比较，夏季槟榔江水牛外周血中C-FOS和C-JUN表达量显著提高，可以推测这与槟榔江公牛在夏季高温高湿环境下发生应激反应，以期对抗不良刺激相关。C-FOS和C-JUN的作用可能与诱导精原细胞增殖、防止生精细胞凋亡、维持成熟精子功能有关。

细胞内的生物学过程，例如代谢、膜运输、信号转导和细胞的生长周期等的完成都是基于各个基因间系统地相互协调而进行的。通过对基因进行KEGG和GO富集分析，可以更深入地了解这些基因的生物学功能[23]。

4结论

本研究结果表明，将HSF1、FOS和JUN基因作为不同季节下影响槟榔江公牛精液品质的候选基因，对于培育抗应激品种，对养牛业持续、健康和稳定发展具有重要的理论和现实意义。

参考文献(References)：

[1]HAMILTON T R，MENDES C M，DE CASTRO L S，et al.Evaluation of lasting effects of heat stress on sperm profile and oxidative status of ram semen and epididymal sperm[J].OxidMedCellLongev，2016，2016：1687657.

[2]KOIVISTO M B，COSTA M T，PERRI S H，et al.The effect of season on semen characteristics and freezability inBosindicusandBostaurusbulls in the southeastern region of Brazil[J].ReprodDomestAnim，2009，44(4)：587-592.

[3]HELBIG L，WOODBURY M R，HAIGH J C，et al.The seasonal fertility of North Ameriman bison (Bisonbison) bulls[J].AnimReprodSci，2007，97(3-4)：265-277.

[4]ROCHA A，RANDEL R D，BROUSSARD J R，et al.High environmental temperature and humidity decrease oocyte quality inBostaurusbut not inBosindicuscows[J].Theriogenology，1998，49(3)：657-665.

[5]王宝全，孔伟，马志远，等.高温对种公牛精液品质的影响研究[J].中国牛业科学，2013，39(5)：13-15.

WANG B Q，KONG W，MA Z Y，et al.Influence of high temperature on semen quality[J].ChinaCattleScience，2013，39(5)：13-15.(in Chinese)

[6]MAYAHI S，MAMOUEI M，TABATABAEI S，et al.Reproductive characteristics and thyroidal function in relation with season in Khuzestan buffalo (Bubalusbubalis) bulls[J].VetResForum，2014，5(3)：201-205.

[7]刘红亮，郑丽明，刘青青，等.非模式生物转录组研究[J].遗传，2013，35(8)：955-970.

LIU H L，ZHENG L M，LIU Q Q，et al.Studies on the transcriptomes of non-model organisms[J].Hereditas，2013，35(8):955-970.(in Chinese)

[8]YAN L，YANG M，GUO H，et al.Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells[J].NatStructMolBiol，2013，20(9)：1131-1139.

[9]贾昌路，张瑶，朱玲，等.转录组测序技术在生物测序中的应用研究进展[J].分子植物育种，2015，13(10)：2388-2394.

JIA C L，ZHANG Y，ZHU L，et al.Application progress of transcriptome sequencing technology in biological sequencing[J].MolecularPlantBreeding，2015，13(10)：2388-2394.(in Chinese)

[10]张汝，赵家才，毛翔光，等.槟榔江水牛精液品质的研究[J].云南畜牧兽医，2013(5)：20-22.

ZHANG R，ZHAO J C，MAO X G，et al.The semen quality study ofBosbubalusbulls[J].YunnanAnimalHusbandry&VeterinaryMedicine，2013(5)：20-22.(in Chinese)

[11]桑润滋，田树军，李俊杰，等.热应激对荷斯坦种公牛精液品质的影响及机理[J].中国兽医学报，2002，22(4)：407-410.

SANG R Z，TIAN S J，LI J J，et al.Influence of heat stress on semen quality of Holstein bull and approach to its mechanism[J].ChineseJournalofVeterinaryScience，2002，22(4)：407-410.(in Chinese)

[12]DURAIRAJANAYAGAM D，AGARWAL A，ONG C.Causes，effects and molecular mechanisms of testicular heat stress[J].ReprodBiomedOnline，2015，30(1)：14-27.

[13]AKERFELT M，MORIMOTO R I，SISTONEN L.Heat shock factors：integrators of cell stress，development and lifespan[J].NatRevMolCellBiol，2010，11(8)：545-555.

[14]XU D，SUN L，LIU S，et al.Molecular cloning of hsfll and hsbpl cDNAs，and the expression of hsfl，hsbpl and hsp70 under heat stress in the sea cucumber Apostichopus japonicus[J].CompBiochemPhysiolBBiochemMolBiol，2016，198：1-9.

[15]ABANE R，MEZGER V.Roles of heat shock factors in gametogenesis and development[J].FEBSJ，2010，277(20)：4150-4172.

[16]HAYASHIDA N，INOUYE S，FUJIMOTO M，et al.A novel HSF1-mediated death pathway that is suppressed by heat shock proteins[J].EMBOJ，2006，25(20)：4773-4783.

[17]NAKAI A，SUZUKI M，TANABE M.Arrest of spermatogenesis in mice expressing an active heat shock transcription factor 1[J].EMBOJ，2000，19(7)：1545-1554.

[18]KORFANTY J，STOKOWY T，WIDLAK P，et al.Crosstalk between HSF1 and HSF2 during the heat shock response in mouse testes[J].IntJBiochemCellBiol，2014，57：76-83.

[19]ROBERTSON K M，O’DONNELL L，JONES M E，et al.Impairment of spermatogenesis in mice lacking a functional aromatase (cyp 19) gene[J].ProcNatlAcadSciUSA，1999，96(14)：7986-7991.

[20]况海斌，方廉.原癌基因与睾丸功能的关系[J].中华男科学，2003，9(5)：377-380.

KUANG H B，FANG L.The relationship between proto-oncogene and testicular function[J].NationalJournalofAndrology，2003，9 (5)：377-380.(in Chinese)

[21]WESTON C R，DAVIS R J.The JNK signal transduction pathway[J].CurrOpinCellBiol，2007，19(2)：142-149.

[22]HALL S H，BERTHELON M C，AVALLET O，et al.Regulation of c-fos，c-jun，jun-B，and c-myc messenger ribonucleic acids by gonadotropin and growth factors in cultured pig Leydig cell[J].Endocrinology，1991，129(3)：1243-1249.

[23]杜琛，付绍印，高鸿雁，等.绒山羊肌内脂肪细胞成熟前后比较转录组分析[J].畜牧兽医学报，2014，45(5)：714-721.

DU C，FU S Y，GAO H Y，et al.Transcriptome analysis of intramuscular preadipocytes and mature adipocyte in cashmere goats[J].ActaVeterinariaetZootechincaSinica，2014，45(5)：714-721.(in Chinese)

(编辑 程金华)

doi:10.11843/j.issn.0366-6964.2016.07.009

收稿日期：2016-03-22

基金项目：国家自然科学基金项目(31372290);“十二五”国家科技计划项目(2012BAD12B10)

作者简介：耿红红(1989-)，女，山西孝义人，硕士生，主要从事动物生殖与生理调控研究，E-mail：15150537075@163.com *通信作者：王根林，教授，E-mail:genlwang@sina.com

中图分类号：S823.8+3

文献标志码：A

文章编号:0366-6964(2016)07-1373-08

Semen Quality Analysis of Betelnut-jiang Buffalo in Different Seasons by RNA-seq

GENG Hong-hong，ZHANG Jing-yang，LI Lian，WANG Gen-lin*

(CollegeofAnimalScienceandTechnology，NanjingAgriculturalUniversity，Nanjing210095，China)

Abstract:To determine the key candidate genes related to semen quality in different seasons，the mRNA profiling of Chinese Betelnut-jiang buffalo were analyzed by RNA-seq.Based on similarities in growth status，5 Betelnut-jiang buffalo bulls were selected and their semen quality were tested (including ejaculate volume，fresh sperm motility，sperm density and abnormality rate).According to the semen quality，the samples of summer (treated as SBL group) and spring (as BL group) were selected for RNA-seq.We collected the whole blood in jugular vein of Betelnut-jiang buffalo bulls，extracted the mRNAs and mixed them within group.Then，the mRNA profiling was used to measure the mRNAs in Betelnut-jiang buffalo bulls.Gene expression levels with >1.5 fold change and P value < 0.05 between the BL and SBL groups were further analyzed.Furthermore，GO and KEGG were used to analyze the genes enrichment.Finally，RT-PCR was used to confirm the RNA-seq result.The results indicated 84 differentially expressed genes.Among these，56 genes were up-regulated and 28 genes were down-regulated in SBL group.The GO search revealed that these differential expression genes were mainly related to stress reaction and immune response.And the KEGG selected the TNF signaling pathway and T cell receptor signaling pathway.Neurotrophin receptor alike death domain(NRADD，also named heat shock factor protein 1(HSF1))，localisation of FOS (FOS) and jun proto-oncogene (JUN)，were strongly associated with semen quality in different seasons which could be used as target genes for further study.

Key words:RNA-seq；Betelnut-jiang buffalo；differently expressed genes；semen quality