成付平, 袁 芳, 席芳貴, 秦凱鑫, 胡小芬, 王子龍,*
(1. 江西農(nóng)業(yè)大學(xué)蜜蜂研究所, 南昌 330045; 2. 江西省蜜蜂生物學(xué)與飼養(yǎng)重點(diǎn)實(shí)驗(yàn)室, 南昌330045;3. 江西省養(yǎng)蜂研究所, 南昌330052; 4. 江西農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院, 南昌 330045)
中華蜜蜂Apisceranacerana(簡(jiǎn)稱中蜂)是東方蜜蜂Apiscerana的指名亞種,由于它在抵御寒冷、采集零散蜜源、抵抗胡蜂及有節(jié)制地產(chǎn)卵等方面有著顯著的優(yōu)勢(shì)(曾志將, 2020),使它成為我國(guó)飼養(yǎng)的重要蜂種。近些年,隨著研究手段的不斷更新與深入,已由描述性研究轉(zhuǎn)到分子水平上來闡述蜜蜂發(fā)育機(jī)理。健康的蜂群由蜂王、工蜂、雄蜂組成,它們的性別決定是一種特殊的單雙倍體模式,雄性是由未受精卵發(fā)育而來的單倍體,雌性是由受精卵發(fā)育而來的二倍體(Beyeetal., 2003)。在一個(gè)蜂群中,蜂王和雄蜂都擁有發(fā)育完全的生殖腺,但它們的性成熟時(shí)間是不同步的。蜂王在羽化一周后達(dá)到性成熟,而雄蜂則需要大約12 d才能達(dá)到性成熟(曾志將, 2017)。
蜜蜂的性別發(fā)育包括兩個(gè)方面,分別是性別決定和性別分化。在蜜蜂中,性別是由csd→fem→dsx這個(gè)級(jí)聯(lián)反應(yīng)來決定的(Wilkins, 1995; Hasselmannetal., 2008, 2010)。在這個(gè)級(jí)聯(lián)調(diào)控過程中,csd是調(diào)節(jié)fem進(jìn)行雌特異性剪接的初始信號(hào)(Hasselmannetal., 2008, 2010),之后fem又控制著dsx的雌特異性剪接。在多細(xì)胞動(dòng)物中,dsx是性別決定級(jí)聯(lián)反應(yīng)中最保守的性別決定基因之一(Wilkins, 1995)。這一級(jí)聯(lián)反應(yīng)的下游基因充當(dāng)性別調(diào)節(jié)因子,形成性別分化途徑,進(jìn)而導(dǎo)致雌雄特異性表型。
雖然目前已經(jīng)鑒定出蜜蜂性別級(jí)聯(lián)的主效基因,但對(duì)蜜蜂生殖活動(dòng)和性腺發(fā)育的研究較少。因此,有必要了解蜜蜂性腺分化和配子發(fā)生的分子機(jī)制。本研究對(duì)中華蜜蜂的卵巢和精巢進(jìn)行了轉(zhuǎn)錄組測(cè)序分析。我們的結(jié)果揭示了雄性和雌性性腺的基因表達(dá)差異。本研究的結(jié)果有助于理解中華蜜蜂性別決定和分化的分子機(jī)制,并為未來的研究提供有價(jià)值的基因表達(dá)信息。
本研究所用的實(shí)驗(yàn)蜂群為江西省養(yǎng)蜂研究所飼養(yǎng)的中華蜜蜂蜂群。選擇兩個(gè)群勢(shì)較強(qiáng)的中華蜜蜂蜂群,采用標(biāo)準(zhǔn)的人工育王方法(曾志將, 2009)培育中華蜜蜂蜂王,待蜂王出房后隨機(jī)挑選個(gè)體發(fā)育優(yōu)良的處女王剪掉翅后放于巢房中培育,到蜂王出房12 d時(shí)進(jìn)行采樣,每2頭蜂王的卵巢作為一個(gè)樣品,每群采集一個(gè)樣品,共采集2個(gè)生物學(xué)重復(fù)。采用同樣的蜂群,將蜂王控制在一張干凈的雄蜂巢脾中產(chǎn)卵,12 h之后放出蜂王,并將該雄蜂巢脾放入繼箱群中繼續(xù)孵化發(fā)育。等雄蜂羽化出房后進(jìn)行標(biāo)記,再放入原實(shí)驗(yàn)群中進(jìn)一步培育,到雄蜂出房12 d時(shí)進(jìn)行采樣,每10頭雄蜂精巢作為一個(gè)樣品,每群采集一個(gè)樣品,共采集2個(gè)生物學(xué)重復(fù)。取樣時(shí),在顯微鏡下對(duì)雄蜂和蜂王腹部分別進(jìn)行解剖,取出雄蜂精巢及蜂王卵巢,將采集好的精巢和卵巢樣品用液氮速凍后裝在1.5 mL EP管中,保存于液氮中。
按照標(biāo)準(zhǔn)的方法分別提取1.1節(jié)每個(gè)精巢和卵巢樣品的總RNA,對(duì)RNA質(zhì)量進(jìn)行檢測(cè),以確保得到合格樣本并且能夠進(jìn)行轉(zhuǎn)錄組測(cè)序。待實(shí)驗(yàn)樣品總RNA檢測(cè)結(jié)果均合格后,進(jìn)行文庫(kù)構(gòu)建,主要操作流程為:用帶有Oligo(dT)的磁珠富集樣品mRNA;加入Fragmentation Buffer隨機(jī)打斷mRNA片段;用六堿基隨機(jī)引物(random hexamers)合成cDNA第1鏈,然后加入酶Buffer, dNTPs, RNase H和DNA 聚合酶I合成cDNA第2鏈,反轉(zhuǎn)錄結(jié)束后,通過AMPure XP Beads純化cDNA;對(duì)純化過后的cDNA進(jìn)行末端修復(fù)、加poly(A)尾并連接測(cè)序接頭,然后用AMPure XP Beads進(jìn)行片段大小選擇;最后通過PCR擴(kuò)增富集得到cDNA文庫(kù)。用Illumina高通量測(cè)序平臺(tái)對(duì)cDNA文庫(kù)進(jìn)行測(cè)序。
利用FastQC軟件對(duì)數(shù)據(jù)進(jìn)行分析篩選,去除含有接頭的讀段和低質(zhì)量的讀段(包括去除N的比例大于10%的讀段;去除質(zhì)量值Q≤10的堿基數(shù)占整條讀段的50%以上的讀段),最終篩選得到高質(zhì)量的過濾后數(shù)據(jù)。Illumina HiSeq堿基質(zhì)量值(quality score或Q-score)的計(jì)算通常使用Phred堿基質(zhì)量值公式(Ewingetal., 1998):Q-score=-10×log10P,P為堿基識(shí)別出錯(cuò)的概率。
利用TopHat2(Kimetal., 2013)軟件將獲得的過濾后讀段與中華蜜蜂最新的基因組序列(ftp:∥ftp.ncbi.nlm.nih.gov/genomes/genbank/invertebrate/Apis_cerana/latest_assembly_versions/GCA_011100585.1_ASM1110058v1)進(jìn)行序列比對(duì),獲取過濾后讀段在參考基因組和基因上的位置信息,以及測(cè)序樣品特有的序列特征信息。
使用Cufflinks軟件的Cuffquant和Cuffnorm組件,通過匹配上的讀段在基因上的位置信息,對(duì)轉(zhuǎn)錄本和基因的表達(dá)水平進(jìn)行定量。采用FPKM(fragments per kilobase of transcript per million fragments mapped)作為衡量轉(zhuǎn)錄本或基因表達(dá)水平的指標(biāo)(Floreaetal., 2013)。
采用DESeq進(jìn)行樣品組間的差異表達(dá)基因分析(Anders and Huber, 2010)。以相對(duì)差異倍數(shù)≥2且FDR<0.05作為篩選標(biāo)準(zhǔn)。在進(jìn)行差異表達(dá)分析過程中,采用Benjamini-Hochberg校正方法對(duì)原有假設(shè)檢驗(yàn)得到的顯著性P值(P-value)進(jìn)行校正,并最終采用FDR作為差異表達(dá)基因篩選的關(guān)鍵指標(biāo)(Westfall, 2008)。將所有的差異表達(dá)基因映射到GO數(shù)據(jù)庫(kù)中進(jìn)行GO富集分析。并在KEGG數(shù)據(jù)庫(kù)中進(jìn)行通路富集分析。
本試驗(yàn)共得到75 304 191條高質(zhì)量讀段,18.95 Gb的過濾后數(shù)據(jù),各組樣品序列累計(jì)長(zhǎng)度介于4 478 945 552~4 893 890 242 bp之間,均達(dá)到4.47 Gb以上,GC含量在39.25%~40.07%之間,屬于正常范圍(表1)。Q30堿基百分比均不小于94.37%,表明測(cè)序的堿基識(shí)別準(zhǔn)確度很高。 與中華蜜蜂基因組序列進(jìn)行比對(duì),每個(gè)樣品比對(duì)到參考基因組唯一位置的讀段數(shù)介于32 601 194~35 445 368之間,其比對(duì)比率介于88.38%~91.60%之間。實(shí)驗(yàn)原始數(shù)據(jù)已上傳至NCBI序列數(shù)據(jù)庫(kù)中,訪問號(hào)為SRR15276363-SRR15276366。
表1 中華蜜蜂精巢和卵巢轉(zhuǎn)錄組測(cè)序數(shù)據(jù)統(tǒng)計(jì)Table 1 Summary of the transcriptome sequencing data of the testis and ovary of Apis cerana cerana
將皮爾遜相關(guān)系數(shù)(Pearson’s correlation coefficient)r作為生物學(xué)重復(fù)相關(guān)性的評(píng)估指標(biāo)(Schulzeetal., 2012)。r2越接近1,說明兩個(gè)重復(fù)樣品相關(guān)性越強(qiáng)。對(duì)同一條件的每一對(duì)生物學(xué)重復(fù)樣品的基因表達(dá)量做相關(guān)性熱圖,如圖1所示,處女王卵巢的2個(gè)不同生物學(xué)重復(fù)樣品基因表達(dá)量之間和成熟雄蜂精巢的2個(gè)不同生物學(xué)重復(fù)樣品之間的r2值均大于0.82,而卵巢和精巢間r2值均小于0.35,說明了測(cè)序數(shù)據(jù)的重復(fù)性很好。
圖1 中華蜜蜂精巢和卵巢轉(zhuǎn)錄組中基因表達(dá)量的相關(guān)性Fig. 1 Correlation of gene expression levels in the transcripome of the testis and ovary of Apis cerana cerana
對(duì)篩選出的差異表達(dá)基因做層次聚類分析(圖2),發(fā)現(xiàn)精巢和卵巢的2個(gè)生物學(xué)重復(fù)分別聚到了一起,而精巢與卵巢間則出現(xiàn)分離,表明本研究所用樣品生物學(xué)重復(fù)性較好,且樣本分組較合理。
圖2 中華蜜蜂精巢和卵巢轉(zhuǎn)錄組中差異表達(dá)基因聚類圖Fig. 2 Cluster map of differentially expressed genes in the transcripome of the testis and ovary of Apis cerana cerana顏色代表了基因在樣品中的表達(dá)量水平[log10(FPKM+0.000001)]。The color represents the gene expression level [log10(FPKM+0.000001)] in the sample.
在精巢和卵巢之間共鑒定出5 312個(gè)差異表達(dá)基因,其中精巢上調(diào)表達(dá)基因有2 668個(gè),卵巢上調(diào)表達(dá)基因有2 644個(gè)。在精巢中上調(diào)倍數(shù)最大的基因是hypotheticalproteinAPICC_07601(GenBank登錄號(hào): 108003596),在卵巢中上調(diào)倍數(shù)最大的基因是myb-likeproteinQ(GenBank登錄號(hào): 108001861)。在這些差異表達(dá)基因中,通過分析找出了11個(gè)與性別決定或精子卵子發(fā)生相關(guān)的基因(表2)。其中包括2個(gè)參與性別決定的基因,6個(gè)與卵子形成相關(guān)的基因,以及3個(gè)參與精子形成的基因。參與性別決定的基因分別是tra2(GenBank登錄號(hào): 107997763)與dsx(GenBank登錄號(hào): 108001352)。其中tra2在卵巢中表達(dá)更高,而另一個(gè)與性別決定相關(guān)的基因dsx在精巢中表達(dá)更高。與卵子形成相關(guān)的基因有卵黃蛋白原受體基因(VgR)(GenBank登錄號(hào): 107994670),Squid(GenBank登錄號(hào): 107993185),Vasa(GenBank登錄號(hào): 108004301),Argonaute-3(GenBank登錄號(hào): 108004055),Aubergine(GenBank登錄號(hào): 108000695)和exuperantia(exu)(GenBank登錄號(hào): 108003169),它們均在卵巢中表達(dá)上調(diào)。與精子形成相關(guān)的基因包括decapentaplegic(Dpp)(GenBank登錄號(hào): 108000002),hedgehog(hh)(GenBank登錄號(hào): 108003950)和Wnt6(GenBank登錄號(hào): 107993075),它們均在精巢中表達(dá)上調(diào)。
將GO數(shù)據(jù)庫(kù)注釋的3 154個(gè)DEGs按生物學(xué)過程、細(xì)胞組分和分子功能三大類進(jìn)行分類,共分為1 458個(gè)功能類別。在生物學(xué)過程中,跨膜轉(zhuǎn)運(yùn)(transmembrane transport)顯著富集;在細(xì)胞組分中,膜的整體構(gòu)成(integral component of membrane)顯著富集;在分子功能類別中,序列特異性DNA結(jié)合轉(zhuǎn)錄因子活性(transcription factor activity, sequence-specific DNA binding)與轉(zhuǎn)運(yùn)活性(transporter activity)顯著富集(圖3)。通過差異表達(dá)基因的KEGG通路分析,發(fā)現(xiàn)這些基因能歸類到132個(gè)生化通路,其中,ECM受體交互作用(KO04512)和溶酶體(KO04142)信號(hào)通路顯著富集(q<0.05),注釋到這兩條通路的差異表達(dá)基因數(shù)目分別為20和37個(gè)。
表2 中華蜜蜂精巢和卵巢轉(zhuǎn)錄組中與性別決定或精子卵子發(fā)生相關(guān)的基因Table 2 Genes involved in sex determination, spermatogenesis and oogenesis in the transcriptomeof the testis and ovary of Apis cerana cerana
圖3 中華蜜蜂精巢和卵巢轉(zhuǎn)錄組中差異表達(dá)基因的GO富集Fig. 3 GO enrichment of differentially expressed genes in the transcriptome of the testis and ovary of Apis cerana cerana
對(duì)在精巢上調(diào)表達(dá)的基因進(jìn)行GO(圖4: A)和KEGG(圖5: A)分析,發(fā)現(xiàn)有10個(gè)GO條目顯著富集(q<0.05),分別是生物學(xué)過程類別中的跨膜轉(zhuǎn)運(yùn)、氧化還原過程(oxidation-reduction process)、蛋白質(zhì)水解(proteolysis)、代謝過程(metabolic process)、ATP水解偶聯(lián)質(zhì)子運(yùn)輸(ATP hydrolysis coupled proton transport);分子功能類別中的轉(zhuǎn)運(yùn)活性、跨膜轉(zhuǎn)運(yùn)活性(transmembrane transporter activity)、底物特異性跨膜轉(zhuǎn)運(yùn)活性(substrate-specific transmembrane transporter activity)、血紅素結(jié)合(heme binding);細(xì)胞組分類別中的膜的整體構(gòu)成。17個(gè)KEGG通路顯著富集(q<0.05),包括溶酶體(lysosome)、丙酮酸代謝(pyruvate metabolism)、色氨酸代謝(tryptophan metabolism)等。
圖4 中華蜜蜂精巢(A)和卵巢(B)轉(zhuǎn)錄組中上調(diào)表達(dá)基因的GO顯著富集圖Fig.4 GO enrichment diagram of up-regulated genes in the transcriptome of the testis (A) and ovary (B) of Apis cerana cerana
對(duì)在卵巢上調(diào)表達(dá)的基因進(jìn)行GO(圖4: B)和KEGG(圖5: B)分析,發(fā)現(xiàn)有26個(gè)GO條目顯著富集(q<0.05),包括生物學(xué)過程類別中的DNA模板的轉(zhuǎn)錄調(diào)控(regulation of transcription, DNA-templated)、細(xì)胞分裂(cell division)、DNA復(fù)制起始(DNA replication initiation)和DNA復(fù)制(DNA replication)等;分子功能類別中的DNA結(jié)合(DNA binding)、核苷酸結(jié)合(nucleotide binding)和微管結(jié)合(microtubule binding)等;細(xì)胞組分類別中的細(xì)胞內(nèi)的核糖核蛋白復(fù)合體(intracellular ribonucleoprotein complex)、核原生質(zhì)(nucleoplasm)和核小體(nucleosome)等。10個(gè)KEGG通路顯著富集(q<0.05),包括DNA復(fù)制(DNA replication)、mRNA監(jiān)測(cè)通路(mRNA surveillance pathway)和堿基切除修復(fù)(base excision repair)等。
圖5 中華蜜蜂精巢(A)和卵巢(B)轉(zhuǎn)錄組中上調(diào)表達(dá)基因的KEGG顯著富集圖Fig.5 KEGG enrichment diagram of up-regulated genes in the transcriptome of the testis (A) and ovary (B) of Apis cerana cerana
本研究采用高通量測(cè)序技術(shù)對(duì)中華蜜蜂蜂王卵巢和雄蜂精巢轉(zhuǎn)錄組差異進(jìn)行了分析。發(fā)現(xiàn)蜜蜂性別決定基因tra2和dsx在精巢和卵巢之間有表達(dá)差異(表2)。在黑腹果蠅Drosophilamelanogaster中,tra2不僅是雌性性別分化所必需的,也是雄性正常精子發(fā)生所必需的(Amreinetal., 1988)。其他昆蟲中tra2的同源物也已陸續(xù)被分離出來,并參與了性別決定,如地中海實(shí)蠅Ceratitiscapitata(Paneetal., 2002; Salveminietal., 2009), 家蠅Muscadomestica(Hedigeretal., 2004; Burghardetal., 2005), 油橄欖實(shí)蠅Bactroceraoleae(Lagosetal., 2007), 銅綠蠅Luciliacuprina(Concha and Scott, 2009)和按實(shí)蠅屬Anastrepha(Ruizetal., 2007; Sarnoetal., 2010; Scheteligetal., 2012)。在西方蜜蜂Apismellifera中,Nissen等(2012)推測(cè)Am-tra2基因可能具有雙重功能,既在雄性中作為一個(gè)調(diào)控者調(diào)控fem雄性拼接,又在雌性中參與增強(qiáng)fem雌性拼接。我們發(fā)現(xiàn)tra2在卵巢的表達(dá)量(表2)顯著高于在精巢的表達(dá)量(表2),表明tra2可能主要在雌蜂中發(fā)揮作用。另一個(gè)參與性別決定的基因dsx最初的功能被認(rèn)為可能是決定性腺的性別身份,然后擴(kuò)展到控制其他組織的性別二態(tài)性。在黑腹果蠅中,由于dsx基因的性別特異性選擇性剪接,雌性體內(nèi)會(huì)產(chǎn)生雌特異性蛋白DSX-F,而在雄性體內(nèi)產(chǎn)生的蛋白是DSX-M,前者蛋白啟動(dòng)雌性發(fā)育,后者調(diào)控雄性發(fā)育(Hedigeretal., 2010)。我們發(fā)現(xiàn)dsx在精巢中的表達(dá)顯著上調(diào)(表2),表明它主要在精巢中發(fā)揮作用。
分析發(fā)現(xiàn)與卵子形成相關(guān)的基因VgR,Vasa,Squid等均在卵巢中表達(dá)上調(diào)(表2)。VgR對(duì)昆蟲的生殖發(fā)育至關(guān)重要。它編碼一個(gè)180~215 kD的卵巢特異性蛋白,VgR參與卵母細(xì)胞攝取卵黃原蛋白(Mizutaetal., 2013)。它通常在雌性生殖發(fā)育過程中的卵巢高度表達(dá)。Boldbaatar等(2008)通過對(duì)蜱VgR的RNA干擾試驗(yàn)得出,相比之下,對(duì)照組顯示出更大的卵子,并且大多數(shù)呈現(xiàn)棕色或褐色,表明攝取了卵黃原蛋白,而實(shí)驗(yàn)組卵子更小,顏色也非常淡。我們的結(jié)果表明,VgR在卵巢中表達(dá)上調(diào)(表2),表明該基因參與了卵子的生成過程。Vasa基因最初在果蠅中被鑒定為生殖細(xì)胞發(fā)育所需的母體效應(yīng)基因(Schüpbach and Wieschaus, 1986)。此外,vasa對(duì)nanosRNA的定位是必需的,也促進(jìn)了nanosRNA的翻譯(Gavisetal., 1996)。VasaRNA在雄性生殖細(xì)胞中的表達(dá)已在許多生物體內(nèi)得到證實(shí)(Kobayashietal., 2000; Xuetal., 2005; Marraccietal., 2007; Nakkrasae and Damrongpho, 2007; Zhouetal., 2010)。Wang等(2012)利用原位雜交檢測(cè)擬穴青蟹Scyllaparamamosain在卵發(fā)生時(shí)Vasa的表達(dá)和分布,結(jié)果表明,在卵發(fā)生過程中,Sp-vasaRNA在Ⅰ期卵漿中、Ⅲ期核區(qū)中、Ⅴ期核周區(qū)中均有表達(dá)。Squid在gurken和oskarmRNA的調(diào)控中起著關(guān)鍵作用,是gurkenmRNA向核背側(cè)移動(dòng)的特異性需要,也在gurkenmRNA的錨定和翻譯中發(fā)揮作用(Bastock and St Johnston, 2008)。
分析發(fā)現(xiàn)與精子形成相關(guān)的基因包括hedgehog和Wnt6等均在精巢中表達(dá)上調(diào)(表2)。Hedgehog信號(hào)分子是由信號(hào)細(xì)胞分泌的一種局域性蛋白配體。在脊椎和無脊椎動(dòng)物的諸多發(fā)育過程中,該通路調(diào)控著細(xì)胞命運(yùn)、增殖與分化。Wnt是一種分泌型糖蛋白,它的不同亞型已在多種動(dòng)物中被找到(Gaoetal., 2014)。Luo等(2015)發(fā)現(xiàn)Wnt6在生殖腺的cap細(xì)胞中特異性表達(dá)。我們的結(jié)果表明,Hedgehog和Wnt6在精巢中的表達(dá)量高于卵巢中的(表2),表明這兩個(gè)基因可能參與了精子的發(fā)生。
將注釋到GO數(shù)據(jù)庫(kù)中的3 154個(gè)精巢和卵巢差異表達(dá)基因進(jìn)行GO分類,與發(fā)育繁殖和性別相關(guān)的生物學(xué)過程包括繁殖(GO:0000003)、生長(zhǎng)(GO:0040007)和發(fā)育過程(GO:0032502),共涉及58個(gè)相關(guān)基因(圖3)。精巢和卵巢的差異表達(dá)基因顯著富集于4個(gè)功能性類別。其中與性別決定相關(guān)的差異表達(dá)基因Dsx顯著富集于特異性序列DNA結(jié)合轉(zhuǎn)錄因子活性類別(GO:0003700),表明這個(gè)基因可能作為一種轉(zhuǎn)錄因子來調(diào)控其他基因的轉(zhuǎn)錄拼接。另外,與精子卵子形成相關(guān)的差異表達(dá)基因VgR,Squid和hedgehog顯著富集于膜的整體構(gòu)成(GO:0016021)(圖4: A),表明這幾個(gè)基因可能參與了精子與卵子膜的形成過程。為進(jìn)一步確定差異表達(dá)基因參與的主要生化代謝途徑和信號(hào)通路,我們對(duì)差異表達(dá)基因進(jìn)行了KEGG分析,結(jié)果顯示ECM受體交互作用(KO04512)和溶酶體(KO04142)顯著富集(圖5: A)。Bello等(2015)指出ECM通過直接或間接的方式影響胚胎和成年生命的所有步驟。具體機(jī)理有待于進(jìn)一步研究論證。