中華蜜蜂精巢和卵巢差異表達基因分析

成付平, 袁 芳, 席芳貴, 秦凱鑫, 胡小芬, 王子龍,*

(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通過直接或間接的方式影響胚胎和成年生命的所有步驟。具體機理有待于進一步研究論證。