0 引言
【研究意义】番鸭属的公番鸭与河鸭属的母家鸭亲缘关系远,其属间杂交后代为无繁殖能力的骡鸭,俗称半番鸭。半番鸭抗逆性强,饲料报酬高,肉质细腻,具有很强的杂种优势,但其公母外型差异不显著,性器官、性行为分化不明显,无法正常繁育后代,无实际的种用价值。实践研究发现,半番鸭并非完全意义上的不育,它们拥有精巢组织,极少数可产生少量精液,少数可正常交配,雄性偶见爬跨,雌性偶见产蛋,这可能涉及到基因表达、调控以及性细胞形成方式等。【前人研究进展】远缘杂交不育是一个很复杂的生物学现象,其遗传基础较复杂,可能受众多的基因及其产物的调控。关于鸭类属间杂交不育的原因,檀俊秩和陈晖等[1]认为半番鸭父母本的染色体核型不一致,减数分裂受到破坏,进而影响性腺轴发育,最终导致其后代不育。刘军须等[2]研究认为大鼠不育性状由常染色体上单一隐性基因控制,呈隐性遗传。张庆波等[3]认为DAZL(Deleted in AZoospermia-like)基因是牛精子发生的重要调控因子,其突变或表达缺乏将导致雄性不育。以Illumina为基础的转录组测序(RNA-Seq)技术被广泛的应用到差异基因的检测及功能注释,在畜禽方面也应用广泛,如鉴定鸭黑、白羽相关基因[4],筛选野鸡和家鸡肉质基因[5],分析绍兴鸭青壳性状相关基因[6]等研究。生殖、发育机制等一直是生物学热点,BAUERSACHS等[7]通过RNA-Seq技术筛选到不同物种间繁殖相关差异基因,FIEDLER等[8]对海兔不同发育阶段进行转录组测序,张伟[9]构建了中华绒螯蟹精巢组织文库,筛选雄性生殖调控关键基因,钟志君[10]研究比较了成年藏猪的精巢和卵巢转录组表达谱,张升利等[11]分析了长尾草金鱼成熟期精巢和卵巢转录组差异表达基因,XU等[12]利用转录组测序方法挖掘出与卵泡发育相关基因,朱志明等[13]探明了山麻鸭开产期和产蛋高峰期卵巢组织的转录组差异。这些成果为研究鸭不育性状和繁殖性能的遗传机制奠定了基础。【本研究切入点】虽然 RNA-Seq技术已被应用到鸭的性状研究,但由于半番鸭不育性状的特殊性,其与番鸭或家鸭精巢组织的转录组比较分析尚未有报道,其不育遗传基础需进一步深入研究。【拟解决的关键问题】本研究通过对半番鸭和番鸭个体精巢组织的转录组比较,对其差异表达基因进行筛选,并进行GO(gene ontology)与KEGG(kyoto encyclopediaof genes and genomes)功能富集,分析通路功能并探索精巢分化相关的差异表达基因,通过转录组数据挖掘半番鸭雄性不育相关的基因,为后续半番鸭雄性不育形成的遗传机制研究提供线索。1 材料与方法
1.1 试验材料
本研究所用的半番鸭和番鸭公鸭由福建省农业科学院畜牧兽医研究所动物房提供,试验于3—6月进行,各试验鸭饲养管理条件一致。180日龄性成熟时,对公鸭进行采精训练。210日龄禁饲12 h后,每个品种分别取2 只个体用于转录组测序,其中半番鸭公鸭为无爬跨行为个体,番鸭公鸭为正常个体。按照国家实验动物处理行为准则屠宰,取精巢组织,置于-80℃备用。1.2 RNA 提取及转录组测序
利用RNA easy Lipid Tissue Mini Kit(QIAGEN, Germany)提取每个个体总RNA,单个建池。采用Nanodrop、Qubit 2.0和Aglient 2100方法检测各RNA样品的纯度、浓度及完整性等。构建文库,Qubit2.0和Agilent 2100分别对文库的浓度和插入片段大小(Insert Size)进行检测,QRT-PCR对文库的有效浓度进行准确定量,以保证文库质量。基于边合成边测序(Sequencing By Synthesis,SBS)技术,利用Illumina HiSeq 2500(Illumina, America)平台进行高通量测序, 测序读长为PE150。1.3 测序数据处理
测序数据去除接头及低质量数据后,利用Trinity软件将数据组装成转录本,进行转录注释及表达量的计算。利用DESeq进行基因的差异表达分析,绘制差异表达基因火山图,并进行聚类分析。1.4 GO注释、KEGG通路分析
利用GO数据库对差异表达基因进行功能注释,采用COG(cluster of orthologous groups of proteins)对差异表达基因进行分类统计,运用KEGG数据库进行通路分析,均以P<0.05作为显著性富集标准。1.5 Real-time PCR检测
采用Primer Premier 6.0和Beacon designer 7.8软件设计荧光定量PCR引物,然后由生物工程(上海)股份有限公司负责合成,引物序列如表1,QRT-PCR扩增体系和反应条件如表2。以甘油醛-3-磷酸脱氢酶基因(glyceraldehyde-3-phosphate dehydrogenase, GAPDH)做内参,采用2-ΔΔCt法计算基因的相对表达量。以P<0.05作为显著性标准。Table 1
表1
表1荧光定量引物序列
Table 1Real-time PCR primers and conditions
基因名称 Gene | 全称 Full name | 引物序列(5'→3')Primer sequences(5' to 3') |
---|---|---|
GAPDH | 甘油醛-3-磷酸脱氢酶基因 Glyceraldehyde-3-phosphate dehydrogenase | GGAGCTGCCCAGAACATTATC |
GCAGGTCAGGTCCACGACA | ||
CPLA2 | 胞质型磷脂酶A2 Cytosolic phospholipase A2 epsilon-like | CACAGATCAGCAGCAGGCACTAA |
CTGCACGAATGAAGGCACCATA | ||
DUSP3 | 双特异性磷酸酶3 Dual specificity phosphatase 3 | CTCCTACCCTGGTCATCGCCTA |
CTTCTGCCGGACAGTGACCAA | ||
FAS | 肿瘤坏死因子受体超家族成员6 Tumor necrosis factor receptor superfamily member 6 | CCCAGGTGAAACACTCCCTGA |
GTCGATAACGGACAAACTTCTTGAC | ||
FGF3 | 成纤维细胞生长因子 Fibroblast growth factor | GGTCCGACTGTTCCACACAAATAC |
CTCCTTCTTCGCTGGTTCTTACTG | ||
GRB2 | 生长因子受体结合蛋白2 Growth factor receptor-bound protein 2 | GGGGCTGGCAAGTATTTCCTCT |
GGCACCTGTTCAATGTCCCGTA |
新窗口打开
Table 2
表2
表2荧光定量PCR反应体系及条件
Table 2Real-time PCR reaction system and conditions
20 µL体系 20 µL system | |
---|---|
SDW | 8.0 µL |
Power SYBR® Green Master Mix | 10.0 µL |
Forward primer (10 μmol·L-1) | 0.5 µL |
Reverse primer (10 μmol·L-1) | 0.5 µL |
cDNA | 1.0 µL |
新窗口打开
2 结果
2.1 半番鸭与番鸭精巢组织转录组测序
本试验共构建半番鸭和番鸭4个精巢组织的转录组文库。测序数据去除接头以及低质量数据后共获得288 008 582个高质量数据。从表3中可以看出,本试验共获得43.84Gb的Clean Data,各样品Clean Data均达到6.29Gb。另外,各样品GC含量均不小于51.88%,Q30(Clean Data质量不小于30的碱基所占的百分比)全在91.36%以上,该结果表明测序质量可靠,构建文库可用于后续分析。数据组装后共获得193 535条Unigene,其中长度在1kb以上的Unigene有46 175条。Table 3
表3
表3转录组数据组装情况
Table 3Summary of the sequencing data assembly
样品 Sample name | 高质量数据 Clean reads | 待分析数据 Clean bases | %≥Q30 | GC含量 GC content | 比对数据 Mapped reads | 比对率 Mapped ratio |
---|---|---|---|---|---|---|
T1 | 29 697 861 | 8 806 233 668 | 91.36% | 54.79% | 20 950 297 | 70.54% |
T2 | 23 809 701 | 7 066 621 636 | 91.36% | 52.72% | 16 466 974 | 69.16% |
T5 | 21 859 316 | 6 489 813 952 | 91.48% | 51.88% | 14 983 171 | 68.54% |
T6 | 21 217 125 | 6 293 982 810 | 91.36% | 52.43% | 14 681 967 | 69.20% |
新窗口打开
2.2 差异表达基因的筛选及聚类分析
FPKM(Fragments Per Kilobase of transcript per Million mapped reads)是每百万Reads中来自比对到某一基因每千碱基长度的Reads数目,是转录组测序数据分析中常用的基因表达水平估算方法。将皮尔逊相关系数r(Pearson’s Correlation Coefficient)作为生物学重复相关性的评估指标[14]。r2越接近1,说明两个重复样品相关性越强。对同一条件的每一对生物学重复样品的基因表达量做相关性图,相关性图见图1。半番鸭组内2个不同重复样品和番鸭组内2个不同重复样品的r2值均大于0.9,而两组间r2值均小于0.6,说明了试验的可靠性和可重复性都很高。显示原图|下载原图ZIP|生成PPT
图1四样品基因表达量相关性图
-->Fig. 1Correlation heat map of gene expression level in 4 samples
-->
以FDR(False Discovery Rate)作为差异表达基因筛选的关键指标, 将FDR小于0.001且差异倍数FC(Fold Change)大于等于2作为两品种鸭个体间差异表达基因显著的筛选标准,共鉴定出3 597个基因,其中包括1 194个表达上调的差异基因,主要有(FDR由小到大排列)c86758.graph_c1(NADH脱氢酶亚基4),c86758.graph_c0(外周型苯二氮卓受体相关蛋白1)等差异基因,以及2 403个表达下调基因,包含c243330.graph_c0(DnaJ同源B亚家族成员8),c277017.graph_c0(凋亡因子BCL-2蛋白14)等。进一步分析发现,显著差异基因中含有一些与繁殖性能相关的基因(表4),例如,成纤维细胞生长因子(fibroblast growth factor,FGF)、c-Jun氨基末端激酶(c-Jun N-terminal kinase,JNK)、细胞外信号调控的蛋白激酶5(ERK5)、蛋白激酶A(protein kinase A, PKA)、生长因子受体结合蛋白2(growth factor receptor-bound protein 2,Grb2)、丝裂原活化蛋白激酶7,部分(mitogen-activated protein kinase 7-like, partial,BMK)等下调基因,肿瘤坏死因子受体超家族成员6(tumor necrosis factor receptor superfamily member 6, FAS)、双特异性磷酸酶3(dual specificity phosphatase 3, DUSP3)、L型电压依赖性钙通道α1c亚单位(alpha-1c-like voltage-dependent L-type calcium channel subunit alpha-1C-like, CACNA1C)、胞质型磷脂酶A2(cytosolic phospholipase A2 epsilon-like, CPLA2)等上调基因。从火山图(图2)能够快速查看两组样品间表达的差异水表平分布情况。通过MA图(图3)可以直观地查看两组样品中基因的表达丰度和差异倍数的整体分布。
Table 4
表4
表4差异表达基因(部分)
Table 4Differentially expressed genes (parts)
基因ID Gene_ID | 基因名称 Genes name | FPKM | 半番/黑番 Regulated | |||
---|---|---|---|---|---|---|
T01 | T02 | T05 | T06 | |||
c267342.graph_c0 | DUSP3 | 1.441 | 1.635 | 11.439 | 8.601 | up |
c237612.graph_c0 | FAS | 0.025 | 0.251 | 0.865 | 0.635 | up |
c268454.graph_c0 | FGF3 | 16.538 | 16.349 | 3.984 | 4.944 | down |
c255957.graph_c0 | CPLA2 | 0.056 | 0.142 | 0.943 | 0.695 | up |
c247033.graph_c0 | GRB2 | 51.821 | 56.988 | 19.466 | 18.698 | down |
新窗口打开
显示原图|下载原图ZIP|生成PPT
图2差异表达基因火山图差异表达火山图中的每一个点表示一个基因。图中绿色和红色的点代表有显著性表达差异的基因,绿色代表基因表达量下调,红色代表基因表达量上调,黑色的点代表无显著性表达差异的基因。下同
-->Fig. 2Volcano plot of differentially expressed genes:Each dot represents a gene. Green dot: Significantly down-regulated genes; Red dot: Significantly up- regulated genes; Black dot: Genes without differential expression. The same as below
-->
显示原图|下载原图ZIP|生成PPT
图3差异表达基因MA图
-->Fig. 3MA plot of differentially expressed genes
-->
对筛选出来的差异表达基因进行聚类分析(图4),发现同一鸭品种的2 个生物学重复聚到了一起,而不同鸭品种间基因表达模式则出现分离,表明本研究所用样本生物学重复性较好,且样本分组较合理。
显示原图|下载原图ZIP|生成PPT
图4显著差异表达基因聚类分析图中不同的列代表不同的样品,不同的行代表不同的基因。颜色代表了基因在样品中的表达量FPKM的对数值
-->Fig. 4Heat map of the differentially expressed genes Columns indicate individual samples, the row represents each differentially expressed genes. The color scale represents log10 (FPKM)
-->
2.3 差异表达基因GO功能富集
本研究运用GO数据库对具有同源比对的差异基因进行生物学过程(biological process)、分子功能(molecular function)与细胞组分(cellular component)三方面的注释。将比对得到的1 381个显著差异基因进行功能注释,382个差异基因在GO分类中有功能意义。上述3个功能被区分为更具体的61个类别,分别包括了22、19和20个功能亚分类。由图5可知,在生物功能的组分中,两组间的差异表达基因在细胞过程(cellular process, GO: 0009987)与单一的生物过程(single-organism process, GO: 0044699)中数目比例最大。在细胞功能中,细胞(cell, GO: 0005623)与细胞部分(cell part, GO: 0044464)数目最多。在分子功能分类中,差异基因在结合(binding,GO: 0005488)中所占的比例最高,催化活性(catalytic activity, GO: 0003824)次之。由图5可知,与发育繁殖相关的生物学过程有繁殖(reproduction, GO:0000003),发育过程(developmental process, GO:0048589)和繁殖过程(reproductive process, GO:0022414),涉及97个相关基因。2.4 差异表达基因COG分类
统计显著富集的Go term中包含的基因数,分析结果见图6。其中一般的功能预测(general function prediction only)基因数目最多,其次为信号传导机制(signal transduction mechanism),转录(Transcription)和复制、重组与修复(replication, recombination and repair)。2.5 差异表达基因KEGG注释
为确定差异基因参与的主要生化代谢途径和信号通路,对差异表达基因进行KEGG通路分析,结果显示差异表达基因共富集到50 条信号通路中。以P<0.05作为差异表达基因在该通路显著富集,共鉴定出17个通路显著富集:包括与丝裂原活化蛋白激酶(MAPK)信号转导通路(MAPK signaling pathway)、钙信号途径(calcium signaling pathway)、甘油酯代谢(glycerolipid metabolism)、紧密连接(gap junction)以及血管平滑肌收缩(Vascular smooth muscle contraction)等信号通路。对差异表达基因的注释结果按照KEGG中通路类型进行分类,分类图如图7所示,其中与生长发育、生殖过程相关的有MAPK信号转导通路和促性腺激素释放激素信号通路(gonadotropinreleasing hormone(GnRH)signaling pathway)等。显示原图|下载原图ZIP|生成PPT
图5差异表达基因GO注释1. 细胞过程Cellular process;2. 单一的生物过程single-organism process;3. 代谢过程Metabolic process;4. 生物调节Biological regulation;5. 刺激反应Response to stimulus;6. 多细胞生物过程Multicellular organismal process;7. 发育过程Developmental process;8. 定位Localization;9. 信号Signaling;10. 组织或生物起源细胞组件Cellular component organization or biogenesis;11. 免疫系统过程Immune system process;12. 多生物体过程Multi-organism process;13. 移动Locomotion;14. 繁殖过程Reproductive process;15. 生物粘附Biological adhesion;16. 繁殖Reproduction;17. 生长Growth;18. 有节奏的过程Rhythmic process;19. 激素分泌hormone secretion;20. biological phase生物相;21. 细胞聚集cell aggregation;22. 细胞杀伤Cell killing1. 细胞 Cell;2. 细胞部分Cell part;3. 细胞器Organelle;4. Membrane膜;5. 细胞器部分Organelle part;6. 高分子复合物Macromolecular complex;7. Membrane part膜部分;8. 膜包围内腔Membrane-enclosed lumen;9. 胞外区Extracellular region;10. 细胞外区域部分Extracellular region part;11. 细胞连接Cell junction;12. 突触Synapse;13. 细胞外基质Extracellular matrix;14. 突出部分Synapse part;15. 胶原三聚体collagen trimer;16. 细胞外基质部分Extracellular matrix part;17. 细胞核Nucleoid;18. 病毒体Virion;19. 病毒体部分Virion part1. 结合Binding;2. 催化活性Catalytic activity;3. 运输活性transporter activity;4. 分子传感器活性Molecular transducer activity;5. 核酸结合转录因子活性Nucleic acid binding transcription factor activity;6. 受体活性Receptor activity;7. 结构分子活性Structural molecule activity;8. 酶调节活性Enzyme regulator activity;9. 蛋白结合转录因子活性Protein binding transcription factor activity;10. 鸟嘌呤核苷酸交换因子活性guanyl-nucleotide exchange factor activity;11. 电子载体活性Electron carrier activity;12. 抗氧化活性Antioxidant activity;13. 通道调节活性Channel regulator activity;14. 受体调节活性Receptor regulator activity;15. 化学引诱物活性Chemoattractant activity;16. 转录调节因子活性Translation regulator activity;17. 金属伴侣活性Metallochaperone activity;18. 成形素活性Morphogen activity;19. 排斥活动Chemorepellent activity ;20. protein tag蛋白标签
-->Fig. 5GO annotation of differentially expressed genes
-->
显示原图|下载原图ZIP|生成PPT
图6差异表达基因COG注释分类统计图
-->Fig. 6COG annotation classification statistics of differentially expressed genes
-->
显示原图|下载原图ZIP|生成PPT
图7差异表达基因显著富集的KEGG 通路
-->Fig. 7List of KEGG pathway for differentially expressed genes
-->
挑选富集显著性最可靠(即Q值最小)的前20个通路以散点图的形式展示(图8),在该图中越靠近右上角的图形代表的通路,参考价值越大;反之亦然。由此可知甘油酯代谢(glycerolipid metabolism),钙信号途径(calcium signaling pathway)信号通路以及血管平滑肌收缩(vascular smooth muscle contraction)富集显著性较为可靠。
2.6 Real-time PCR验证转录组测序数据
为验证转录组测序结果,本研究选择了GRB2、CPLA2、FGF3、DUSP3及FAS基因进行Real-time PCR试验。结果表明除FAS基因整体表达低,数据不准确外,其他4个基因在番鸭和半番鸭个体间表达变化模式与转录组测序结果一致(表5),表明本研究利用转录组测序获得的数据较为准确。显示原图|下载原图ZIP|生成PPT
图8差异表达基因KEGG通路富集散点图
-->Fig. 8Enriched scatter map of differential expression gene KEGG pathway
-->
Table 5
表5
表5Real-time PCR 验证转录组测序数据
Table 5Validation of the RNA-seq expression data by Real-time PCR for selected genes
基因名称 Genes name | 转录组测序 RNA-seq | 荧光定量PCR Real-time PCR | ||
---|---|---|---|---|
表达倍数 (半番鸭/番鸭) log2FC (Mule duck/Muscovy duck) | FDR | 表达倍数 (半番鸭/番鸭) log2FC (Mule duck/ Muscovy duck) | P | |
FGF3 | -2.103 | 1.81×10-8 | -0.448 | 0.217 |
GRB2 | -1.734 | 1.80×10-6 | -0.334 | 0.124 |
FAS | 2.298 | 6.36×10-4 | 1.333 | 0.200 |
DUSP3 | 2.475 | 2.64×10-9 | 0.855 | 0.010 |
CPLA2 | 2.853 | 5.48×10-5 | 0.340 | 0.445 |
新窗口打开
3 讨论
近来, 高通量转录组测序技术取得了较大的研究进展[15-20]。雄性生殖系统的发育、分化过程是一个复杂的生理过程,涉及多种代谢过程,目前已有研究表明p53信号通路[21]、Wnt代谢通路[22]以及TGFβ通路[23-24]在雄性生殖发育过程中起重要作用。为了解鸭精巢的系统机制和半番鸭不育的特性,本研究对半番鸭和番鸭精巢组织进行转录组学测序及比对分析,共筛选获得3 597个差异表达基因,同时为进一步确定差异基因参与的主要生化代谢途径和信号通路进行KEGG分析,结果发现MAPK、甘油酯代谢以及钙信号途径等17个通路富集显著,其中与生殖过程相关的有GnRH信号通路和MAPK信号转导通路。MAPK信号通路是细胞间信号传递的重要通路[25],可参与细胞增殖、分化[26]、精子成熟[27]及凋亡[28-29]等动物生殖过程,被认为是精子发育的重要影响因素之一。本研究筛选到参与此通路的FGF、c- JNK、ERK5等10个下调基因,FAS、DUSP3等4个上调基因,以及CPLA2、CACN 2个混合型调节基因。在此通路筛选的差异表达基因中,FGF属于可促进成纤维细胞生长的多肽家族。VALVE等[30]研究发现FGF8基因在成年鼠睾丸精子发生的特定阶段表达,在鼠、牛卵巢中卵子发生的特定阶段也有表达。DVORAK等[31]研究表明FGF信号参与包括细胞增殖、迁移、分化等多种细胞应答,调控广泛的生理和病理过程。JNK在细胞分化、凋亡和疾病的发生均有重要作用,是MAPK信号转导系统的重要效应因子[32]。目前发现的JNK的底物有转录因子c-jun、ATF-2及P53等。陈丽莉等[33]发现JNK可调控黄鳝卵巢发育、凋亡及雄性发育的启动。ERK基因则主要参与动物生殖细胞增殖、分化以及凋亡等多个发育过程[34]。
GnRH信号通路可调控动物体内性腺轴生殖激素的分泌,下丘脑分泌产生的神经激素能促进垂体分泌促性腺物质的释放,并参与动物生殖调控[35]。本研究筛选到了参与此通路上的PKA、Grb2和BMK等下调基因,以及CPLA2、CACNA1C 2个混合调节基因。在此通路筛选的差异表达基因中,已有研究证实Grb2是信号转导途径中的一个重要成分,它最初是作为表皮生长因子受体(epidermal growth factor receptor, EGFR)和MAPK通路之间的衔接蛋白被发现的[36]。Grb2作为一种信号接头蛋白,参与细胞信号转导过程,具有促进细胞增殖、细胞生长、细胞分化等功能[37]。PLA2参与雄性生殖过程, 其与精子获能、顶体反应以及精卵融合等过程密切相关,并受到各种信号通路的调节,不同通路相互协调。精子PLA2活性的激活及其调控机制受G蛋白受体可介导[38]。
另外,JNK、Crb2、CACN等基因同时参与上述两个通路,如PKA、蛋白激酶C(protein kinase C, PKC)及Ca2+也在生理性顶体反应的精子信号转导中调节PLA2 的活性,此外,钙信号途径以及甘油酯代谢、紧密连接以及血管平滑肌收缩等信号通路也富集显著,说明各种发挥不同作用的信号传导通路交织在一起, 组成一个复杂庞大的细胞信号传导网络系统。关于钙信号途径,Ca2+作为细胞内重要第二信使,通过精子内Ca2+浓度调节精子生理活动。SANTI等[39]也发现 Ca2+信号通路与精子的形成有着密切的关系。甘油酯代谢、紧密连接以及血管平滑肌等信号通路在鸭精巢组织中的作用及其相互调节机制将进一步探究。
4 结论
本研究通过半番鸭和番鸭转录组比较,首次在转录组水平上筛选出与繁殖性能相关的差异表达基因,如PKA、Grb2、BMK、CPLA2、FGF、JNK及ERK5等,并进一步证实了GnRH 信号通路和MAPK信号通路在公鸭生殖活动中发挥了重要作用。信号通路的分析为鸭精巢组织的信号调控提供了线索,但还需对涉及这些通路的相关基因进行深入的生物信息学分析和验证,进一步明确繁殖性状的基因表达谱和调控模式,该研究结果可为今后探索半番鸭生殖系统的分化机理提供可靠的参考依据。The authors have declared that no competing interests exist.