中华蜜蜂精巢和卵巢差异表达基因分析

成付平, 袁 芳, 席芳贵, 秦凯鑫, 胡小芬, 王子龙,*

(1. 江西农业大学蜜蜂研究所, 南昌 330045; 2. 江西省蜜蜂生物学与饲养重点实验室, 南昌330045;3. 江西省养蜂研究所, 南昌330052; 4. 江西农业大学动物科学技术学院, 南昌 330045)

中华蜜蜂Apisceranacerana(简称中蜂)是东方蜜蜂Apiscerana的指名亚种，由于它在抵御寒冷、采集零散蜜源、抵抗胡蜂及有节制地产卵等方面有着显著的优势(曾志将, 2020)，使它成为我国饲养的重要蜂种。近些年，随着研究手段的不断更新与深入，已由描述性研究转到分子水平上来阐述蜜蜂发育机理。健康的蜂群由蜂王、工蜂、雄蜂组成，它们的性别决定是一种特殊的单双倍体模式，雄性是由未受精卵发育而来的单倍体，雌性是由受精卵发育而来的二倍体(Beyeetal., 2003)。在一个蜂群中，蜂王和雄蜂都拥有发育完全的生殖腺，但它们的性成熟时间是不同步的。蜂王在羽化一周后达到性成熟，而雄蜂则需要大约12 d才能达到性成熟(曾志将, 2017)。

蜜蜂的性别发育包括两个方面，分别是性别决定和性别分化。在蜜蜂中，性别是由csd→fem→dsx这个级联反应来决定的(Wilkins, 1995; Hasselmannetal., 2008, 2010)。在这个级联调控过程中，csd是调节fem进行雌特异性剪接的初始信号(Hasselmannetal., 2008, 2010)，之后fem又控制着dsx的雌特异性剪接。在多细胞动物中，dsx是性别决定级联反应中最保守的性别决定基因之一(Wilkins, 1995)。这一级联反应的下游基因充当性别调节因子，形成性别分化途径，进而导致雌雄特异性表型。

虽然目前已经鉴定出蜜蜂性别级联的主效基因，但对蜜蜂生殖活动和性腺发育的研究较少。因此，有必要了解蜜蜂性腺分化和配子发生的分子机制。本研究对中华蜜蜂的卵巢和精巢进行了转录组测序分析。我们的结果揭示了雄性和雌性性腺的基因表达差异。本研究的结果有助于理解中华蜜蜂性别决定和分化的分子机制，并为未来的研究提供有价值的基因表达信息。

1 材料与方法

1.1 供试蜜蜂及采样

本研究所用的实验蜂群为江西省养蜂研究所饲养的中华蜜蜂蜂群。选择两个群势较强的中华蜜蜂蜂群，采用标准的人工育王方法(曾志将, 2009)培育中华蜜蜂蜂王，待蜂王出房后随机挑选个体发育优良的处女王剪掉翅后放于巢房中培育，到蜂王出房12 d时进行采样，每2头蜂王的卵巢作为一个样品，每群采集一个样品，共采集2个生物学重复。采用同样的蜂群，将蜂王控制在一张干净的雄蜂巢脾中产卵，12 h之后放出蜂王，并将该雄蜂巢脾放入继箱群中继续孵化发育。等雄蜂羽化出房后进行标记，再放入原实验群中进一步培育，到雄蜂出房12 d时进行采样，每10头雄蜂精巢作为一个样品，每群采集一个样品，共采集2个生物学重复。取样时，在显微镜下对雄蜂和蜂王腹部分别进行解剖，取出雄蜂精巢及蜂王卵巢，将采集好的精巢和卵巢样品用液氮速冻后装在1.5 mL EP管中，保存于液氮中。

1.2 cDNA文库的创建与测序

按照标准的方法分别提取1.1节每个精巢和卵巢样品的总RNA，对RNA质量进行检测,以确保得到合格样本并且能够进行转录组测序。待实验样品总RNA检测结果均合格后，进行文库构建，主要操作流程为：用带有Oligo(dT)的磁珠富集样品mRNA；加入Fragmentation Buffer随机打断mRNA片段；用六碱基随机引物(random hexamers)合成cDNA第1链，然后加入酶Buffer, dNTPs, RNase H和DNA 聚合酶I合成cDNA第2链，反转录结束后，通过AMPure XP Beads纯化cDNA；对纯化过后的cDNA进行末端修复、加poly(A)尾并连接测序接头，然后用AMPure XP Beads进行片段大小选择；最后通过PCR扩增富集得到cDNA文库。用Illumina高通量测序平台对cDNA文库进行测序。

1.3 测序数据与参考基因组比对

利用FastQC软件对数据进行分析筛选，去除含有接头的读段和低质量的读段(包括去除N的比例大于10%的读段；去除质量值Q≤10的碱基数占整条读段的50%以上的读段)，最终筛选得到高质量的过滤后数据。Illumina HiSeq碱基质量值(quality score或Q-score)的计算通常使用Phred碱基质量值公式(Ewingetal., 1998)：Q-score=-10×log10P，P为碱基识别出错的概率。

利用TopHat2(Kimetal., 2013)软件将获得的过滤后读段与中华蜜蜂最新的基因组序列(ftp:∥ftp.ncbi.nlm.nih.gov/genomes/genbank/invertebrate/Apis_cerana/latest_assembly_versions/GCA_011100585.1_ASM1110058v1)进行序列比对，获取过滤后读段在参考基因组和基因上的位置信息，以及测序样品特有的序列特征信息。

1.4 基因表达水平确定及标准化

使用Cufflinks软件的Cuffquant和Cuffnorm组件，通过匹配上的读段在基因上的位置信息，对转录本和基因的表达水平进行定量。采用FPKM(fragments per kilobase of transcript per million fragments mapped)作为衡量转录本或基因表达水平的指标(Floreaetal., 2013)。

1.5 精巢和卵巢差异表达基因的筛选

采用DESeq进行样品组间的差异表达基因分析(Anders and Huber, 2010)。以相对差异倍数≥2且FDR<0.05作为筛选标准。在进行差异表达分析过程中，采用Benjamini-Hochberg校正方法对原有假设检验得到的显著性P值(P-value)进行校正，并最终采用FDR作为差异表达基因筛选的关键指标(Westfall, 2008)。将所有的差异表达基因映射到GO数据库中进行GO富集分析。并在KEGG数据库中进行通路富集分析。

2 结果

2.1 精巢和卵巢转录组测序

本试验共得到75 304 191条高质量读段，18.95 Gb的过滤后数据，各组样品序列累计长度介于4 478 945 552～4 893 890 242 bp之间，均达到4.47 Gb以上，GC含量在39.25%～40.07%之间，属于正常范围(表1)。Q30碱基百分比均不小于94.37%，表明测序的碱基识别准确度很高。 与中华蜜蜂基因组序列进行比对，每个样品比对到参考基因组唯一位置的读段数介于32 601 194～35 445 368之间，其比对比率介于88.38%～91.60%之间。实验原始数据已上传至NCBI序列数据库中，访问号为SRR15276363-SRR15276366。

表1 中华蜜蜂精巢和卵巢转录组测序数据统计Table 1 Summary of the transcriptome sequencing data of the testis and ovary of Apis cerana cerana

2.2 生物学重复相关性评估

将皮尔逊相关系数(Pearson’s correlation coefficient)r作为生物学重复相关性的评估指标(Schulzeetal., 2012)。r2越接近1，说明两个重复样品相关性越强。对同一条件的每一对生物学重复样品的基因表达量做相关性热图，如图1所示，处女王卵巢的2个不同生物学重复样品基因表达量之间和成熟雄蜂精巢的2个不同生物学重复样品之间的r2值均大于0.82，而卵巢和精巢间r2值均小于0.35，说明了测序数据的重复性很好。

图1 中华蜜蜂精巢和卵巢转录组中基因表达量的相关性Fig. 1 Correlation of gene expression levels in the transcripome of the testis and ovary of Apis cerana cerana

对筛选出的差异表达基因做层次聚类分析(图2)，发现精巢和卵巢的2个生物学重复分别聚到了一起，而精巢与卵巢间则出现分离，表明本研究所用样品生物学重复性较好，且样本分组较合理。

图2 中华蜜蜂精巢和卵巢转录组中差异表达基因聚类图Fig. 2 Cluster map of differentially expressed genes in the transcripome of the testis and ovary of Apis cerana cerana颜色代表了基因在样品中的表达量水平[log10(FPKM+0.000001)]。The color represents the gene expression level [log10(FPKM+0.000001)] in the sample.

2.3 精巢和卵巢之间差异表达基因

在精巢和卵巢之间共鉴定出5 312个差异表达基因，其中精巢上调表达基因有2 668个，卵巢上调表达基因有2 644个。在精巢中上调倍数最大的基因是hypotheticalproteinAPICC_07601(GenBank登录号: 108003596)，在卵巢中上调倍数最大的基因是myb-likeproteinQ(GenBank登录号: 108001861)。在这些差异表达基因中，通过分析找出了11个与性别决定或精子卵子发生相关的基因(表2)。其中包括2个参与性别决定的基因，6个与卵子形成相关的基因，以及3个参与精子形成的基因。参与性别决定的基因分别是tra2(GenBank登录号: 107997763)与dsx(GenBank登录号: 108001352)。其中tra2在卵巢中表达更高，而另一个与性别决定相关的基因dsx在精巢中表达更高。与卵子形成相关的基因有卵黄蛋白原受体基因(VgR)(GenBank登录号: 107994670),Squid(GenBank登录号: 107993185),Vasa(GenBank登录号: 108004301),Argonaute-3(GenBank登录号: 108004055),Aubergine(GenBank登录号: 108000695)和exuperantia(exu)(GenBank登录号: 108003169)，它们均在卵巢中表达上调。与精子形成相关的基因包括decapentaplegic(Dpp)(GenBank登录号: 108000002),hedgehog(hh)(GenBank登录号: 108003950)和Wnt6(GenBank登录号: 107993075)，它们均在精巢中表达上调。

2.4 差异表达基因的GO和KEGG富集分析

将GO数据库注释的3 154个DEGs按生物学过程、细胞组分和分子功能三大类进行分类，共分为1 458个功能类别。在生物学过程中，跨膜转运(transmembrane transport)显著富集；在细胞组分中，膜的整体构成(integral component of membrane)显著富集；在分子功能类别中，序列特异性DNA结合转录因子活性(transcription factor activity, sequence-specific DNA binding)与转运活性(transporter activity)显著富集(图3)。通过差异表达基因的KEGG通路分析，发现这些基因能归类到132个生化通路，其中，ECM受体交互作用(KO04512)和溶酶体(KO04142)信号通路显著富集(q<0.05)，注释到这两条通路的差异表达基因数目分别为20和37个。

表2 中华蜜蜂精巢和卵巢转录组中与性别决定或精子卵子发生相关的基因Table 2 Genes involved in sex determination, spermatogenesis and oogenesis in the transcriptomeof the testis and ovary of Apis cerana cerana

图3 中华蜜蜂精巢和卵巢转录组中差异表达基因的GO富集Fig. 3 GO enrichment of differentially expressed genes in the transcriptome of the testis and ovary of Apis cerana cerana

对在精巢上调表达的基因进行GO(图4: A)和KEGG(图5: A)分析，发现有10个GO条目显著富集(q<0.05)，分别是生物学过程类别中的跨膜转运、氧化还原过程(oxidation-reduction process)、蛋白质水解(proteolysis)、代谢过程(metabolic process)、ATP水解偶联质子运输(ATP hydrolysis coupled proton transport)；分子功能类别中的转运活性、跨膜转运活性(transmembrane transporter activity)、底物特异性跨膜转运活性(substrate-specific transmembrane transporter activity)、血红素结合(heme binding)；细胞组分类别中的膜的整体构成。17个KEGG通路显著富集(q<0.05)，包括溶酶体(lysosome)、丙酮酸代谢(pyruvate metabolism)、色氨酸代谢(tryptophan metabolism)等。

图4 中华蜜蜂精巢(A)和卵巢(B)转录组中上调表达基因的GO显著富集图Fig.4 GO enrichment diagram of up-regulated genes in the transcriptome of the testis (A) and ovary (B) of Apis cerana cerana

对在卵巢上调表达的基因进行GO(图4: B)和KEGG(图5: B)分析，发现有26个GO条目显著富集(q<0.05)，包括生物学过程类别中的DNA模板的转录调控(regulation of transcription, DNA-templated)、细胞分裂(cell division)、DNA复制起始(DNA replication initiation)和DNA复制(DNA replication)等；分子功能类别中的DNA结合(DNA binding)、核苷酸结合(nucleotide binding)和微管结合(microtubule binding)等；细胞组分类别中的细胞内的核糖核蛋白复合体(intracellular ribonucleoprotein complex)、核原生质(nucleoplasm)和核小体(nucleosome)等。10个KEGG通路显著富集(q<0.05)，包括DNA复制(DNA replication)、mRNA监测通路(mRNA surveillance pathway)和碱基切除修复(base excision repair)等。

图5 中华蜜蜂精巢(A)和卵巢(B)转录组中上调表达基因的KEGG显著富集图Fig.5 KEGG enrichment diagram of up-regulated genes in the transcriptome of the testis (A) and ovary (B) of Apis cerana cerana

3 讨论

本研究采用高通量测序技术对中华蜜蜂蜂王卵巢和雄蜂精巢转录组差异进行了分析。发现蜜蜂性别决定基因tra2和dsx在精巢和卵巢之间有表达差异(表2)。在黑腹果蝇Drosophilamelanogaster中，tra2不仅是雌性性别分化所必需的，也是雄性正常精子发生所必需的(Amreinetal., 1988)。其他昆虫中tra2的同源物也已陆续被分离出来，并参与了性别决定，如地中海实蝇Ceratitiscapitata(Paneetal., 2002; Salveminietal., 2009), 家蝇Muscadomestica(Hedigeretal., 2004; Burghardetal., 2005), 油橄榄实蝇Bactroceraoleae(Lagosetal., 2007), 铜绿蝇Luciliacuprina(Concha and Scott, 2009)和按实蝇属Anastrepha(Ruizetal., 2007; Sarnoetal., 2010; Scheteligetal., 2012)。在西方蜜蜂Apismellifera中，Nissen等(2012)推测Am-tra2基因可能具有双重功能，既在雄性中作为一个调控者调控fem雄性拼接，又在雌性中参与增强fem雌性拼接。我们发现tra2在卵巢的表达量(表2)显著高于在精巢的表达量(表2)，表明tra2可能主要在雌蜂中发挥作用。另一个参与性别决定的基因dsx最初的功能被认为可能是决定性腺的性别身份，然后扩展到控制其他组织的性别二态性。在黑腹果蝇中，由于dsx基因的性别特异性选择性剪接，雌性体内会产生雌特异性蛋白DSX-F，而在雄性体内产生的蛋白是DSX-M，前者蛋白启动雌性发育，后者调控雄性发育(Hedigeretal., 2010)。我们发现dsx在精巢中的表达显著上调(表2)，表明它主要在精巢中发挥作用。

分析发现与卵子形成相关的基因VgR,Vasa,Squid等均在卵巢中表达上调(表2)。VgR对昆虫的生殖发育至关重要。它编码一个180～215 kD的卵巢特异性蛋白，VgR参与卵母细胞摄取卵黄原蛋白(Mizutaetal., 2013)。它通常在雌性生殖发育过程中的卵巢高度表达。Boldbaatar等(2008)通过对蜱VgR的RNA干扰试验得出，相比之下，对照组显示出更大的卵子，并且大多数呈现棕色或褐色，表明摄取了卵黄原蛋白，而实验组卵子更小，颜色也非常淡。我们的结果表明，VgR在卵巢中表达上调(表2)，表明该基因参与了卵子的生成过程。Vasa基因最初在果蝇中被鉴定为生殖细胞发育所需的母体效应基因(Schüpbach and Wieschaus, 1986)。此外，vasa对nanosRNA的定位是必需的，也促进了nanosRNA的翻译(Gavisetal., 1996)。VasaRNA在雄性生殖细胞中的表达已在许多生物体内得到证实(Kobayashietal., 2000; Xuetal., 2005; Marraccietal., 2007; Nakkrasae and Damrongpho, 2007; Zhouetal., 2010)。Wang等(2012)利用原位杂交检测拟穴青蟹Scyllaparamamosain在卵发生时Vasa的表达和分布，结果表明，在卵发生过程中，Sp-vasaRNA在Ⅰ期卵浆中、Ⅲ期核区中、Ⅴ期核周区中均有表达。Squid在gurken和oskarmRNA的调控中起着关键作用，是gurkenmRNA向核背侧移动的特异性需要，也在gurkenmRNA的锚定和翻译中发挥作用(Bastock and St Johnston, 2008)。

分析发现与精子形成相关的基因包括hedgehog和Wnt6等均在精巢中表达上调(表2)。Hedgehog信号分子是由信号细胞分泌的一种局域性蛋白配体。在脊椎和无脊椎动物的诸多发育过程中，该通路调控着细胞命运、增殖与分化。Wnt是一种分泌型糖蛋白，它的不同亚型已在多种动物中被找到(Gaoetal., 2014)。Luo等(2015)发现Wnt6在生殖腺的cap细胞中特异性表达。我们的结果表明，Hedgehog和Wnt6在精巢中的表达量高于卵巢中的(表2)，表明这两个基因可能参与了精子的发生。

将注释到GO数据库中的3 154个精巢和卵巢差异表达基因进行GO分类，与发育繁殖和性别相关的生物学过程包括繁殖(GO:0000003)、生长(GO:0040007)和发育过程(GO:0032502)，共涉及58个相关基因(图3)。精巢和卵巢的差异表达基因显著富集于4个功能性类别。其中与性别决定相关的差异表达基因Dsx显著富集于特异性序列DNA结合转录因子活性类别(GO:0003700)，表明这个基因可能作为一种转录因子来调控其他基因的转录拼接。另外，与精子卵子形成相关的差异表达基因VgR,Squid和hedgehog显著富集于膜的整体构成(GO:0016021)(图4: A)，表明这几个基因可能参与了精子与卵子膜的形成过程。为进一步确定差异表达基因参与的主要生化代谢途径和信号通路，我们对差异表达基因进行了KEGG分析，结果显示ECM受体交互作用(KO04512)和溶酶体(KO04142)显著富集(图5： A)。Bello等(2015)指出ECM通过直接或间接的方式影响胚胎和成年生命的所有步骤。具体机理有待于进一步研究论证。