不同發育階段斯氏副柔線蟲比較轉錄組學分析

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

?

不同發育階段斯氏副柔線蟲比較轉錄組學分析

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

（1內蒙古農業大學獸醫學院/農業部動物疾病臨床診療技術重點實驗室，呼和浩特 010018；2內蒙古農業大學生命科學學院，呼和浩特 010018）

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

斯氏副柔線蟲；轉錄組；差異表達基因；發育相關基因

0 引言

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

1 材料與方法

試驗于2014—2016年在內蒙古農業大學獸醫學院完成。

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

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

1.2 方法

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

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

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

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

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

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

2 結果

2.1 斯氏副柔線蟲Unigenes聚類后數據可靠性分析

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

2.2 不同發育階段斯氏副柔線蟲的差異表達基因相關分析結果

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

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

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

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

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

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

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

1：下調差異表達基因 Up-regulated DEGs；2：上調差異表達基因 Down-regulated DEGs

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

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

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

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

2.5 差異表達基因功能聚類

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

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

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

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代謝通路中差異表達基因數量；KEGG：在KEGG數據庫中涉及到該通路的基因數

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

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

2.6 不同發育階段斯氏副柔線蟲生長發育相關基因分析

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

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

表4 蟲卵期、幼蟲期和雌蟲期生長發育相關的重要基因

3 討論

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

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

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

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

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

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

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

3.4 不同發育階段斯氏副柔線蟲生長發育相關基因分析

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

4 結論

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

[1] 宋銘忻, 張龍現. 獸醫寄生蟲學. 北京: 科學出版社, 2009: 183-184.

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

[2] 黃兵, 沈杰. 中國畜禽寄生蟲形態分類圖譜. 北京: 中國農業科學技術出版社, 2006: 462-463.

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

[3] 趙治國. 我國駱駝斯氏副柔線蟲病傳播媒介的研究[D]. 內蒙古農業大學, 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] 劉嬌, 張建珍, 李大琪, 張婷婷, 馬恩波, 張建琴. 中華稻蝗羧酸酯酶家族基因生物信息學及組織表達特異性分析. 中國農業科學, 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] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 張曌楠, 黃枳腱, 張璐, 王鴻權, 解彥玲, 童新宇. 中華蜜蜂幼蟲腸道響應球囊菌早期脅迫的轉錄組學. 中國農業科學, 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.

（責任編輯 林鑒非）

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

國家自然科學基金（31260603）

聯系方式：王文龍，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