不同發(fā)育階段斯氏副柔線蟲比較轉(zhuǎn)錄組學(xué)分析

王文龍，馮陳晨，紅梅，岳建偉，呼和巴特爾，劉春霞

?

不同發(fā)育階段斯氏副柔線蟲比較轉(zhuǎn)錄組學(xué)分析

王文龍1，馮陳晨1，紅梅1，岳建偉1，呼和巴特爾1，劉春霞2

（1內(nèi)蒙古農(nóng)業(yè)大學(xué)獸醫(yī)學(xué)院/農(nóng)業(yè)部動物疾病臨床診療技術(shù)重點實驗室，呼和浩特 010018；2內(nèi)蒙古農(nóng)業(yè)大學(xué)生命科學(xué)學(xué)院，呼和浩特 010018）

探明不同發(fā)育階段駱駝斯氏副柔線蟲的轉(zhuǎn)錄組差異，了解不同發(fā)育階段蟲體在功能分類和代謝通路等方面的生物學(xué)特征，挖掘生長發(fā)育相關(guān)功能基因，豐富寄生性線蟲轉(zhuǎn)錄組學(xué)信息。采用Illumina HiSeq2000TM高通量測序技術(shù)對斯氏副柔線蟲蟲卵、第三期幼蟲和雌蟲進行轉(zhuǎn)錄組測序，構(gòu)建3個樣本的cDNA文庫，評估建庫質(zhì)量；利用Trinity軟件對所得序列進行De novo組裝及組裝效率評估；之后對獲得的有效序列進行功能注釋及相關(guān)生物信息學(xué)分析。測序及組裝后蟲卵、第三期幼蟲和雌蟲分別獲得47 717、76 342和54 624個Unigenes，在蟲卵和第三期幼蟲比對中共有33 579個差異表達基因，其中表達上調(diào)的基因有20 477個，表達下調(diào)的基因有13 102個；在雌蟲和第三期幼蟲比對中共有32 199個差異表達基因，其中表達上調(diào)的基因有9 293個，表達下調(diào)的基因有22 906個。將成對比較的差異表達基因分別進行GO功能分類，其中蟲卵和第三期幼蟲比對中有6 617、3 891和8 755個Unigenes，雌蟲和第三期幼蟲比對中有7 043、3 686和10 177個Unigenes分別注釋到生物過程、細胞組成和分子功能三大類中；通過KEGG pathway數(shù)據(jù)庫分析，蟲卵和第三期幼蟲比對中有6 521個差異表達基因參與到251條通路中，顯著富集MAPK信號通路、Wnt信號通路和氧化磷酸化等通路，雌蟲和第三期幼蟲比對中有6 528個參與到251條通路中，顯著富集新陳代謝通路、DNA復(fù)制和細胞周期等通路。差異表達基因功能聚類分析中，雌蟲和蟲卵在調(diào)控生長速率和生殖發(fā)育等方面高表達；第三期幼蟲在防御機制和糖類代謝等方面高表達；且3個發(fā)育階段均在胚胎發(fā)育，胚后發(fā)育中富集高表達，其中在胚胎發(fā)育和胚后發(fā)育中分別有242和202個相關(guān)功能基因在3個發(fā)育階段都表達，這些基因可能在胚胎（后）發(fā)育過程中起核心作用。此外，獲得了9種異時性基因（如：LIN-28、LIN-14和RHEB-1等）、48個核激素受體（NHRs），包括NHR-49，NHR-48，NHR-40，NHR-1等以及36個鋅金屬蛋白酶（NAS），包括NAS-36、NAS-33和NAS-14等，并對其在斯氏副柔線蟲蟲卵、第三期幼蟲和雌蟲中的富集程度進行分析，發(fā)現(xiàn)這些基因在維持線蟲正常生長發(fā)育過程中發(fā)揮關(guān)鍵作用，利用RNA-seq技術(shù)對3個發(fā)育階段的斯氏副柔線蟲進行測序和生物信息學(xué)分析，研究了不同發(fā)育階段蟲體的差異表達基因在GO功能分類、KEGG代謝通路和基因功能聚類等方面的生物學(xué)特性，篩選出多種異時性基因和發(fā)育相關(guān)重要基因，為后續(xù)開展斯氏副柔線蟲功能基因組學(xué)研究，蟲體與宿主互作、致病機制、免疫逃避等研究提供了理論依據(jù)。

斯氏副柔線蟲；轉(zhuǎn)錄組；差異表達基因；發(fā)育相關(guān)基因

0 引言

【研究意義】駱駝斯氏副柔線蟲（）屬于旋尾目、副柔屬[1-2]，是一種寄生于偶蹄反芻動物真胃的吸血性線蟲，駱駝是其最適宜終末宿主。大量感染斯氏副柔線蟲后，可引起駱駝腹瀉，貧血甚至死亡。據(jù)資料顯示，在中國內(nèi)蒙古巴彥淖爾市雙峰駝斯氏副柔線蟲的感染率高達91.7%，感染強度最高可達1 315條，嚴(yán)重威脅駱駝的健康[3]。2009年，趙治國等在吸血蠅體內(nèi)發(fā)現(xiàn)斯氏副柔線蟲的第三期幼蟲，并首次明確了截脈角蠅和西方角蠅是駱駝斯氏副柔線蟲病的傳播媒介[3]。但是有關(guān)駱駝斯氏副柔線蟲在傳播媒介與終末宿主體內(nèi)不同發(fā)育階段蟲體的代謝水平差異、發(fā)育相關(guān)重要基因表達及致病機理等方面的研究并未見報道，嚴(yán)重阻礙了駱駝斯氏副柔線蟲病的防控與治療。因此，開展駱駝斯氏副柔線蟲不同發(fā)育階段蟲體的差異表達基因研究對于從根本上解決斯氏副柔線蟲對駱駝的危害具有十分重要的意義。【前人研究進展】隨著高通量測序技術(shù)的快速發(fā)展，轉(zhuǎn)錄組測序已被廣泛應(yīng)用到不同發(fā)育階段生物個體的基因差異表達研究中。目前，寄生蟲轉(zhuǎn)錄組學(xué)研究也越來越受到重視。秀麗隱桿線蟲[4]（）基因組測序組織對不同發(fā)育階段的模式生物進行基因組學(xué)研究，構(gòu)建了基因組學(xué)圖譜并詳細注釋出發(fā)育相關(guān)的功能基因。Fu等[5]對犬惡絲蟲的轉(zhuǎn)錄組進行研究，組裝出20 810個轉(zhuǎn)錄本，并發(fā)現(xiàn)有1 101個是犬惡絲蟲特有的基因，為免疫抗原的發(fā)現(xiàn)提供幫助。Li等[6]對不同發(fā)育階段的馬來絲蟲進行轉(zhuǎn)錄組測序，詳細闡明了馬來絲蟲在不同發(fā)育階段的轉(zhuǎn)錄表達模式，及差異表達基因的功能。Laing等[7-8]繪制出捻轉(zhuǎn)血矛線蟲（）基因組及不同發(fā)育期轉(zhuǎn)錄本的草圖，發(fā)掘出重要的疫苗和藥物靶點。【本研究切入點】前人在多種寄生蟲發(fā)育轉(zhuǎn)錄組學(xué)研究中取得了階段性成果，但是對斯氏副柔線蟲的研究還很匱乏，尤其是對其不同發(fā)育階段轉(zhuǎn)錄組基因表達研究仍屬于未知。【擬解決的關(guān)鍵問題】本研究通過對不同發(fā)育階段斯氏副柔線蟲進行轉(zhuǎn)錄組測序，試圖揭示斯氏副柔線在傳播媒介與終末宿主體內(nèi)不同發(fā)育階段蟲體在代謝水平上的差異、參與的調(diào)控機制；發(fā)掘出吸血性線蟲特有功能基因及生長發(fā)育相關(guān)基因，為斯氏副柔線蟲病的相關(guān)基礎(chǔ)理論研究、診斷方法及防治研究奠定基礎(chǔ)。

1 材料與方法

試驗于2014—2016年在內(nèi)蒙古農(nóng)業(yè)大學(xué)獸醫(yī)學(xué)院完成。

1.1 3個不同發(fā)育階段斯氏副柔線蟲樣本采集

選取斯氏副柔線蟲蟲卵、第三期幼蟲（L3s）、雌蟲3個階段的蟲體進行轉(zhuǎn)錄組測序。雌蟲成蟲樣本為2014年和2015年11—12月間采集于內(nèi)蒙古巴盟烏拉特后旗的雙峰駝?wù)嫖钢?，在顯微鏡下鑒定出雌蟲，分裝標(biāo)記后在液氮中保存。第三期幼蟲為2014年和2015年的7—8月間采集于內(nèi)蒙古巴盟烏拉特后旗駝群環(huán)境中的吸血角蠅體內(nèi)，在顯微鏡下鑒定后保存于液氮中。蟲卵的收集通過將成年雌蟲置于37℃生理鹽水中過夜產(chǎn)卵，次日收集、鏡下鑒定和計數(shù)后液氮保存，用于RNA提取。

1.2 方法

1.2.1 RNA提取 參照Invitrogen 公司的Trizol Reagent說明書分別對不同發(fā)育階段斯氏副柔線蟲總RNA進行提取。使用RQ1酶消解RNA中的DNA，純化后將樣本稀釋進行UV檢測及1.5%普通瓊脂糖凝膠電泳質(zhì)檢合格后備用。

1.2.2 建立cDNA文庫及Illumina測序 利用oligodT-磁珠富集捕獲帶有polyA尾巴的mRNA。在高溫鹽離子作用下，mRNA被隨機打斷，修復(fù)并連接5′Adaptor，然后用帶有3’Adaptor和隨機六聚體的RT引物反轉(zhuǎn)錄合成cDNA。對反轉(zhuǎn)錄后的cDNA進行PCR擴增，擴增時引入barcode序列，最終選取片段大小為300—500 bp的PCR產(chǎn)物，利用Illumina HiSeq2000TM測序平臺進行測序。

1.2.3 測序數(shù)據(jù)分析 將測序所得的原始數(shù)據(jù)進行質(zhì)量評估和可信度分析，并去除測序過程中低質(zhì)量的序列和不確定的序列（Q＜20），將得到的Clean reads合并，利用Trinity軟件做轉(zhuǎn)錄組重頭組裝，對樣品組裝得到的Unigenes做進一步序列拼接、去冗余處理和同源聚類，最終得到轉(zhuǎn)錄本。將3個樣本Unigenes的表達量進行RPKM值歸一化處理。

在差異表達基因（differentially expressed genes，DEG）的篩選中，使用edgeR軟件進行兩兩樣本間的差異表達分析。檢測過程中，將差異倍數(shù)（fold change，F(xiàn)C）≥2 或≤0.5且值≤0.01作為篩選標(biāo)準(zhǔn)，并利用logCMP模型對兩個樣本進行標(biāo)準(zhǔn)化，這樣可以避免不明確的值和不明確的少數(shù)趨向于零的logFC，使兩樣本之間的比較更加詳細和全面。

1.2.4 差異表達基因注釋分析 斯氏副柔線蟲蟲卵和L3s，L3s和雌蟲相比較，將獲得的差異表達基因分別進行GO功能注釋，然后將其按細胞組分、分子功能和生物過程3個GO數(shù)據(jù)庫做功能分類及富集分析；同時將差異表達基因進行 KEGG 富集分析，把差異顯著的通路進行富集，找到不同發(fā)育階段蟲體內(nèi)顯著性差異

變化的生物學(xué)調(diào)控通路；將兩組比較的差異表達基因與基因編碼的蛋白分別進行比對，并注釋。然后用能夠被注釋的差異基因所對應(yīng)的蛋白GI號在DAVID平臺進行功能聚類分析，獲得可信度高的功能通路和更細致全面的功能聚類。

2 結(jié)果

2.1 斯氏副柔線蟲Unigenes聚類后數(shù)據(jù)可靠性分析

本研究對斯氏副柔線蟲蟲卵、第三期幼蟲（L3s）和成年雌蟲3個發(fā)育階段分別進行轉(zhuǎn)錄組測序，將得到的數(shù)據(jù)進行質(zhì)量控制、拼接組裝和聚類，總共獲得99 481個Cluster Unigenes用于差異基因分析。使用RPKM法計算各個發(fā)育階段Unigenes表達量。結(jié)果顯示，聚類后的Unigenes在47%—76%之間都有reads的分布，蟲卵、L3s、雌蟲中RPKM值大于10的Unigenes占所有表達Unigenes的9%—22%左右,說明聚類后數(shù)據(jù)可靠性較高（表1）。

2.2 不同發(fā)育階段斯氏副柔線蟲的差異表達基因相關(guān)分析結(jié)果

2.2.1 不同發(fā)育階段斯氏副柔線蟲的差異表達基因篩選結(jié)果 選用edgeR軟件分別對斯氏副柔線蟲L3s和蟲卵、雌蟲和L3s比較的基因表達量做差異表達分析，以-value≤0.01及Fold Change≥2或≤0.5為標(biāo)準(zhǔn)篩選差異表達基因（圖1）。圖1中紅色的點表示顯著差異表達基因（DEGs），縱坐標(biāo)為logFC表示某一個基因在兩樣本中表達量差異倍數(shù)的對數(shù)值，且絕對值越大，表明基因表達量變化的倍數(shù)越大；橫坐標(biāo)為logCPM表示兩樣本之間比較同一個基因時的總聚集點，且坐標(biāo)值越大，表明篩選的差異表達基因越可靠。

表1 每個樣本表達的Unigene聚類

a表達的Unigenes數(shù)占總參考基因組Unigenes總數(shù)的比例；bRPKM值大于等于10的Unigenes數(shù)占RPKM值大于0的Unigenes的比例

圖1 edgeR方法鑒定蟲卵和L3s、雌蟲和L3s的差異表達基因

2.2.2 不同發(fā)育階段斯氏副柔線蟲差異基因表達結(jié)果 不同發(fā)育階段斯氏副柔線蟲成對比較后，將得到的差異表達基因進行統(tǒng)計，結(jié)果表明，L3s和蟲卵相比，差異表達基因共有33 579個，表達上調(diào)的基因有20 477個，其中有7 561個上調(diào)基因注釋出同源蛋白；表達下調(diào)的基因有13 102個，其中有2 645個下調(diào)基因注釋出同源蛋白。雌蟲和L3s比對時，差異表達基因共有32 199個，表達上調(diào)的基因有9 293個，其中3 874個上調(diào)基因注釋出同源蛋白，表達下調(diào)的基因有22 906個，其中注釋出6 384個下調(diào)基因注釋出同源蛋白（圖2）。

2.3 不同發(fā)育階段斯氏副柔線蟲差異表達基因Gene Ontology分析

2.3.1 不同發(fā)育階段斯氏副柔線蟲差異表達基因GO-生物過程富集分析 在L3s和蟲卵的差異表達基因GO-生物過程富集分析中，有6 617個差異表達基因比對到BP數(shù)據(jù)庫中，顯著富集在蛋白質(zhì)氨基酸磷酸化過程中的差異表達基因（DEG）有119個，占該生物過程基因數(shù)（EG）的51.74%，占比對到BP數(shù)據(jù)庫所有差異表達基因數(shù)的1.80%；顯著富集在翻譯過程中的差異表達基因有130個，占該生物過程基因數(shù)的46.76%，占比對到BP數(shù)據(jù)庫所有差異表達基因數(shù)的1.97%。雌蟲和L3s比對中，有7 043個差異表達基因比對到BP數(shù)據(jù)庫中，除了在翻譯、肌肉組織發(fā)育等過程中顯著富集，在染色體結(jié)構(gòu)中有9個差異表達基因，占該生物學(xué)過程基因數(shù)的81.82%，占比對到BP數(shù)據(jù)庫所有差異表達基因數(shù)的0.13%；有8個差異基因顯著富集在性腺發(fā)育過程中，占該生物學(xué)過程基因數(shù)的80.00%，占比對到BP數(shù)據(jù)庫所有差異表達基因數(shù)的0.11%（表2）。

1：下調(diào)差異表達基因 Up-regulated DEGs；2：上調(diào)差異表達基因 Down-regulated DEGs

2.3.2 不同發(fā)育階段差異表達基因GO-細胞組成富集分析 L3s和蟲卵GO-細胞組分富集分析中，有3 891個差異表達基因比對到CC數(shù)據(jù)庫，其中，有90個顯著富集在細胞溶質(zhì)中，占該細胞組分基因數(shù)的61.64%，占比對到CC數(shù)據(jù)庫所有差異表達基因數(shù)的2.31%；顯著富集在細胞核中的差異表達基因有215個，占該細胞組分基因數(shù)的47.46%，占比對到CC數(shù)據(jù)庫所有差異表達基因數(shù)的5.53%；顯著富集在線粒體中的差異表達基因有58個，占該細胞組分基因數(shù)的55.77%，占比對到 CC數(shù)據(jù)庫所有差異表達基因數(shù)的1.49%。雌蟲和L3s比對中，有3 686個差異表達基因比對到CC數(shù)據(jù)庫。除了基本細胞組分的富集之外，有6個差異表達基因在細胞內(nèi)無膜細胞器中顯著富集，占該細胞組分基因數(shù)的100%，占比對到CC數(shù)據(jù)庫所有差異表達基因數(shù)的0.16%。有20個差異表達基因顯著富集在膠原蛋白中，占該細胞組分基因數(shù)的58.82%，占比對到CC數(shù)據(jù)庫所有差異表達基因數(shù)的0.54%。

2.3.3 不同發(fā)育階段差異表達基因GO-分子功能富集分析 L3s和蟲卵GO-分子功能富集分析中，有8 755個差異表達基因比對到MF數(shù)據(jù)庫。其中，有172個差異基因顯著富集在蛋白結(jié)合功能中，占該分子功能基因數(shù)的47.78%，占比對到MF數(shù)據(jù)庫所有差異表達基因數(shù)的1.97%；有77個顯著富集在蛋白絲氨酸/蘇氨酸激酶活性中，占該細胞組分基因數(shù)的52.38%，占比對到MF數(shù)據(jù)庫所有差異表達基因數(shù)的0.88%；顯著富集在亞鐵血紅素結(jié)合中的差異表達基因有38個，占該細胞組分基因數(shù)的48.10%，占比對到MF數(shù)據(jù)庫所有差異表達基因數(shù)的0.43%。雌蟲和L3s比對中，有10 177個差異表達基因比對到MF數(shù)據(jù)庫，在蛋白結(jié)合、鋅離子結(jié)合和亞鐵血紅素結(jié)合等分子功能中顯著富集的同時，有18個差異表達基因在氫離子跨膜轉(zhuǎn)運活性中顯著富集，占該分子功能基因數(shù)的69.23%，占比對到MF數(shù)據(jù)庫所有差異表達基因數(shù)的0.18%。有18個顯著富集在表皮結(jié)構(gòu)組成中，占該細胞組分基因數(shù)的69.23%，占比對到MF數(shù)據(jù)庫所有差異表達基因數(shù)的0.18%（表2）。

2.4 不同發(fā)育階段斯氏副柔線蟲差異基因KEGG Pathway分析

將不同發(fā)育階段斯氏副柔線蟲的差異基因注釋到KEGG數(shù)據(jù)庫，結(jié)果顯示，L3s和蟲卵中比對到差異表達基因12 803個，其中6 521個差異表達基因有具體的定義，并顯著富集到255條通路。其中氧化磷酸化通路顯著富集221個差異表達基因，占該通路基因數(shù)的64.62%；MAPK信號通路顯著富集100個差異表達基因，占該通路基因數(shù)的63.69%；Wnt信號通路顯著富集90個差異表達基因，占該通路基因數(shù)的58.82%。雌蟲和L3s中比對到13 153個差異表達基因，其中有6 528個差異表達基因有具體的定義，同時涉及251條通路。其中代謝通路顯著富集有1 679個差異表達基因，占該通路基因數(shù)的51.61%，細胞循環(huán)通路顯著基因有60個，占該通路的63.16%（表3）。

2.5 差異表達基因功能聚類

將斯氏副柔線蟲L3s和蟲卵、雌蟲和L3s相比的差異表達基因利用DAVID平臺進行功能聚類分析。結(jié)果顯示，有關(guān)防御機制；己糖代謝、葡萄糖代謝、糖酵解/糖異生等功能聚類在L3s中富集程度明顯高于蟲卵和雌蟲階段，這可能與L3s期幼蟲在感染哺乳動物宿主時，采取自身免疫保護和免疫逃避有關(guān)。雌蟲和蟲卵階段在生長調(diào)控速率和生殖發(fā)育等相關(guān)功能聚類中富集性高表達，而L3s中富集不明顯。雌蟲在性別分化、生殖器發(fā)育等功能聚類中顯著性高表達。蟲卵、三期幼蟲和雌蟲在胚胎發(fā)育、胚后發(fā)育以及幼蟲發(fā)育等功能聚類中都富集性高表達。在胚胎發(fā)育功能聚類中有242個相關(guān)功能基因在3個發(fā)育階段都表達，而在蟲卵、L3s和雌蟲中分別有63個、571個和248個特異性表達的功能基因。在胚后發(fā)育功能聚類中有202個相關(guān)功能基因富集在3個發(fā)育階段，詳情見圖3。

表2 L3s和Egg，雌蟲和L3s差異表達基因GO-ontology富集分析

1)CF=DEG/TDEG×100%,CF：在該ontology中差異表達基因簇的頻率；DEG：在該ontology中差異表達基因的數(shù)量；TDEG:比對到ontology中所有差異表達基因數(shù)。2)GF=EG/TEG×100%,GF:比對到該ontology中基因的頻率；EG：比對到該ontology中所有基因數(shù)量；TEG：比對到ontology中所有基因數(shù)

1)CF=DEG/TDEG×100%, CF is cluster frequency of differentially expressed gene annotated to each ontology; DEG is the numbers of differentially expressed gene annotated to each ontology; TDEG is numbers of all differentially expressed genes annotated to GO ontology.2)GF=EG/TEG×100%,GF is the genome frequency of all genes annotated to the ontology; EG is the numbers of genes annotated to each ontology; TEG is the numbers of all genes annotated to GO ontology

表3 L3s和蟲卵，雌蟲和L3s中差異表達基因顯著性富集的通路

DEGs：在每個KEGG代謝通路中差異表達基因數(shù)量；KEGG：在KEGG數(shù)據(jù)庫中涉及到該通路的基因數(shù)

DEGs: The numbers of Differentially expressed genes in each KEGG pathway; KEGG: Gene numbers in this pathway in KEGG database

圖3 蟲卵、幼蟲和雌蟲在胚胎發(fā)育和胚后發(fā)育過程中相關(guān)功能基因的韋恩圖分布

2.6 不同發(fā)育階段斯氏副柔線蟲生長發(fā)育相關(guān)基因分析

本文以影響生長發(fā)育相關(guān)的基因為參考，分析出不同發(fā)育階段斯氏副柔線蟲中生長發(fā)育相關(guān)的重要基因并對其進行富集分析。結(jié)果顯示，在蟲卵、L3s和雌蟲3個發(fā)育階段中都出現(xiàn)了不同表達程度的異時性相關(guān)基因: 核受體DAF-12、 LIN-12和BLMP-1基因在L3s中高表達；DRE-1/FBXO11和LIN-42基因在雌蟲中高表達；LIN-29基因在蟲卵中較高表達，在L3s和雌蟲中幾乎不表達；RHEB-1和LIN-28在蟲卵和雌蟲中高表達，在L3中表達量較低；LIN-14在3個發(fā)育階段中都有表達，其中在L3s中表達量相對較高。

本研究分析了核激素受體（nuclear hormone receptor，NHRs）和鋅金屬蛋白酶（NAS）兩類在線蟲生長發(fā)育過程中重要的發(fā)育基因。NHRs是生物體內(nèi)發(fā)育和代謝過程中重要的調(diào)節(jié)者，在斯氏副柔線蟲3個發(fā)育階段中共注釋出48個核激素受體同源蛋白，其中有26個具有詳細的功能分類，主要包括：NHR-49、NHR-48、NHR-40、NHR-1等，其中NHR-40和NHR-49在L3s期中較高表達，NHR-48在雌蟲和L3s中相對高表達，NHR-1在蟲卵和雌蟲期中較高表達。鋅金屬蛋白酶，又可稱作線蟲蝦紅素（nematode astacin，NAS）對線蟲表皮合成和表皮膠原蛋白酶裂解有重要作用。在斯氏副柔線蟲3個發(fā)育階段中注釋出36個NAS，其中有26個有詳細的功能分類，主要包括NAS-36、NAS-33和NAS-14等。其中NAS- 36在L3s期中相對高表達，在蟲卵期和雌蟲期未表達；NAS-15在蟲卵和雌蟲期的表達量高于L3s期，NAS-14在L3s期的表達量高于蟲卵和雌蟲期，具體詳情見表4。

表4 蟲卵期、幼蟲期和雌蟲期生長發(fā)育相關(guān)的重要基因

3 討論

新一代高通量測序技術(shù)的不斷發(fā)展，已徹底改變了轉(zhuǎn)錄組學(xué)的研究，使RNA-Seq無需預(yù)先設(shè)計探針即可對特定條件下任意生物生長發(fā)育階段整體轉(zhuǎn)錄活動進行測序，并且探測各種條件下的基因表達情況，發(fā)現(xiàn)了許多未知的研究領(lǐng)域[9-10]。本研究針對蟲卵、第三期幼蟲和雌蟲3個發(fā)育階段的斯氏副柔線蟲進行轉(zhuǎn)錄組測序分析，探索不同發(fā)育階段中基因表達的差異、代謝通路的差異及生長發(fā)育過程中功能聚類的顯著差異。樣本采集的代表性對轉(zhuǎn)錄組數(shù)據(jù)的準(zhǔn)確性和代表性起至關(guān)重要作用，所以該研究中樣本采集是重要的質(zhì)控過程。由于不同采集時間和樣本自身特點等因素，導(dǎo)致不同時間采集的相同發(fā)育階段蟲體基因表達量可能存在一定的差異。因此，本研究為了使獲取的轉(zhuǎn)錄本數(shù)據(jù)更加全面，將三個樣本的采集次數(shù)均增加至兩年中的9—10次；同時采集樣本數(shù)量也增加到蟲卵80萬只左右，第三期幼蟲5 000只左右，雌蟲200只左右。通過將多次采集時間的大量樣本進行混合上機測序，以保證轉(zhuǎn)錄本測序數(shù)據(jù)的全面性和代表性，提高數(shù)據(jù)的重復(fù)性和多樣性。

3.1 差異表達基因的GO功能富集分析

不同發(fā)育階段斯氏副柔線蟲的差異表達基因在GO功能富集分析中顯示，L3s和蟲卵相比，顯著富集到負調(diào)控細胞程序性死亡通路；而幼蟲到雌蟲富集明顯減少，可能與其為了生長發(fā)育，加快代謝通路有關(guān)。雌蟲和L3s相比，富集到肌肉器官發(fā)育、生殖腺發(fā)育、膠原蛋白和基于角質(zhì)的表皮發(fā)育信號通路；而L3s和蟲卵中富集明顯減少，說明在幼蟲到雌蟲開始攝血與膠原蛋白、表皮發(fā)育、生殖腺發(fā)育相關(guān)的基因表達顯著增加，這與SCHWARZ 等[8]對捻轉(zhuǎn)血矛線蟲的研究一致。在GO-細胞組分?jǐn)?shù)據(jù)庫中，L3s和蟲卵相比，顯著富集到細胞溶質(zhì)、細胞核和膠原蛋白中；雌蟲和L3s相比，顯著富集到細胞核、核糖體及細胞內(nèi)有膜細胞器中，而細胞核、核糖體與蟲體生長發(fā)育有關(guān)。在GO-分子功能數(shù)據(jù)庫中，L3s和蟲卵相比，顯著富集到蛋白質(zhì)結(jié)合功能、鈣離子結(jié)合功能和蛋白激酶活性中；雌蟲和幼蟲相比，顯著富集到鋅離子結(jié)合功能、氫離子跨膜轉(zhuǎn)運蛋白活性及ATP結(jié)合功能，而蛋白激酶活性和氫離子跨膜轉(zhuǎn)運蛋白活性與生長發(fā)育過程中能量消耗有關(guān)，表明在不同發(fā)育階段蟲體代謝耗能的主要方式不同，而ATP結(jié)合功能在幼蟲到雌蟲中顯著富集，在蟲卵到L3s中富集明顯減少，可能與其發(fā)育階段不同所需能量不同有關(guān)。從蟲卵到L3s、L3s到雌蟲的發(fā)育過程中亞鐵血紅素結(jié)合功能都顯著富集，根據(jù)報道秀麗隱桿線蟲中存在與脊椎動物（如鴿子、豬等）SCS-β亞基（丁二酰輔酶A連接酶）同源性較高的基因，GTP依賴型SCS在脊椎動物中參與三羧酸循環(huán)中的可逆反應(yīng)，反向反應(yīng)激活酮體亞鐵血紅素的合成[11-12]，這與GO富集在亞鐵血紅素結(jié)合功能的結(jié)果一致。

3.2 差異表達基因的KEGG pathway富集分析

通過KEGG pathway對差異表達基因進行顯著性富集分析，將差異顯著的 pathway 進行富集，有助于找到不同發(fā)育階段蟲體中顯著性差異變化的生物學(xué)調(diào)控通路。KEGG注釋和聚類分析結(jié)果顯示，L3s和蟲卵相比富集到Wnt/MAPK信號通路，其通路中起重要作用的lit-1基因上調(diào)表達；根據(jù)報道Wnt/MAPK信號通路參與秀麗線蟲側(cè)線細胞的時序分化調(diào)控，如細胞命運特化和對稱/不對稱分裂[13-14]，線蟲中l(wèi)it-1基因的缺失和增強會引起發(fā)育遲緩和發(fā)育過早的異時性缺陷[15-16]，表明Wnt/MAPK信號通路富集在L3s發(fā)育階段維持正常發(fā)育速度。雌蟲和L3s相比主要富集在與生化代謝有關(guān)的代謝通路中，如嘌呤代謝、嘧啶代謝；與遺傳信息有關(guān)的DNA復(fù)制、RNA聚合酶和細胞周期通路中；與糖異生信號有關(guān)通路中。已報道捻轉(zhuǎn)血矛線蟲L1到L3時期以及秀麗隱桿線蟲Dauer時期糖異生作用顯著增強[17]，而斯氏副柔線蟲在雌蟲和幼蟲相比中有顯著富集，可能由于物種差異導(dǎo)致代謝通路在時空性上存在差異。

3.3 差異表達基因功能聚類分析

本研究利用DAVID平臺進行功能聚類分析，DAVID的功能聚類數(shù)據(jù)庫整合了Gene Ontology、Interpro、KEGG等基因功能數(shù)據(jù)庫，獲得的功能通路可信度高并且全面細致。同時將繁多的聚類獲得的Clusters 進行翻譯和歸類，以便進行不同樣品之間的比較。功能聚類結(jié)果分析顯示，雌蟲在生殖發(fā)育、性別分化、生殖器發(fā)育、雌雄同體的生殖器發(fā)育中相關(guān)基因高表達，而在L3s中低表達，其中雌雄同體生殖器發(fā)育相關(guān)基因?qū)τ诰€蟲的性別分化發(fā)揮重要的功能，與性別有關(guān)的基因小窩蛋白（CAV-1）在雌蟲中上調(diào)表達，這與秀麗隱桿線蟲發(fā)現(xiàn)CAV-1在胚胎和生殖細胞中高表達[18]以及CAV-1在旋毛蟲中高表達的研究一致[19]。在L3s和蟲卵中，發(fā)現(xiàn)ACT-4基因在L3s和卵中富集性高表達，而在雌蟲和L3s中富集性低表達，根據(jù)報道在秀麗隱桿線蟲中，ACT-4基因在低氧脅迫條件下表達量上升，與LEV-11、MLC-1一起發(fā)揮細胞骨架的功能[20]，所以在卵和L3s中可能由于低氧環(huán)境ACT-4基因表達量上升。

3.4 不同發(fā)育階段斯氏副柔線蟲生長發(fā)育相關(guān)基因分析

以模式生物秀麗隱桿線蟲為參考，發(fā)現(xiàn)斯氏副柔線三個發(fā)育階段中表達出不同程度的異時性相關(guān)基因。經(jīng)分析發(fā)現(xiàn)，LIN-14基因存在于斯氏副柔線蟲3個發(fā)育階段中，這可能與LIN-14在胚胎發(fā)育后期分開調(diào)控發(fā)育時期特異性活動有關(guān)，且LIN-14過早或遲緩發(fā)育的突變可使雌蟲期外陰發(fā)育遲緩，影響產(chǎn)卵系統(tǒng)[21]。LIN-14抑制LIN-12活性，LIN-14通過LIN-12對外陰前體細胞（VPCs）的時間控制來作用于VPCs的空間發(fā)育模式，兩者共同調(diào)控外陰前體細胞的命運[22]。RHEB-1基因在斯氏副柔線蟲蟲卵和雌蟲期表達量較高，這可能與其調(diào)控線蟲壽命長度的功能有關(guān)[23]；LIN-42基因在斯氏副柔線蟲雌蟲期相對高表達，這可能與LIN-42調(diào)控線蟲性腺發(fā)育，外陰及性肌母細胞的發(fā)育時序功能有關(guān)[24]；BLMP-1和 DRE-1/FBOX11在斯氏副柔線蟲3個發(fā)育階段的表達量呈負相關(guān)，這與2014年Horn等發(fā)現(xiàn)通過鋅指蛋白BLMP-1 與DRE-1/FBOX11結(jié)合而被降解來控制秀麗隱桿線蟲的時空發(fā)育、dauer形成、蛻皮、性腺的成熟和壽命的研究一致[25]，說明兩個蛋白之間呈負相關(guān)相互調(diào)控，在后生生物成熟過程中發(fā)揮作用。2015年WANG等[26]對線蟲類固醇激素受體DAF-12基因研究發(fā)現(xiàn)，在適宜的環(huán)境條件下，DAF-12被DA（dafachronic acids）激活，DAF-12的活化作用能夠協(xié)調(diào)動能的儲存，從而使幼蟲生長繁殖；在斯氏副柔線蟲L3s期中DAF-12基因高表達，由于L3s在中間宿主體內(nèi)適應(yīng)寄生環(huán)境后，L3s中DAF-12被活化，存儲動能使幼蟲生長發(fā)育。HUANG等[27]研究發(fā)現(xiàn)頂端細胞（DTCs）對性腺的形成有重要的作用，BLMP-1負調(diào)控DTCs背向彎曲，而DAF-12，DRE-1，LIN-29功能冗余，正調(diào)控背向彎曲。DAF-12和LIN-29抑制BLMP-1轉(zhuǎn)錄，而DRE-1結(jié)合到BLMP-1促進其降解，使BLMP-1負調(diào)控獲得正確調(diào)控時間。在器官形成期間，不同的異時性基因?qū)r間的和空間的信號整合成一個回路共同協(xié)調(diào)生物體整體發(fā)育。這些異時性相關(guān)基因形成了復(fù)雜的調(diào)控網(wǎng)絡(luò)，調(diào)控線蟲整個發(fā)育過程的時間性和空間性，保證各種器官的準(zhǔn)確定位和適時發(fā)育。核激素受體（NHRs）是調(diào)控應(yīng)對發(fā)育、環(huán)境及營養(yǎng)狀況信號的基因表達蛋白，NHR對很多發(fā)育過程有重要作用，其中NHR-49 具有調(diào)節(jié)和控制脂肪代謝、維持脂肪酸飽和的正常平衡功能[28]；在腸道中抑制溶酶體脂質(zhì)的累積，與PKG信號通路共同調(diào)控短期禁食誘導(dǎo)溶酶體脂質(zhì)的積累[29]；并且NHR-49調(diào)控胚胎發(fā)育、幼蟲發(fā)育及成蟲壽命。NHR-49在斯氏副柔線蟲L3s期中較高表達，可能與其調(diào)控幼蟲發(fā)育有關(guān)。NHR-48參與轉(zhuǎn)錄本調(diào)控并具有序列特異的DNA結(jié)合活性和轉(zhuǎn)錄因子活性等功能。根據(jù)BROZOVA等[30]研究發(fā)現(xiàn)NHR-40主要調(diào)控晚期胚胎和早期幼蟲的延長和形態(tài)發(fā)生，且NHR-40功能表型缺失會導(dǎo)致不規(guī)律的體壁肌肉細胞發(fā)育，損害運動和神經(jīng)肌肉協(xié)調(diào)。在線蟲的不同幼蟲期發(fā)育中鋅金屬蛋白酶對矯正發(fā)育、表皮適當(dāng)蛻皮起著至關(guān)重要的作用[31]。NAS-36是控制蛻皮的關(guān)鍵基因，其突變后可擾亂蛻皮過程，使線蟲無法進入正常的發(fā)育階段，并且NAS-36還可以通過加工特定的細胞外基質(zhì)蛋白調(diào)控新的表皮；在寄生性線蟲期間NAS-36也是關(guān)鍵的靶基因，可用來控制殺蟲效應(yīng)。在斯氏副柔線蟲中NAS-36主要在L3s期高表達，根據(jù)Stepek等[32-33]研究發(fā)現(xiàn)，NAS-36基因無論是在自生性線蟲還是在寄生性線蟲中基因功能都較保守，所以NAS-36在L3s期中也是控制蛻皮的關(guān)鍵基因。

4 結(jié)論

本研究繪制出斯氏副柔線蟲在蟲卵、第三期幼蟲和雌蟲3個發(fā)育階段的差異基因表達譜；分析了斯氏副柔線蟲在傳播媒介和終末宿主體內(nèi)不同發(fā)育階段代謝水平上的差異；探索了不同發(fā)育階段蟲體之間重要的功能聚類，并注釋出生長發(fā)育相關(guān)的功能基因，為斯氏副柔線蟲相關(guān)理論研究、藥物靶點發(fā)掘、免疫學(xué)診斷和防治研究奠定了基礎(chǔ)。

[1] 宋銘忻, 張龍現(xiàn). 獸醫(yī)寄生蟲學(xué). 北京: 科學(xué)出版社, 2009: 183-184.

SONG M X, ZHANG L X.Beijing: Science Press, 2009: 183-184. (in Chinese)

[2] 黃兵, 沈杰. 中國畜禽寄生蟲形態(tài)分類圖譜. 北京: 中國農(nóng)業(yè)科學(xué)技術(shù)出版社, 2006: 462-463.

HUANG B, SHEN J.Beijing: China Agriculture Press, 2006: 462-463. (in Chinese)

[3] 趙治國. 我國駱駝斯氏副柔線蟲病傳播媒介的研究[D]. 內(nèi)蒙古農(nóng)業(yè)大學(xué), 2010.

Zhao Z G. Study on the vector of camel parabronemosis in China[D]. Inner Mongolia: Inner Mongolia Agricultural University, 2010. (in Chinese)

[4] elegans Sequencing ConsortiumC. Genome sequence of the nematode: a platform for investigating biology., 1998, 282: 2012–2018.

[5] Fu Y, Lan J, Zhang Z, Hou R, Wu X, Yang D, Zhang R, Zheng W, Nie H, Xie Y, Yan N, Yang Z, Wang C, Luo L, Liu L, Gu X, Wang S, Peng X, Yang G. Novel insights into the transcriptome of, 2012, 7(7): e41639.

[6] Li B W, Wang Z Y, Rush C A, Mitreva M, Weil J G. Transcription profiling reveals stage-and function-dependent expression patterns in the filarial nematode, 2012, 13: 184.

[7] Laing R, Kikuchi T, Martinelli A, Tsai IJ, Beech R N, Redman E, Holroyd N, Bartley D J, Beasley H, Britton C, Curran D, Devaney E, Gilabert A, Hunt M, Jackson F, Johnston SL, Kryukov I, Li K, Morrison A A, Reid A J, Sargison N, Saunders G I, Wasmuth J D, Wolstenholme A, Berriman M, Gilleard J S, Cotton J A. The genome and transcriptome of, a key model parasite for drug and vaccine discovery., 2013, 14: R88.

[8] Schwarz E M, Korhonen P K, Campbell B E, Young N D, Jex A R, Jabbar A, Hall R S, Mondal A, Howe A C, Pell J, Hofmann A, Boag P R, Zhu X Q, Gregory T, Loukas A, Williams B A, Antoshechkin I, Brown C, Sternberg P W, Gasser R B. The genome and developmental transcriptome of the strongylid nematode., 2013, 14(8): R89.

[9] 劉嬌, 張建珍, 李大琪, 張婷婷, 馬恩波, 張建琴. 中華稻蝗羧酸酯酶家族基因生物信息學(xué)及組織表達特異性分析. 中國農(nóng)業(yè)科學(xué), 2015, 48(21): 4272-4284.

LIU J, ZHANG J Z, LI D Q, ZHANG T T, MA E B, ZHANG J Q. Bioinformatics and tissue-specific expression analysis of carboxylesterase genes from., 2015, 48(21): 4272-4284. (in Chinese)

[10] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 張曌楠, 黃枳腱, 張璐, 王鴻權(quán), 解彥玲, 童新宇. 中華蜜蜂幼蟲腸道響應(yīng)球囊菌早期脅迫的轉(zhuǎn)錄組學(xué). 中國農(nóng)業(yè)科學(xué), 2017, 50(13): 2614-2623.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, ZHANG Z N, HUANG Z J, ZHANG L, WANG H Q, XIE Y L, TONG X Y. Transcriptome oflarval gut under the stress of., 2017, 50(13): 2614-2623. (in Chinese)

[11] Johnson J D, Mehus J G, Tews K, Milavetz B I, Lambeth D O. Genetic evidence for the expression of ATP- and GTP-specific succinyl-CoA synthetases in multicellular eucaryotes.1998, 273(42): 27580-27586.

[12] Przybyla Z B, Dennis R A, Zakharkin S O, McCammon M T. Genes of succinyl-CoA ligase from., 1998, 258(2): 736-743.

[13] Delawary M, Nakazawa T, Tezuka T, Sawa M, Iino Y, Takenawa T, Yamamoto T. Molecular characterization of a novel RhoGAP, RRC-1 of the nematode., 2007, 357(2): 377-382.

[14] Yang X D, Karhadkar T R, Medina J, Robertson S M, Lin R. β-Catenin-related protein WRM-1 is a multifunctional regulatory subunit of the LIT-1 MAPK complex., 2015, 112(2): E137-146.

[15] Chuang M H, Chiou S H, Huang C H, Yang W B, Wong C H. The lifespan-promoting effect of acetic acid and Reishi polysaccharide., 2009, 17(22): 7831-7840.

[16] Jackson B M, Abete L P, Krause M W, Eisenmann D M. Use of an activated beta-catenin to identify Wnt pathway target genes in, including a subset of collagen genes expressed in late larval development., 2014, 4(4): 733-747.

[17] John W, Stuart K K. Global analysis of dauer gene expression in., 2003, 130: 1621-1634.

[18] Parker S, Baylis H A. Overexpression of caveolins ininduces changes in egg-laying and fecundity., 2009, 2(5): 382-384.

[19] Romel H B, Ricardo R N, Saé M H, Karen N C, Lenin P, Ana G S, Jorge M M. Sex steroids effects on the molting process of the helminth human parasite

[20] Li H, Ren C, Shi J, Hang X, Zhang F, Gao Y, Wu Y, Xu L, Chen C, Zhang C. A proteomic view ofcaused by short-term hypoxic stress., 2010, 8: 49.

[21] Johnson R W, Liu L Y, Hanna R W, Chamberlin H M. Theheterochronic gene lin-14 coordinates temporal progression and maturation in the egg-laying system., 2009, 238(2): 394-404.

[22] Li J, Greenwald I. LIN-14 inhibition of LIN-12 contributes to precision and timing ofvulval fate patterning., 2010, 20(20): 1875-1879.

[23] Honjoh S, Yamamoto T, Uno M, Nishida E. Signalling through RHEB-1 mediates intermittent fasting-induced longevity in., 2009, 457(7230): 726-730.

[24] Tennessen J M, Gardner H F, Volk M L, Rouqvie A E. Novel heterochronic functions of theperiod- related protein LIN-42., 2006, 289(1): 30-43.

[25] Horn M, Geisen C, Cermak L, Becker B, Nakamura S, Klein C, Pagano M, Antebi A. DRE-1/FBXO11-dependent degradation of BLMP-1/BLIMP-1 Governsdevelopmental timing and maturation., 2014, 28(6): 697-710.

[26] Wang Z, Jonathan S, You Y J, Ranjit N, Tang H, Xie Y, Lok J B, Mangelsdorf D J, Kliewer S A. The nuclear receptor DAF-12 regulates nutrient metabolism and reproductive growth in nematodes., 2015, 11(3): e1005027.

[27] Huang T F, Cho C Y, Cheng Y T, Huang J W, Wu Y Z, Yeh A Y, Nishiwaki K, Chang S C, Wu Y C. BLMP-1/Blimp-1 Regulates the spatiotemporal cell migration pattern in, 2014, 10(6): e1004428.

[28] Van G M R, Hadjivassiliou H, Jolly A, Yamamoto K R. Nuclear hormone receptor NHR-49 controls fat consumption and fatty acid composition in.2005, 3(2): e53.

[29] Huang W M, Li Z Y, Xu Y J, Wang W, Zhou M G, Zhang P, Liu P S, Xu T, Wu Z X. PKG and NHR-49 signalling co-ordinately regulate short-term fasting-induced lysosomal lipid accumulation in., 2014, 461(3): 509-520.

[30] Brozova Z, Simeckova K, Kostrouch Z, Rall J E, Kostrouchova M. NHR-40, asupplementary nuclear receptor, regulates embryonic and early larval development., 2006, 123(9): 689-701.

[31] Suzuki M, Sagoh N, Iwasaki H, Inoue H, Takahashi K. Metalloproteases with EGF, CUB, and thrombospondin-1 domains function in molting of., 2004, 385(6): 565-568.

[32] Stepek G, McCormack G, Birnie A J, Page A P. The astacin metalloprotease moulting enzyme NAS-36 is required for normal cuticle ecdysis in free-living and parasitic nematode., 2011, 138(2): 237-248.

[33] Sharma O P, Agrawal S, Kumar M S. Physicochemical properties of the modeled structure of astacin metalloprotease moulting enzyme NAS-36 and mapping the druggable allosteric space of,andvia molecular dynamics simulation., 2013, 5(4): 312-323.

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

The Comparative Transcriptome Analysis ofat Different Developmental Stages

WANG WenLong1, FENG ChenChen1, HONG Mei1, YUE JianWei1, Huhebateer1, LIU ChunXia2

(1College of Veterinary Medicine, Inner Mongolia Agricultural University/Key Laboratory of Clinical Diagnosis and Treatment Technology in Animal Disease, Ministry of Agriculture, Hohhot 010018 ;2College of Life Sciences, Inner Mongolia Agricultural University, Hohhot 010018)

The objective of this study was to identify the differentially expressed genes(DEG) and describe biological characteristics involved in functional classifications and metabolic pathways at different developmental stages ofinfecting camel, which is necessary to better understand functional genes involved in growth and development and enrich the transcriptome data of parasitical nematodes.Eggs, the third-stage larvae(L3s) and females ofwere sequenced by Illumina HiSeq2000TMsequencing platform and constructed their cDNA libraries after quality filtering. De novo assembling and assembly efficiency assessment were carried out using Trinity, a short-read assembly program. Then all the effective sequential data obtained wereassigned to the relevant databases to perform functional annotation and bioinformatic analysis.The results showed that 47 717, 76 342 and 54 624 unigenes were obtained respectively in eggs，the third-stage larvae and female stages. 33 579 differentially expressed genes(DEGs) were identified by comparing the unigenes obtained from eggs and L3s stages, of which 20 477 were up-regulated and 13 102 were down-regulated. There were 32 199 differentially expressed genes between L3s and female stages, of these genes, 9 293 were up-regulated genes and 22 906 were down-regulated genes. The differentially expressed genes of two pairwise comparisons were respectively enriched in Gene Ontology. 6 617, 3 891 and 8 755 differentially expressed genes comparing eggs and L3s stages were annotated into database of biological process, cellular component and molecular function respectively, while the number by comparing L3s and female stages were 7 043, 3 686 and 10 177 respectively. In KEGG pathways identification, 6 521 differentially expressed genes comparing eggs and L3s stages were assigned to 251 KEGG pathways, and clustered significantly in MARK, Wnt signaling pathways and oxidative phosphorylation. In comparison of L3s and female stages, 6 528 differentially expressed genes were enriched in metabolic pathways, DNA replication and cell cycle. Functional cluster analysis indicated that the regulation related differential genes of growth rate and the reproductive and genital development were highly expressed in eggs and females stages, then the genes of the defense and carbohydrate metabolism were enriched in exclusive L3s stage. Differential genes of embryonic and post-embryonic development were highly expressed in all three stages, and 196 function genes from embryonic development and 166 function genes from post-embryonic development were co-expressed in all three stages, implied that these genes played a crucial role in (post-)embryonic development. Moreover, we identified 9 heterochronic genes, such as LIN-28, LIN-14 and RHEB-1, 48 nuclear hormone receptors (NHRs) genes, including NHR-49, NHR-48, NHR-40 and NHR-1, and 36 zinc metalloproteinase(NAS) genes, including NAS-36, NAS-33 and NAS-14. Then the analysis of enrichment capacity in three stages showed these genes were necessary to regulate the different development stages ofThe transcriptomic research ofat three developmental stages using RNA-seq revealed biological characteristics of differentially expressed genes involved in development- related GO functional classification, KEGG pathway and functional cluster and identified many kinds of heterochronic genes and developmental genes, which provides a foundation and reference for further investigation of the whole genome sequence analysis, interaction betweenand host, pathogenic mechanism and immune evasion.

; transcriptome; differentially expressed genes; development-related genes

2016-06-20；

2017-10-30

國家自然科學(xué)基金（31260603）

聯(lián)系方式：王文龍，E-mail：wwl.imau@163.com。馮陳晨，E-mail：410580576@qq.com。王文龍和馮陳晨為同等貢獻作者。通信作者呼和巴特爾，Tel：0471-4303726；E-mail：hhbte@163.com。通信作者劉春霞，Tel：0471-4309240；E-mail：lcx.imau@163.com