半番鴨與番鴨精巢組織差異表達(dá)轉(zhuǎn)錄組測(cè)序分析

李麗, 繆中緯，辛清武，朱志明, 章琳俐，莊曉東，鄭嫩珠,3

（1福建省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所，福州 350013；2福建昌龍集團(tuán)，福建漳州 363000;3福建農(nóng)林大學(xué)食品科學(xué)學(xué)院,福州351100）

畜牧·獸醫(yī)·資源昆蟲(chóng)

半番鴨與番鴨精巢組織差異表達(dá)轉(zhuǎn)錄組測(cè)序分析

李麗1, 繆中緯1，辛清武1，朱志明1, 章琳俐1，莊曉東2，鄭嫩珠1,3

（1福建省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所，福州 350013；2福建昌龍集團(tuán)，福建漳州 363000;3福建農(nóng)林大學(xué)食品科學(xué)學(xué)院,福州351100）

【目的】研究半番鴨與番鴨精巢組織轉(zhuǎn)錄組差異表達(dá)基因，為進(jìn)一步闡明半番鴨不育的遺傳機(jī)制奠定理論基礎(chǔ)?！痉椒ā坷棉D(zhuǎn)錄組測(cè)序方法對(duì)半番鴨和番鴨的精巢組織進(jìn)行研究，篩選其差異表達(dá)基因，并對(duì)其功能進(jìn)行注釋和熒光定量PCR(quantitative real-time PCR, QRT-PCR)驗(yàn)證?！窘Y(jié)果】測(cè)序共獲得43.84Gb Clean Data，組裝后共獲得193 535條Unigene。DESeq分析發(fā)現(xiàn)3 597個(gè)基因在兩個(gè)鴨品種間差異表達(dá)，其中上調(diào)基因1 194個(gè)和下調(diào)基因2 403個(gè)，包括與生殖功能相關(guān)的基因，如成纖維細(xì)胞生長(zhǎng)因子（fibroblast growth factor, FGF）、蛋白激酶A（protein kinaseA, PKA）、絲裂原活化蛋白激酶7（mitogen-activated protein kinase 7-like, partial, BMK）、生長(zhǎng)因子受體結(jié)合蛋白2（growth factor receptor-bound protein 2, GRB2）、c-Jun氨基末端激酶(c-Jun N-terminal kinase, JNK)和腫瘤壞死因子受體超家族成員6（tumor necrosis factor receptor superfamily member 6, FAS）等。GO (gene Ontology)分析發(fā)現(xiàn)382個(gè)差異基因獲得功能注釋，其中97個(gè)基因涉及發(fā)育繁殖生物學(xué)過(guò)程。KEGG (kyoto Encyclopedia of Genes and Genomes)通路分析表明差異表達(dá)基因共富集到 50 條信號(hào)通路中，其中 17個(gè)通路顯著富集，包括絲裂原活化蛋白激酶信號(hào)轉(zhuǎn)導(dǎo)通路(mitogen-activated protein kinase signaling pathway, MAPK)、甘油酯代謝（glycerolipid metabolism）以及鈣信號(hào)途徑 (calcium signaling pathway)信號(hào)通路等，與生殖功能密切相關(guān)的有促性腺激素釋放激素信號(hào)通路（gonadotropin releasing hormone (GnRH) signaling pathway）和MAPK信號(hào)轉(zhuǎn)導(dǎo)通路。經(jīng)QRT-PCR驗(yàn)證，差異基因表達(dá)變化模式與轉(zhuǎn)錄組測(cè)序結(jié)果一致，測(cè)序結(jié)果可靠?！窘Y(jié)論】在轉(zhuǎn)錄組水平上篩選出半番鴨和番鴨精巢組織差異表達(dá)基因，揭示了GnRH和MAPK信號(hào)通路在鴨的生殖活動(dòng)中發(fā)揮了重要作用，為進(jìn)一步探索半番鴨生殖系統(tǒng)的分化機(jī)理提供可靠的參考依據(jù)。

半番鴨；番鴨；精巢組織；差異表達(dá)基因；轉(zhuǎn)錄組

Abstract:【Objective】 The purpose of this study is to analyze the transcriptome differential gene expression of mule duck and muscovy duck testis, results of the study will lay a theoretical foundation for the further elucidation of the mechanisms of genetic sterility of mule duck. 【Method】 Transcriptome sequencing of testis from mule duck and muscovy duck was performed using the Illumina HiSeq 2500 platform with 2 biological replicates per duck breed, and verified by quantitative real-time PCR(QRT-PCR). 【Result】 After removing sequencing adaptors and the low-quality reads, a total of 43.84 Gb clean reads wereobtained, the Q30 base percentage at 91.36% and above, and clean reads were assembled into 193 535 Unigene. The r2Differential expression analysis showed that 3 597 differentially expressed genes were found between two duck breeds, including 1 194 up regulated genes and 2 403 down-regulated genes. Several genes were related with reproductive function, such as fibroblast growth factor, FGF, protein kinase A (PKA), mitogen-activated protein kinase 7-like, partial (BMK),growth factor receptor-bound protein 2 (GRB2), c-Jun N-terminal kinase (JNK )and tumor necrosis factor receptor superfamily member 6(FAS),and so on. With Gene Ontology (GO) annotation, 382 differentially expressed genes were identified including 97 related annotation genes involving development and reproduction biological process. Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis showed that differentially expression genes annotated to 50 metabolic pathway, and 17 pathways were enriched significantly, such as calcium signaling pathway, glycerolipid metabolism and MAPK signaling pathway, and gonadotropin releasing hormone (GnRH)signaling pathway and MAPK signaling pathway associated with testis physiology and reproduction activities. Verified by QRT-PCR, the pattern of differential gene expression was consistent with the results of transcriptome sequencing, which showed the sequencing results were reliable. 【Conclusion】 Differentially expressed genes of mule duck and muscovy duck testis were screened by transcriptional analysis, revealed that the GnRH and MAPK signaling pathway play an important role in duck reproductive activities, which provide reliable reference for exploring the differentiation mechanism of reproductive system in mule duck.

Key words:mule duck; muscovy duck; testis tissue; differentially expressed genes; transcriptome

0 引言

【研究意義】番鴨屬的公番鴨與河鴨屬的母家鴨親緣關(guān)系遠(yuǎn)，其屬間雜交后代為無(wú)繁殖能力的騾鴨，俗稱半番鴨。半番鴨抗逆性強(qiáng)，飼料報(bào)酬高，肉質(zhì)細(xì)膩，具有很強(qiáng)的雜種優(yōu)勢(shì)，但其公母外型差異不顯著，性器官、性行為分化不明顯，無(wú)法正常繁育后代，無(wú)實(shí)際的種用價(jià)值。實(shí)踐研究發(fā)現(xiàn)，半番鴨并非完全意義上的不育，它們擁有精巢組織，極少數(shù)可產(chǎn)生少量精液，少數(shù)可正常交配，雄性偶見(jiàn)爬跨，雌性偶見(jiàn)產(chǎn)蛋，這可能涉及到基因表達(dá)、調(diào)控以及性細(xì)胞形成方式等。【前人研究進(jìn)展】遠(yuǎn)緣雜交不育是一個(gè)很復(fù)雜的生物學(xué)現(xiàn)象，其遺傳基礎(chǔ)較復(fù)雜，可能受眾多的基因及其產(chǎn)物的調(diào)控。關(guān)于鴨類屬間雜交不育的原因，檀俊秩和陳暉等[1]認(rèn)為半番鴨父母本的染色體核型不一致，減數(shù)分裂受到破壞，進(jìn)而影響性腺軸發(fā)育，最終導(dǎo)致其后代不育。劉軍須等[2]研究認(rèn)為大鼠不育性狀由常染色體上單一隱性基因控制，呈隱性遺傳。張慶波等[3]認(rèn)為 DAZL（Deleted in AZoospermia-like）基因是牛精子發(fā)生的重要調(diào)控因子，其突變或表達(dá)缺乏將導(dǎo)致雄性不育。以 Illumina為基礎(chǔ)的轉(zhuǎn)錄組測(cè)序（RNA-Seq）技術(shù)被廣泛的應(yīng)用到差異基因的檢測(cè)及功能注釋，在畜禽方面也應(yīng)用廣泛，如鑒定鴨黑、白羽相關(guān)基因[4]，篩選野雞和家雞肉質(zhì)基因[5]，分析紹興鴨青殼性狀相關(guān)基因[6]等研究。生殖、發(fā)育機(jī)制等一直是生物學(xué)熱點(diǎn)，BAUERSACHS等[7]通過(guò)RNA-Seq技術(shù)篩選到不同物種間繁殖相關(guān)差異基因，F(xiàn)IEDLER等[8]對(duì)海兔不同發(fā)育階段進(jìn)行轉(zhuǎn)錄組測(cè)序，張偉[9]構(gòu)建了中華絨螯蟹精巢組織文庫(kù)，篩選雄性生殖調(diào)控關(guān)鍵基因，鐘志君[10]研究比較了成年藏豬的精巢和卵巢轉(zhuǎn)錄組表達(dá)譜，張升利等[11]分析了長(zhǎng)尾草金魚成熟期精巢和卵巢轉(zhuǎn)錄組差異表達(dá)基因，XU等[12]利用轉(zhuǎn)錄組測(cè)序方法挖掘出與卵泡發(fā)育相關(guān)基因，朱志明等[13]探明了山麻鴨開(kāi)產(chǎn)期和產(chǎn)蛋高峰期卵巢組織的轉(zhuǎn)錄組差異。這些成果為研究鴨不育性狀和繁殖性能的遺傳機(jī)制奠定了基礎(chǔ)。【本研究切入點(diǎn)】雖然 RNA-Seq技術(shù)已被應(yīng)用到鴨的性狀研究，但由于半番鴨不育性狀的特殊性，其與番鴨或家鴨精巢組織的轉(zhuǎn)錄組比較分析尚未有報(bào)道，其不育遺傳基礎(chǔ)需進(jìn)一步深入研究?！緮M解決的關(guān)鍵問(wèn)題】本研究通過(guò)對(duì)半番鴨和番鴨個(gè)體精巢組織的轉(zhuǎn)錄組比較，對(duì)其差異表達(dá)基因進(jìn)行篩選，并進(jìn)行GO（gene ontology）與 KEGG（kyoto encyclopediaof genes and genomes）功能富集，分析通路功能并探索精巢分化相關(guān)的差異表達(dá)基因，通過(guò)轉(zhuǎn)錄組數(shù)據(jù)挖掘半番鴨雄性不育相關(guān)的基因，為后續(xù)半番鴨雄性不育形成的遺傳機(jī)制研究提供線索。

1 材料與方法

1.1 試驗(yàn)材料

本研究所用的半番鴨和番鴨公鴨由福建省農(nóng)業(yè)科學(xué)院畜牧獸醫(yī)研究所動(dòng)物房提供，試驗(yàn)于3—6月進(jìn)行，各試驗(yàn)鴨飼養(yǎng)管理?xiàng)l件一致。180日齡性成熟時(shí)，對(duì)公鴨進(jìn)行采精訓(xùn)練。210日齡禁飼12 h后，每個(gè)品種分別取2 只個(gè)體用于轉(zhuǎn)錄組測(cè)序，其中半番鴨公鴨為無(wú)爬跨行為個(gè)體，番鴨公鴨為正常個(gè)體。按照國(guó)家實(shí)驗(yàn)動(dòng)物處理行為準(zhǔn)則屠宰，取精巢組織，置于-80℃?zhèn)溆谩?/p>

1.2 RNA 提取及轉(zhuǎn)錄組測(cè)序

利用RNA easy Lipid Tissue Mini Kit（QIAGEN,Germany）提取每個(gè)個(gè)體總 RNA，單個(gè)建池。采用Nanodrop、Qubit 2.0和Aglient 2100方法檢測(cè)各RNA樣品的純度、濃度及完整性等。構(gòu)建文庫(kù)，Qubit2.0和 Agilent 2100分別對(duì)文庫(kù)的濃度和插入片段大?。↖nsert Size）進(jìn)行檢測(cè)，QRT-PCR對(duì)文庫(kù)的有效濃度進(jìn)行準(zhǔn)確定量，以保證文庫(kù)質(zhì)量。基于邊合成邊測(cè)序（Sequencing By Synthesis，SBS）技術(shù)，利用 Illumina HiSeq 2500（Illumina, America）平臺(tái)進(jìn)行高通量測(cè)序,測(cè)序讀長(zhǎng)為PE150。

1.3 測(cè)序數(shù)據(jù)處理

測(cè)序數(shù)據(jù)去除接頭及低質(zhì)量數(shù)據(jù)后，利用Trinity軟件將數(shù)據(jù)組裝成轉(zhuǎn)錄本，進(jìn)行轉(zhuǎn)錄注釋及表達(dá)量的計(jì)算。利用DESeq進(jìn)行基因的差異表達(dá)分析，繪制差異表達(dá)基因火山圖，并進(jìn)行聚類分析。

1.4 GO注釋、KEGG通路分析

利用GO數(shù)據(jù)庫(kù)對(duì)差異表達(dá)基因進(jìn)行功能注釋，采用COG（cluster of orthologous groups of proteins）對(duì)差異表達(dá)基因進(jìn)行分類統(tǒng)計(jì)，運(yùn)用KEGG數(shù)據(jù)庫(kù)進(jìn)行通路分析，均以P＜0.05作為顯著性富集標(biāo)準(zhǔn)。

1.5 Real-time PCR檢測(cè)

采用Primer Premier 6.0和Beacon designer 7.8軟件設(shè)計(jì)熒光定量PCR引物，然后由生物工程（上海）股份有限公司負(fù)責(zé)合成，引物序列如表1，QRT-PCR擴(kuò)增體系和反應(yīng)條件如表2。以甘油醛-3-磷酸脫氫酶基因（glyceraldehyde-3-phosphate dehydrogenase, GAPDH）做內(nèi)參，采用2-ΔΔCt法計(jì)算基因的相對(duì)表達(dá)量。以P＜0.05作為顯著性標(biāo)準(zhǔn)。

表1 熒光定量引物序列Table 1 Real-time PCR primers and conditions

表2 熒光定量PCR反應(yīng)體系及條件Table 2 Real-time PCR reaction system and conditions

2 結(jié)果

2.1 半番鴨與番鴨精巢組織轉(zhuǎn)錄組測(cè)序

本試驗(yàn)共構(gòu)建半番鴨和番鴨4個(gè)精巢組織的轉(zhuǎn)錄組文庫(kù)。測(cè)序數(shù)據(jù)去除接頭以及低質(zhì)量數(shù)據(jù)后共獲得288 008 582個(gè)高質(zhì)量數(shù)據(jù)。從表3中可以看出，本試驗(yàn)共獲得43.84Gb的Clean Data，各樣品Clean Data均達(dá)到6.29Gb。另外，各樣品GC含量均不小于51.88%，Q30（Clean Data質(zhì)量不小于30的堿基所占的百分比）全在 91.36%以上，該結(jié)果表明測(cè)序質(zhì)量可靠，構(gòu)建文庫(kù)可用于后續(xù)分析。數(shù)據(jù)組裝后共獲得193 535條Unigene，其中長(zhǎng)度在1kb以上的Unigene有46 175條。

表3 轉(zhuǎn)錄組數(shù)據(jù)組裝情況Table 3 Summary of the sequencing data assembly

圖1 四樣品基因表達(dá)量相關(guān)性圖Fig. 1 Correlation heat map of gene expression level in 4 samples

2.2 差異表達(dá)基因的篩選及聚類分析

FPKM（Fragments Per Kilobase of transcript per Million mapped reads）是每百萬(wàn)Reads中來(lái)自比對(duì)到某一基因每千堿基長(zhǎng)度的Reads數(shù)目，是轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析中常用的基因表達(dá)水平估算方法。將皮爾遜相關(guān)系數(shù)r（Pearson’s Correlation Coefficient）作為生物學(xué)重復(fù)相關(guān)性的評(píng)估指標(biāo)[14]。r2越接近1，說(shuō)明兩個(gè)重復(fù)樣品相關(guān)性越強(qiáng)。對(duì)同一條件的每一對(duì)生物學(xué)重復(fù)樣品的基因表達(dá)量做相關(guān)性圖，相關(guān)性圖見(jiàn)圖 1。半番鴨組內(nèi) 2個(gè)不同重復(fù)樣品和番鴨組內(nèi) 2個(gè)不同重復(fù)樣品的 r2值均大于 0.9，而兩組間 r2值均小于 0.6，說(shuō)明了試驗(yàn)的可靠性和可重復(fù)性都很高。

以FDR（False Discovery Rate）作為差異表達(dá)基因篩選的關(guān)鍵指標(biāo), 將FDR小于0.001且差異倍數(shù)FC（Fold Change）大于等于2作為兩品種鴨個(gè)體間差異表達(dá)基因顯著的篩選標(biāo)準(zhǔn)，共鑒定出3 597個(gè)基因，其中包括1 194個(gè)表達(dá)上調(diào)的差異基因，主要有（FDR由小到大排列）c86758.graph_c1（NADH脫氫酶亞基4），c86758.graph_c0（外周型苯二氮卓受體相關(guān)蛋白1）等差異基因，以及2 403個(gè)表達(dá)下調(diào)基因，包含c243330.graph_c0（DnaJ同源B亞家族成員8），c277017.graph_c0（凋亡因子BCL-2蛋白14）等。進(jìn)一步分析發(fā)現(xiàn)，顯著差異基因中含有一些與繁殖性能相關(guān)的基因（表 4），例如，成纖維細(xì)胞生長(zhǎng)因子（fibroblast growth factor，F(xiàn)GF）、c-Jun氨基末端激酶（c-Jun N-terminal kinase，JNK）、細(xì)胞外信號(hào)調(diào)控的蛋白激酶5（ERK5）、蛋白激酶A（protein kinase A，PKA）、生長(zhǎng)因子受體結(jié)合蛋白 2（growth factor receptor-bound protein 2，Grb2）、絲裂原活化蛋白激酶 7，部分（mitogen-activated protein kinase 7-like,partial，BMK）等下調(diào)基因，腫瘤壞死因子受體超家族成員 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）、胞質(zhì)型磷脂酶 A2（cytosolic phospholipase A2 epsilon-like,CPLA2）等上調(diào)基因。從火山圖（圖 2）能夠快速查看兩組樣品間表達(dá)的差異水表平分布情況。通過(guò) MA圖（圖 3）可以直觀地查看兩組樣品中基因的表達(dá)豐度和差異倍數(shù)的整體分布。

表4 差異表達(dá)基因（部分）Table 4 Differentially expressed genes (parts)

圖2 差異表達(dá)基因火山圖Fig. 2 Volcano plot of differentially expressed genes

圖3 差異表達(dá)基因MA圖Fig. 3 MA plot of differentially expressed genes

對(duì)篩選出來(lái)的差異表達(dá)基因進(jìn)行聚類分析（圖4），發(fā)現(xiàn)同一鴨品種的2 個(gè)生物學(xué)重復(fù)聚到了一起，而不同鴨品種間基因表達(dá)模式則出現(xiàn)分離，表明本研究所用樣本生物學(xué)重復(fù)性較好，且樣本分組較合理。

2.3 差異表達(dá)基因GO功能富集

本研究運(yùn)用GO數(shù)據(jù)庫(kù)對(duì)具有同源比對(duì)的差異基因進(jìn)行生物學(xué)過(guò)程（biological process）、分子功能（molecular function）與細(xì)胞組分（cellular component）三方面的注釋。將比對(duì)得到的1 381個(gè)顯著差異基因進(jìn)行功能注釋，382個(gè)差異基因在GO分類中有功能意義。上述3個(gè)功能被區(qū)分為更具體的61個(gè)類別，分別包括了22、19和20個(gè)功能亞分類。由圖5可知，在生物功能的組分中，兩組間的差異表達(dá)基因在細(xì)胞過(guò)程（cellular process, GO: 0009987）與單一的生物過(guò)程（single-organism process, GO: 0044699）中數(shù)目比例最大。在細(xì)胞功能中，細(xì)胞（cell, GO: 0005623）與細(xì)胞部分（cell part, GO: 0044464）數(shù)目最多。在分子功能分類中，差異基因在結(jié)合（binding,GO: 0005488）中所占的比例最高，催化活性（catalytic activity, GO:0003824）次之。由圖5可知，與發(fā)育繁殖相關(guān)的生物學(xué)過(guò)程有繁殖（reproduction, GO:0000003）,發(fā)育過(guò)程（developmental process, GO:0048589）和繁殖過(guò)程（reproductive process, GO:0022414），涉及97個(gè)相關(guān)基因。

圖4 顯著差異表達(dá)基因聚類分析Fig. 4 Heat map of the differentially expressed genes

2.4 差異表達(dá)基因COG分類

統(tǒng)計(jì)顯著富集的Go term中包含的基因數(shù)，分析結(jié)果見(jiàn)圖6。其中一般的功能預(yù)測(cè)（general function prediction only）基因數(shù)目最多，其次為信號(hào)傳導(dǎo)機(jī)制（signal transduction mechanism），轉(zhuǎn)錄（Transcription）和復(fù)制、重組與修復(fù)（replication, recombination and repair）。

2.5 差異表達(dá)基因KEGG注釋

為確定差異基因參與的主要生化代謝途徑和信號(hào)通路，對(duì)差異表達(dá)基因進(jìn)行KEGG通路分析，結(jié)果顯示差異表達(dá)基因共富集到50 條信號(hào)通路中。以P＜0.05作為差異表達(dá)基因在該通路顯著富集，共鑒定出17個(gè)通路顯著富集：包括與絲裂原活化蛋白激酶（MAPK）信號(hào)轉(zhuǎn)導(dǎo)通路（MAPK signaling pathway）、鈣信號(hào)途徑（calcium signaling pathway）、甘油酯代謝（glycerolipid metabolism）、緊密連接（gap junction）以及血管平滑肌收縮（Vascular smooth muscle contraction）等信號(hào)通路。對(duì)差異表達(dá)基因的注釋結(jié)果按照 KEGG中通路類型進(jìn)行分類，分類圖如圖7所示，其中與生長(zhǎng)發(fā)育、生殖過(guò)程相關(guān)的有 MAPK信號(hào)轉(zhuǎn)導(dǎo)通路和促性腺激素釋放激素信號(hào)通路（gonadotropinreleasing hormone（GnRH）signaling pathway）等。

圖5 差異表達(dá)基因GO注釋Fig. 5 GO annotation of differentially expressed genes

圖6 差異表達(dá)基因COG注釋分類統(tǒng)計(jì)圖Fig. 6 COG annotation classification statistics of differentially expressed genes

圖7 差異表達(dá)基因顯著富集的KEGG 通路Fig. 7 List of KEGG pathway for differentially expressed genes

挑選富集顯著性最可靠（即 Q值最小）的前 20個(gè)通路以散點(diǎn)圖的形式展示（圖 8），在該圖中越靠近右上角的圖形代表的通路，參考價(jià)值越大；反之亦然。由此可知甘油酯代謝（glycerolipid metabolism），鈣信號(hào)途徑（calcium signaling pathway）信號(hào)通路以及血管平滑肌收縮（vascular smooth muscle contraction）富集顯著性較為可靠。

2.6 Real-time PCR驗(yàn)證轉(zhuǎn)錄組測(cè)序數(shù)據(jù)

為驗(yàn)證轉(zhuǎn)錄組測(cè)序結(jié)果，本研究選擇了GRB2、CPLA2、FGF3、DUSP3及FAS基因進(jìn)行Real-time PCR試驗(yàn)。結(jié)果表明除FAS基因整體表達(dá)低，數(shù)據(jù)不準(zhǔn)確外，其他4個(gè)基因在番鴨和半番鴨個(gè)體間表達(dá)變化模式與轉(zhuǎn)錄組測(cè)序結(jié)果一致（表 5），表明本研究利用轉(zhuǎn)錄組測(cè)序獲得的數(shù)據(jù)較為準(zhǔn)確。

3 討論

近來(lái), 高通量轉(zhuǎn)錄組測(cè)序技術(shù)取得了較大的研究進(jìn)展[15-20]。雄性生殖系統(tǒng)的發(fā)育、分化過(guò)程是一個(gè)復(fù)雜的生理過(guò)程，涉及多種代謝過(guò)程，目前已有研究表明p53信號(hào)通路[21]、Wnt代謝通路[22]以及TGFβ通路[23-24]在雄性生殖發(fā)育過(guò)程中起重要作用。為了解鴨精巢的系統(tǒng)機(jī)制和半番鴨不育的特性，本研究對(duì)半番鴨和番鴨精巢組織進(jìn)行轉(zhuǎn)錄組學(xué)測(cè)序及比對(duì)分析，共篩選獲得3 597個(gè)差異表達(dá)基因，同時(shí)為進(jìn)一步確定差異基因參與的主要生化代謝途徑和信號(hào)通路進(jìn)行KEGG分析，結(jié)果發(fā)現(xiàn)MAPK、甘油酯代謝以及鈣信號(hào)途徑等17個(gè)通路富集顯著，其中與生殖過(guò)程相關(guān)的有GnRH信號(hào)通路和MAPK信號(hào)轉(zhuǎn)導(dǎo)通路。

圖8 差異表達(dá)基因KEGG通路富集散點(diǎn)圖Fig. 8 Enriched scatter map of differential expression gene KEGG pathway

表5 Real-time PCR 驗(yàn)證轉(zhuǎn)錄組測(cè)序數(shù)據(jù)Table 5 Validation of the RNA-seq expression data by Real-time PCR for selected genes

MAPK信號(hào)通路是細(xì)胞間信號(hào)傳遞的重要通路[25]，可參與細(xì)胞增殖、分化[26]、精子成熟[27]及凋亡[28-29]等動(dòng)物生殖過(guò)程，被認(rèn)為是精子發(fā)育的重要影響因素之一。本研究篩選到參與此通路的FGF、c- JNK、ERK5等10個(gè)下調(diào)基因，F(xiàn)AS、DUSP3等4個(gè)上調(diào)基因，以及CPLA2、CACN 2個(gè)混合型調(diào)節(jié)基因。在此通路篩選的差異表達(dá)基因中，F(xiàn)GF屬于可促進(jìn)成纖維細(xì)胞生長(zhǎng)的多肽家族。VALVE等[30]研究發(fā)現(xiàn) FGF8基因在成年鼠睪丸精子發(fā)生的特定階段表達(dá)，在鼠、牛卵巢中卵子發(fā)生的特定階段也有表達(dá)。DVORAK等[31]研究表明FGF信號(hào)參與包括細(xì)胞增殖、遷移、分化等多種細(xì)胞應(yīng)答，調(diào)控廣泛的生理和病理過(guò)程。JNK在細(xì)胞分化、凋亡和疾病的發(fā)生均有重要作用，是MAPK信號(hào)轉(zhuǎn)導(dǎo)系統(tǒng)的重要效應(yīng)因子[32]。目前發(fā)現(xiàn)的JNK的底物有轉(zhuǎn)錄因子c-jun、ATF-2及P53等。陳麗莉等[33]發(fā)現(xiàn)JNK可調(diào)控黃鱔卵巢發(fā)育、凋亡及雄性發(fā)育的啟動(dòng)。ERK基因則主要參與動(dòng)物生殖細(xì)胞增殖、分化以及凋亡等多個(gè)發(fā)育過(guò)程[34]。

GnRH信號(hào)通路可調(diào)控動(dòng)物體內(nèi)性腺軸生殖激素的分泌，下丘腦分泌產(chǎn)生的神經(jīng)激素能促進(jìn)垂體分泌促性腺物質(zhì)的釋放，并參與動(dòng)物生殖調(diào)控[35]。本研究篩選到了參與此通路上的PKA、Grb2和BMK等下調(diào)基因，以及CPLA2、CACNA1C 2個(gè)混合調(diào)節(jié)基因。在此通路篩選的差異表達(dá)基因中，已有研究證實(shí)Grb2是信號(hào)轉(zhuǎn)導(dǎo)途徑中的一個(gè)重要成分，它最初是作為表皮生長(zhǎng)因子受體（epidermal growth factor receptor,EGFR）和MAPK通路之間的銜接蛋白被發(fā)現(xiàn)的[36]。Grb2作為一種信號(hào)接頭蛋白，參與細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)過(guò)程，具有促進(jìn)細(xì)胞增殖、細(xì)胞生長(zhǎng)、細(xì)胞分化等功能[37]。PLA2參與雄性生殖過(guò)程, 其與精子獲能、頂體反應(yīng)以及精卵融合等過(guò)程密切相關(guān)，并受到各種信號(hào)通路的調(diào)節(jié)，不同通路相互協(xié)調(diào)。精子PLA2活性的激活及其調(diào)控機(jī)制受G蛋白受體可介導(dǎo)[38]。

另外，JNK、Crb2、CACN等基因同時(shí)參與上述兩個(gè)通路，如PKA、蛋白激酶C（protein kinase C, PKC）及 Ca2+也在生理性頂體反應(yīng)的精子信號(hào)轉(zhuǎn)導(dǎo)中調(diào)節(jié)PLA2 的活性，此外，鈣信號(hào)途徑以及甘油酯代謝、緊密連接以及血管平滑肌收縮等信號(hào)通路也富集顯著，說(shuō)明各種發(fā)揮不同作用的信號(hào)傳導(dǎo)通路交織在一起, 組成一個(gè)復(fù)雜龐大的細(xì)胞信號(hào)傳導(dǎo)網(wǎng)絡(luò)系統(tǒng)。關(guān)于鈣信號(hào)途徑，Ca2+作為細(xì)胞內(nèi)重要第二信使，通過(guò)精子內(nèi)Ca2+濃度調(diào)節(jié)精子生理活動(dòng)。SANTI等[39]也發(fā)現(xiàn) Ca2+信號(hào)通路與精子的形成有著密切的關(guān)系。甘油酯代謝、緊密連接以及血管平滑肌等信號(hào)通路在鴨精巢組織中的作用及其相互調(diào)節(jié)機(jī)制將進(jìn)一步探究。

4 結(jié)論

本研究通過(guò)半番鴨和番鴨轉(zhuǎn)錄組比較，首次在轉(zhuǎn)錄組水平上篩選出與繁殖性能相關(guān)的差異表達(dá)基因，如PKA、Grb2、BMK、CPLA2、FGF、JNK及ERK5等，并進(jìn)一步證實(shí)了GnRH 信號(hào)通路和MAPK信號(hào)通路在公鴨生殖活動(dòng)中發(fā)揮了重要作用。信號(hào)通路的分析為鴨精巢組織的信號(hào)調(diào)控提供了線索，但還需對(duì)涉及這些通路的相關(guān)基因進(jìn)行深入的生物信息學(xué)分析和驗(yàn)證，進(jìn)一步明確繁殖性狀的基因表達(dá)譜和調(diào)控模式，該研究結(jié)果可為今后探索半番鴨生殖系統(tǒng)的分化機(jī)理提供可靠的參考依據(jù)。

[1]檀俊秩, 陳暉, 宋建捷, 劉玉濤. 半番鴨繁殖性狀的研究. 福建農(nóng)業(yè)學(xué)報(bào), 1998, 13(2): 41-45.TAN J Z, CHEN H, SONG J J, LIU Y T. Studies on reproductive character of mule duck. Journal of Fujian Agricultural Sciences, 1998,13(2):41-45. (in Chinese)

[2]劉軍須, 蔡月花, 張敬各, 張福英. 雄性不育大鼠近交系MIJ的建立及其遺傳特征觀察. 動(dòng)物學(xué)雜志, 2008, 43(5) :37-44.LIU J X, CAI Y H, ZHANG J G, ZHANG F Y. Establishment and genetic characteristics of rat inbred strain MIJ with spontaneous male infertility. Journal of Animal Science, 2008, 43(5):37-44. (in Chinese)

[3]張慶波, 李齊發(fā), 李家璜, 李新福, 劉振山, 潘增祥, 宋大偉, 謝莊.牛精子發(fā)生相關(guān)新基因 b-DAZL的克隆、生物信息學(xué)分析與組織表達(dá)研究. 自然科學(xué)進(jìn)展, 2008, 18(5):493-504.ZHANG Q B, LI Q F, LI J H, LI X F, LIU Z S, PAN Z X, SONG D W,XIE Z. Cloning, bioinformatics analysis and tissue expression of a novel b-DAZL gene related to spermatogenesis in cattle. Progress in Natural Science, 2008, 18(5):493-504. (in Chinese)

[4]LI S J, WANG C, YU W H, ZHAO S H, GONG Y Z. Identification of genes related to white and black plumage formation by RNA-Seq from white and black feather bulbs in ducks. PloS One, 2012, 7:e36592.

[5]LI Q H, WANG N, DU Z, HUX X, CHENL, FEI J, WANGY Y, LI.N. Gastrocnemius transcriptome analysis reveals domestication induced gene expression changes between wild and domestic chickens.Genomics, 2012, 100: 314-319.

[6]陳黎, 黃學(xué)濤, 田勇, 陶爭(zhēng)榮, 盧立志. 利用轉(zhuǎn)錄組測(cè)序篩選與鴨青殼性狀形成相關(guān)的基因. 農(nóng)業(yè)生物技術(shù)學(xué)報(bào), 2016, 24(7):1064-1072.CHEN L, HUANG X T, TIAN Y, TAO Z R, LU L Z. Identifying genes associated with blue eggshell in ducks (Anasplatyrhynchos domesticus) by transcriptome analysis. Journal of Agricultural Biotechnology, 2016, 24(7): 1064-1072. (in Chinese)

[7]BAUERSACHS S, WOLF E. Transcriptome analyses of bovine,porcine and equine endometrium during the pre-implantation phase[J].Animal reproduction science. 2012,134:84-94.

[8]FIEDLER T J, HUDDER A, MCKAY S J, SHIVKUMAR S, CAPO T R, SCHMALE M C, WALSH P J. The transcriptome of the early life history stages of the California Sea Hare Aplysia californica,comparative biochemistry and physiology Part D. Genomics &Proteomics, 2010, 5:165-170.

[9]張偉. 中華絨螯蟹精巢組織文庫(kù)構(gòu)建和基因的克隆與序列分析[D].上海: 華東師范大學(xué), 2012.ZHANG W. Construction of cDNA library of Chinese mitten crab,Eriocheir sinensis and cloeing and molecular characterization ofgene[D]. Shanghai：East China Normal University, 2012. (in Chinese)

[10]鐘志君. 豬睪丸和卵巢組織轉(zhuǎn)錄組的差異分析[D]. 雅安：四川農(nóng)業(yè)大學(xué), 2012 .ZHONG Z J. Repertoire of porcine microRNA in adult ovary and testis by deep sequencing[D]. Sichuan: Sichuan Agricultural University, 2012 (in Chinese)

[11]張升利, 付成東, 梁擁軍, 李文通, 孫硯勝, 史東杰, 張 欣. 長(zhǎng)尾草金魚成熟期雌雄性腺 RNA-Seq轉(zhuǎn)錄組分析. 水產(chǎn)科學(xué), 2014,33(12): 750-753.ZHANG S L, FU C D, LIANG Y J, LI W T, SUN Y S, SHI D J,ZHANG X. The RNA-Seq transcriptome analysis in male gonads of Long-tailed goldfish Rassius auratus. Fisheries Science, 2014, 33(12):750-753. (in Chinese)

[12]XU Q, ZHAO W M, CHEN Y, TONG Y Y, RONG G H, HUANG Z Y,ZHANG Y, CHANG G B, WU X S, CHEN G H. Transcriptome profiling of the goose (Anser cygnoides) ovaries identify laying and broodiness phenotypes. PloS One, 2013, 8:e55496.

[13]朱志明, 陳紅萍, 林如龍, 繆中緯, 辛清武, 李麗, 張丹青, 鄭嫩珠.山麻鴨開(kāi)產(chǎn)期和產(chǎn)蛋高峰期卵巢組織轉(zhuǎn)錄組分析. 中國(guó)農(nóng)業(yè)科學(xué),2016, 49(5):998-1007.ZHU Z M, CHEN H P, LIN R L, MIAO Z W, XIN Q W, LI L,ZHANG D Q, ZHENG N Z. Transcriptome analysis of ovary tissue in early laying period and egg laying peak period of shanma Ducks.Scientia Agricultura Sinica, 2016, 49(5):998-1007. (in Chinese)

[14]SCHULZE S K, KANWAR R, G?LZENLEUCHTER M,THERNEAU T M, BEUTLER A S. SERE: Single-parameter quality control and sample comparison for RNA-Seq. [BMC genomics Italic],2012, 13(1): 524.

[15]XIANG L X, HE D, DONG W R, ZHANG Y W, SHAO J Z.Deep sequencing-based transcriptome profiling analysis of bacteria challenged Lateolabrax japonicus reveals insight into the immune-relevant genes in marine fish. BMC Genomics, 2010, 11: 472.

[16]WANG B, GUO G W, WANG C, LIN Y, WANG X M, ZHAO M M,GUO Y, HE M H, ZHANG Y, PAN L. Survey of the transcriptome of Aspergillus oryzae via massively parallel mRNA sequencing. Nucleic Acids Research, 2010, 38(15): 5075-5087.

[17]LI M Y, TAN H W, WANG F, JIANG Q, XU Z S, TIAN C, XIONG A S. De novo transcriptome sequence assembly and identification of AP2/ERF transcription factor related to abiotic stress in parsley(Petroselinum crispum). PLoS One, 2014, 9(9): e108977.

[18]ZHU S Y, TANG S M, TANG Q M, LIU T M. Genome-wide transcriptional changes of ramie (Boehmeria nivea L. Gaud) in response to root-lesion nematode infection. Gene, 2014, 552: 67-74.

[19]CONG F, LIU X L, HAN Z X, SHAO Y H, KONG X P, LIU S W.Transcriptome analysis of chicken kidney tissues following coronavirus avian infectious bronchitis virus infection. Genomics, 2013, 14:743-756.

[20]FRASER B A, WEADICK C J, JANOWITZ I, RODD F H, HUGHES K A. Sequencing and characterization of the guppy (Poecilia reticulata) transcriptome. BMC Genomics, 2011, 12:202.

[21]LI G Y, XIE P, LI H Y, HAO L, XIONG Q, QIU T. Involment of p53,Bax, and Bcl-2 pathway in microcystins-induced apoptosis in rat testis.Environmental Toxicology, 2011, 26(2):111-117.

[22]鄭煒, 趙宗勝, 李青峰, 班謙, 李洪濤, 梁耀偉. Solexa測(cè)序技術(shù)分析雞與鵪鶉屬間雜交雌性和雄性胚胎的差異 microRNAs. 中國(guó)獸醫(yī)學(xué)報(bào), 2014, 34(1):116-122.ZHENG W, ZHAO Z S, LI Q F, BAN Q, LI H T, LIANG H T.Analysis of difference microRNAs in female and male embryos of chicken-quail hybrid by Solexa sequencing. Chinese Journal of Veterinary Science, 2014, 34(1):116-122. (in Chinese)

[23]DRUMMOND A E. TGFβ signaling in the development of ovarian function. Cell and Tissue Research, 2005, 322(l): 107-115.

[24]KONRAD L, KEILANI M M, LAIBLE L, NOTTELMANN U,HOFMANN R. Effects of TGF-betas and a specific antagonist on apoptosis of immature rat male germ cells in vitro. Apoptosis, 2006,11(5): 739-748.

[25]SUN Q Y, BREITBART H, SCHATTEN H. Role of the MAPK cascade in mammalian germ cells. Reproduction, Fertility, and Development, 1999, 11(8):443-450.

[26]HARRIS V K, COTICCHIA C M, KAGAN B L, AHMAD S,WELLSTEIN A, RIEGE A T, HARRIS V K, COTICCHIA C M.Induction of the angiogenic modulator fibroblast growth factorbinding protein by epidermal growth factor is mediated through both MEK/ERK and p38 signal transduction path ways. Journal Biology Chemistry, 2000, 275(15): 10802-10811.

[27]ALMOG T, NAOR Z. Mitogen activated protein kinases (MAPKs) as regulators of spermatogenesis and spermatozoa functions. Molecular and Cellular Endocrinology, 2008, 282:39-44.

[28]JOHNSON C, JIA Y, WANG C, LUE Y H, SWERDLOFF R S,ZHANG X S, HU Z Y, LI Y C, LIU Y X, AMIYA P. Sinha Hikim.Role of Caspase 2 in apoptotic signaling in primate and murine germ cells. Biology of Reproduction, 2008, 79 (5): 806.

[29]SHOW M D, HILL C M, ANWAY M D, WRIGHT W W, ZIRKIN B R. Phosphorylation of mitogen-activated protein kinase 8 (MAPK8) is associated with germ cell apoptosis and redistribution of the Bcl2-Modifying Factor (BMF). Journal of Andrology, 2008，29(3): 338-344.

[30]VALVE E, PENTTILA T L, PARANKO J, HA¨RKO¨NEN P. FGF-8 is expressed during specific phases of rodent oocyte and spermatogonium development. Biochemical and Biophysical Research Communications, 1997, 232(1) :173-177.

[31]DVORAK P, DVORAKOV D, HAMPL A. Fibroblast growth factor signaling in embryonic and cancer stem cells. FEBS Letters , 2006,580: 2869-2874.

[32]MIALON A, SANKINEN M, SO¨DERSTRO¨M H, JUNTTILA T T,HOLMSTRO¨M T, KOIVUSALO R, PAPAGEORGIOU A C,JOHNSON R S, HIETANEN S, ELENIUS K, WESTERMARCK J.DNA topoisomerase I is a cofactor for c-Jun in the regulation of epidermal growth factor receptor expression and cancer cell proliferation.Molecular and Cellular Biology, 2005, 12(25): 5040-5051.

[33]陳麗莉,肖亞梅,劉文彬,趙如榕,劉姣,劉筠,李萬(wàn)程.ERK/JNK 在黃鱔雌、雄發(fā)育階段生殖腺中的表達(dá)和定位.動(dòng)物學(xué)報(bào)，2007, 53(2):325-331.CHEN L L, XIAO Y M, LIU W B, ZHAO R R, LIU J, LI W C.Expression and location of ERK/JNK in ovary and spermary in rice field eels Monopterus albus. Journal of Animal Science, 2007,53(2):325-331. (in Chinese)

[34]SHEIKH A Q, TAGHIAN T, HEMINGWAY B, CHO H, KOGAN A B,NARMONEVA1 D A. Regulation of endothelial MAPK/ERK signalling and capillary morphogenesis by low-amplitude electric field.Journal of the Royal Society Interface, 2013, 10(78): 1-13.

[35]Vien H Y Lee, Leo T O Lee, Billy K C Chow. Gonadotropin- releasing hormone: regulation of the GnRH gene. FEBS Journal, 2008, 275:5458-5478.

[36]BELOV A A, MOHAMMADI M. Grb2, a double-edged sword of receptor tyrosine kinase signaling. Science Signaling, 2012, 5(249):49.

[37]HAINES E, MINOO P, FENG Z Q, RESALATPANAH N, NIE X M,CAMPIGLIO M, ALVAREZ L, COCOLAKIS E, RIDHA M,FULVIO M D, GOMEZ-CAMBRONERO J, LEBRUN J J, ALI S.Tyrosine phosphorylation of Grb2:Role in prolactin /epidermal growth factor cross talk in mammary epithelial cell growth and differentiation.Molecular and Cellular Biology, 2009, 29(10):2505-2520.

[38]YUAN Y Y, CHEN W Y, SHI Q X, MAO L Z,YU S Q, FANG X,ROLDAN E R S. Zona pellucida induces activation of phospholipase A2 during acrosomal exocytosis in guinea pig spermatozoa. Biology of Reproduction, 2003, 68(3): 904-913.

[39]SANTI C M, SANTOS T, HERNáNDEZ-CRUZ A, DARSZON A.Properties of a novel pH-dependent Ca2+permeation pathway present in male germ cells with possible roles in spermatogenesis and mature sperm function. The Journal of General Physiology, 1998,112: 33-53.

（責(zé)任編輯 林鑒非）

Transcriptome Analysis of Differential Gene Expression Associated with Testis Tissue in Mule Duck and Muscovy Duck

LI Li1, MIAO ZhongWei1, XIN QingWu1, ZHU ZhiMing1, ZHANG LinLi1,ZHUANG XiaoDong2, ZHENG NenZhu1,3
(1Institute of Animal Science and Veterinary Medicine, Fujian Academy of Agricultural Sciences, Fuzhou 350013;2Fujian Changlong Group, Zhangzhou 363000, Fujian;3Food College of Fujian Agriculture and Forestry University, Fuzhou 350002)

2017-02-04；接受日期：2017-06-13

福建省農(nóng)科院青年人才創(chuàng)新基金(YC2017-7)、福建省省屬公益類科研院所基本科研專項(xiàng)(2017R1023-5)、福建省農(nóng)科院所青年基金(MYQJ2015-5)

聯(lián)系方式：李麗，Tel：13960985616；E-mail：576801792@qq.com。通信作者鄭嫩珠，Tel：0591-83815170；E-mail：zhengnz@163.com