青海草原毛蟲轉錄組分析及SSR位點開發

南彥斌, 許嘉誠, 何啟玥, 潘學能, 周淵濤

(青海大學農牧學院, 青海 西寧 810016)

微衛星(Microsatellite)又名簡單重復序列(Single sequence reperts,SSR),它由1～6個核苷酸堿基串聯重復而形成,在真核生物基因組中分布比較廣泛[1]。微衛星具有幾個非常顯著的特征,如可重復性高,共顯性,多態性高,數量豐富,對基因組有較好的覆蓋性等[2]。因此,其被廣泛應用于群體遺傳結構[3]、生物遺傳多樣性[4]、系統發生和遺傳圖譜繪制[5]等研究領域。雖然微衛星標記優點突出但是無法掩蓋傳統的SSR引物開發方法成本極高且成功率低的缺點[6]。目前,高通量測序技術和生物信息學的的快速發展,使得SSR標記開發更加高效便捷[7]。基于轉錄組數據已成功挖掘出多種昆蟲的SSR位點。譬如李敏等[8]從中華大仰蝽(Notonectachinensis) 37 801條unigenes中搜索到3124個SSR位點,利用半定量實驗篩選出16個SSR位點,他們認為SSR新標記的開發在中華大仰蝽的遺傳多樣性分析基因圖譜構建中具有重要作用。王定鋒等[9]基于轉錄組數據挖掘灰茶尺蠖(Ectropisgrisescens)微衛星標記,發現灰茶尺蠖轉錄組中SSR位點數量大、出現頻率高、基元類型十分豐富,這些信息將為今后更好地利用SSR標記研究灰茶尺蠖種群間的遺傳分化和制定灰茶尺蠖綜合防治措施在分子生物學方面提供借鑒。韓小紅[10]等通過對星天牛(Anoplophorachinensis)轉錄組SSR位點特征分析,發現其轉錄組中SSR的主要重復類型為單堿基重復,其次是三堿基重復。SSR位點信息和位點特征是昆蟲分子標記研究的重要內容。然而截至目前,對于青海草原毛蟲微衛星位點的研究還處于空白狀態。

草原毛蟲(Gynaephora)俗稱草原毒蛾,隸屬于鱗翅目毒蛾科。在世界范圍內,目前有記載的草原毛蟲屬昆蟲共有15種,而我國獨占8種且大部分分布在青藏高原上,對青藏高原高寒甸草產生了嚴重的危害[11]。該蟲通常一年發生一代,其1齡幼蟲在土塊、牛糞和草根下越冬,與國外的草原毛蟲生活史周期相比國內的生活史周期較短[12-13]。由于青海草原毛蟲對生境的要求比較特殊,故而在國內對它的研究較為有限,國外的相關研究也僅限于毒蛾科的其他蟲種[14-15]。近十年從分子生物學的角度出發關于青海草原毛蟲的研究,主要是從腸道菌、線粒體基因、序列變化和轉錄組基因表達等方面解析青海草原毛蟲適應高海拔環境機制[15-18]。近年來有研究者從生物防治的角度開始探究青海草原毛蟲的防治。如楊忠岐等研究金小蜂對青海草原毛蟲的寄生作用[19],陳青紅等研究4種生物藥劑對青海草原毛蟲的防控效果[20],白海濤等利用昆蟲病原線蟲防治青海草原毛蟲21等。目前,隨著國家對生態環境的重視程度的增加,未來借助于分子生物學研究,推動青海草原毛蟲的生物防治是大勢所趨。本研究通過高通量測序技術獲得青海草原毛蟲的轉錄組數據,結合生物信息學方法,對青海草原毛蟲的SSR位點進行挖掘并進行序列特征的統計與分析,設計和篩選能夠穩定擴增的SSR引物,為青海草原毛蟲的分子鑒定與遺傳多樣性分析提供參考。

1 材料與方法

1.1 供試昆蟲

本試驗中的青海草原毛蟲樣品于2020年7月采自青海省海北州(36°59′ N,100°52′ E),設3個生物學重復,每組3頭4齡幼蟲及雌雄成蟲,使用液氮快速冷凍后于—80℃超低溫冰箱中保存。取5只單頭4齡幼蟲用于SSR引物驗證,使用TIANamp Genomic DNA Kit基因組DNA試劑盒提取總DNA,經過2%的瓊脂糖凝膠電泳和微量核酸分光光度計(北京凱奧科技發展有限公司)檢測,結果滿足進一步實驗要求,置于-20℃保存備用。

1.2 總RNA提取、建庫及數據評估

委托北京諾禾致源科技股份有限公司通過Trizol法提取青海草原毛蟲樣品RNA,經純度和濃度檢測合格后等物質的量混合,采用Illumina HiSeqTM2000高通量測序平臺對青海草原毛蟲雌雄成蟲和4齡幼蟲進行全長轉錄組測序、文庫構建和數據評估。在得到原始數據后通過SMRTlink軟件獲得未組裝的Subreads,進行校正,之后產生環形一致性序列(Circular-consensus sequence,CCS)。根據3′和5′引物及PloyA特征在CCS是否存在,從而確定全長序列和非全長序列,并分類找出全長非嵌合序列,通過同種類型聚類法(Iterative isoform-clustering,ICE)可將相似的全長非嵌合序列(Full length non-chimeric,FLNC)聚成一簇,得到Cluster,每個Cluster獲得一條一致性序列,結合非全長序列對其進行校對(Polishing),篩選得到高質量序列用于后續分析[22]。

1.3 青海草原毛蟲轉錄組unigenes的功能注釋

在Trinity v2.0.2拼接的基礎上,使用Corset v4.6軟件進行轉錄本聚類,去除冗余轉錄本,得到非冗余unigenes。對聚類得到的unigenes使用blast2 go v2.5軟件同GO,KOG和KEGG數據庫比對,進行基因功能注釋。

1.4 SSR的篩選

利用MISA軟件(Micro satellite identification tool,www.pgrc.ipk-gatresleben.de/misa)對青海草原毛蟲轉錄組中1 Kb以上的unigene進行SSR位點分析。對應的各個unitsize的最少重復次數分別為10,6,5,5,5和5,重復單元的長度大小分別是1,2,3,4,5和6 bp。兩個相鄰的微衛星重復單元間的間隔長度至少為100 bp。

1.5 SSR引物設計與驗證

在1.4節篩選到的SSR位點的基礎上,運用Primer3軟件[23]對引物進行批量設計。目標擴增片段必須設置為包含SSR起始-3 bp,終止為+6 bp,擴增片段的長度為700～300 bp。引物的長度設為18～25 bp,最適長度為22 bp,引物最大能允許的不能識別的堿基數為1個,引物的最適退火溫度為60℃,退火溫度設置為55～65℃,上下游引物間的退火溫度差異不能超過3℃,引物末端穩定性不能超過250[24]。

為了驗證引物的穩定性,從批量設計的引物中隨機選出25對引物,交由生工生物工程(上海)股份有限公司。PCR反應體系為25 μL:青海草原毛蟲幼蟲DNA模板1 μL,上下游引物(10 μmol·L-1)各加1 μL,2×EasyTaq?PCR SuperMix 12.5 μL,無核酸酶水9.5 μL。反應程序:94℃預變性5 min;94℃變性30 s,57℃退火30 s,72℃延伸40 s,總共35個循環;72℃延伸7 min,保存在4℃下。擴增后的產物用2%瓊脂糖凝膠電泳檢測(100V,30 min)。

2 結果與分析

2.1 青海草原毛蟲轉錄組測序分析

采用Illumina HiSeqTM2000高通量測序平臺對青海草原毛蟲雌雄成蟲、4齡幼蟲轉錄組數據測序,獲得71 038 419條原始reads;測序獲得unigenes序列總長度為69 070 017 bp,堿基錯誤率為0.02%;Q20(質量值≥20堿基所占百分比)含量平均約為98.18%,Q30(質量值≥30堿基所占百分比)含量平均約為94.60%,測序質量較好。采用Trinity軟件對clean reads數據進行拼接組裝,組裝后獲得63 335條all unigenes。其中,N90(拼接轉錄本按照長度從長到短排序,累加轉錄本的長度,到不小于總長90%的拼接轉錄本的長度)和N50(拼接轉錄本按照長度從長到短排序,累加轉錄本的長度,到不小于總長50%的拼接轉錄本的長度)的長度分別為440 bp和1 714 bp,平均長度為1 091 bp。所有63 335條unigenes序列用于后續的SSR搜索。

2.2 青海草原毛蟲轉錄組的GO注釋

利用Blast2GO v2.5進行GO注釋,將注釋成功的基因按照GO三個大類(生物學過程,細胞組成,分子功能)的下一層級進行分類。結果表明(圖1),18 588條unigenes共得到78 461條功能注釋,分布在43個功能區。其中,獲得注釋序列數最多的是生物學過程分類,有39474條,分布在26個功能區;其次是分子功能分類的注釋序列,有21 443條,分布在12個功能區;注釋序列最少的是細胞組分類別,為17 544條,分布在5個功能區。在生物學過程類別中,細胞進程和代謝進程的注釋序列數占主導地位,分別為10 846和9 355條。在分子功能分類中,注釋到結合和催化活性的序列數最多,分別為9 932和7 295條。在細胞組分分類中,注釋到細胞結構體和胞內部分的序列占主導地位,分別為8 267和4 877條。有11個功能區的注釋序列超過2 000條,分別是細胞進程(10 846),結合(9 932),代謝過程(9 355),細胞結構體(8 267),催化活性(7 295),生物調節(4 085),生物調控進程(3 803),含蛋白復合物(3 490),定位(3 002)和應激響應(2 607)。

圖1 青海草原毛蟲unigenes GO功能注釋

2.3 青海草原毛蟲Unigenes的KEGG分類

KEGG代謝通路富集分析結果表明(圖2),共有11 846條unigenes參與到細胞進程、環境信息處理、遺傳信息處理、新陳代謝和有機體系統這五大類生化代謝通路中。其中代謝系統途徑中基因最多(3 346條),它們主要參與碳水化合物代謝和脂質代謝等過程,分別為566條和436條。在這34組代謝通路子類別中,信號轉導(Signal transduction)通路的基因最多,為1 302條。這些代謝通路相關的基因分布于已知的284個代謝通路,其中富集最多的10條通路分別是核糖體,內質網蛋白加工,嘌呤代謝,RNA轉運,剪接體,氧化磷酸化,碳代謝,嘧啶代謝,cAMP信號通路,真核生物的核糖體生物發生,真核生物中的核糖體生物合成(表1)。

表1 KEGG數據庫中單基因的十大代謝途徑富集

2.4 青海草原毛蟲Unigenes的KOG分類

對有注釋的基因進行直系同源分類(圖3),9 521條unigenes共得到10 628條注釋,分為26類,其中,一般功能預測和翻譯后修飾的基因最多,分別為1 487條(13.99%) 和1 114條(10.48%);其余依次是T類信號轉導機制,有1 005條(9.46%);J類翻譯、核糖體結構和生物合成,有864條(8.13%);K類轉錄,有644條(6.06%);U類細胞內運輸、分泌和囊泡轉運,有639條(6.01%);S類功能未知,有630條(5.93%);A類核酸處理與修飾,有598條(5.63%)。其余的都低于500條,X類未命名蛋白質,有1條(0.01%)。

圖3 青海草原毛蟲unigenes KOG分類注釋

2.5 青海草原毛蟲轉錄組中SSR分布特征

對青海草原毛蟲轉錄組中的63 335條unigenes的數據利用MISA軟件進行簡單重復序列(SSR)位點搜索。總共找到12 597個微衛星位點,分布于9 851條unigenes中,其微衛星出現頻率(含有SSR的unigenes數量與總unigenes數量之比)是15.55%,其中含有1個以上SSR位點的獨立基因序列有1 852條,以復合型形式存在的SSR有801個,SSR的分布頻率(SSR個數與總unigene數量之比)為19.89%(表2)。這些SSR基序包含1～5 bp的串聯重復序列。

表2 青海草原毛蟲轉錄組SSR分布及特征

2.6 青海草原毛蟲轉錄組序列中SSR位點重復類型分析

從青海草原毛蟲轉錄組數據中可以得到,單核苷酸重復是最主要的重復,占SSR總數的71.83%。二核苷酸重復居其次,占SSR總數的14.27%。三、四核苷酸重復,分別占SSR總數的10.59%和2.99%。五核苷酸和六核苷酸重復相較于前幾種核苷酸的含量很少,分別占SSR總數的0.29%和0.0.3%。出現頻率最高的是10次重復基元,共5 265個SSR位點,占總SSR位點的41.77%(表3)。

表3 青海草原毛蟲轉錄組數據的SSR重復類型分布

對5～99次SSR基元重復次數進行分析結果表明(圖4),SSR位點重復在5～10次的有8 335個,占總SSR位點的66.17%;重復次數在11～20次的有3 896個,占總SSR位點的30.93%;重復次數在21～30次的有267個,占總SSR位點的2.12%;重復次數在31～40次的有53個,占總SSR位點的0.42%;重復次數在41～50次的有20個,占總SSR位點的0.16%;重復次數在51～99次的有26個,占總SSR位點的0.21%。對單核苷酸重復,二核苷酸重復和三核苷酸重復的分析發現,SSR的數量與重復次數成反比出現,隨著重復次數的增加,SSR的數量反而在減少。在圖中未呈現的四、五、六核苷酸重復具有與單、二、三核苷酸重復相同的趨勢表現。

圖4 青海草原毛蟲轉錄組中SSR基元重復次數

2.7 青海草原毛蟲轉錄組基元類型分析

青海草原毛蟲轉錄組SSR基元類型結果顯示(圖5和表4),有60個SSR重復基元類型,共獲得12 597個SSR位點。其中,2種類型的基元在單核苷酸重復中呈現,(A/T)n是最主要的基元類型,形成8990個SSR位點,占同類基元SSR位點的99.36%,占總SSR位點的71.37%。4種類型的基元在二核苷酸重復中呈現,(AC/GT)n為主要優勢基元,形成763個SSR位點(42.46%,6.06%),其次為(AT/AT)n和(AG/CT)n基元,分別形成619個(34.45%,4.91%)和370個SSR位點(20.59%,2.94%)。10種類型的基元在三核苷酸重復中呈現,(AAT/ATT)n是最主要的優勢基元,有533個SSR位點(39.96%,4.23%);其次為(ATC/ATG)n,(AAC/GTT)n,(AAG/CTT)n。分別形成236個(17.69%,1.87%),142個(10.64%,1.13%),114個(8.55%,0.91%)SSR位點。16種類型的基元在四核苷酸重復中呈現,(AAAT/ATTT)n是最主要的優勢基元,有170個SSR位點(45.09%,1.35%);其他優勢基元分別為(ACAT/ATGT)n、(AGAT/ATCT)n,分別形成57個(15.12%,0.45%)和44個(11.67%,0.35%)SSR位點。在五核苷酸重復中有22個重復基元類型,最主要的優勢基元是(AAAAT/ATTTT)n,(AACGT/ACGTT)n,在這兩個基元中都形成了4個(11.11%,0.03%)SSR位點。六核苷酸重復有5個重復基元類型,全都形成1個SSR位點(20%,0.01%)。

表4 青海草原毛蟲轉錄組不同重復基元SSR位點數及頻率

圖5 青海草原毛蟲轉錄組不同重復類型優勢基元SSR位點數分布

2.8 青海草原毛蟲SSRs引物驗證

本研究從被篩選出的微衛星引物中隨機選取25對,對5個不同樣品(海北州同一地點的5個青海草原毛蟲個體重復)的青海草原毛蟲DNA樣本進行PCR擴增,結果11對引物在每個樣品中擴增取得的片段與預期的片段幾乎一致。有效擴增率為44%(表5,圖6)。

表5 青海草原毛蟲部分EST-SSR篩選引物信息

圖6 青海草原毛蟲5個樣品25個SSR位點引物的PCR擴增產物電泳結果

3 討論

伴隨著高通量測序技術的快速發展,對昆蟲轉錄組學的研究變得非常廣泛[25]。本研究測序獲得unigenes序列總長度為69 070 017 bp,平均長度為1 091 bp,N50長度為1 714 bp,N90長度為440 bp。由于N50的值較大,因此其組裝完整性較高,保證了轉錄組分析的準確性。

本研究發現,青海草原毛蟲轉錄組數據庫中存在大量的微衛星,且種類豐富,核苷酸的重復類型出現了6種。單核苷酸重復中最主要的重復類型是A/T,這與對印度谷螟(Plodiainterpunctella)[2]、黏蟲(Mythimnaseparata)[26]、褐飛虱(Nilaparuatalugens)[27]SSR位點特征的研究結果基本相同。但本次試驗結果與對黃粉蟲(Tenebriomolitor)[28]、麥紅吸漿蟲(Sitodiplosismosellana)[29]、沙蔥葉甲(Galerucadaurica)[30]等昆蟲的SSR位點特征研究結果相比較存在差異。這些研究發現核苷酸的重復類型以單核苷酸為主,三核苷酸次之。本研究的結果是單核苷酸重復是最主要的重復類型,二核苷酸次之。該結果與意大利微蝗(Calliptamusitalicus)[31]、梨小食心蟲(Grapholithamolesta)[32]、沙棘木蠹蛾(Holcocerushippophaecolus)的研究結果類似。通常情況下,在ESTs里,單核苷酸為主要重復類型,三核苷酸重復居于其次。原因在于,與編碼區其他的重復基元類型相比,三核苷酸核心基元更加穩定,基本不會出現編碼框滑動突變現象[33]。本研究中出現差異,可能是不同樣品之間本身存在差異,從而導致結果不同。

為驗證青海草原毛蟲SSR的穩定性,本研究利用其轉錄組數據設計出2 714對引物,從中隨機選出25對引物進行擴增試驗,結果表明11對引物在每個單個樣品中都能成功擴增,并得到預期條帶。這些位點的多態性目前還不明確,需要進一步試驗評估,但它的數量足以開展青海草原毛蟲種群遺傳學相關研究。導致不同樣品的EST-SSR位點擴增失敗的原因有很多,如:由于基因組DNA中的內含子在轉錄后沒表現出來,對擴增過程造成干擾從而導致預期的條帶與實際擴增的條帶相比,會出現擴增條帶遠遠大于預期的現象;在基因組中EST-SSR的表達豐富度太低,導致產物的濃度低,無法檢測[2]。

微衛星的重要作用體現在基因調控、蛋白功能、種群遺傳結構、種群擴散及遷移等方面[34]。引物的開發是微衛星應用的前提,傳統的方法開發成本太高,不合時宜。通過在轉錄組高通量數據中直接發掘微衛星并設計引物的方法,大大減少了開發成本,很大程度上對微衛星的開發進程起到推動作用[30]。本研究在青海草原毛蟲轉錄組數據庫中共篩出12 597個SSR位點,并對其分布情況及可用性進行了詳細敘述,這些SSR的發掘,對于之后研究青海草原毛蟲的遺傳多樣性和功能基因組學奠定了堅實的基礎,也更進一步表明通過高通量測序技術發掘昆蟲的SSR引物是一種符合人類需求的研究方法。

4 結論

本研究通過KOG,GO注釋和KEGG通路三大數據庫分析,發現unigenes注釋主要集中于一般功能預測、細胞進程和碳水化合物代謝。本研究共篩出12 597個SSR位點,發生頻率為19.89%,SSR位點數量豐富,(A/T)n,(AC/GT)n分別是單核苷酸重復和二核苷酸的優勢基元,SSR數量隨重復次數增加而降低,重復次數類型隨基元序列長度增長而減少。從1 556條青海草原毛蟲unigenes中成功設計出2 714對青海草原毛蟲SSR引物,隨機挑出25對引物進行PCR驗證,有11對引物擴增出目的DNA片段。本研究基于轉錄組數據成功篩選出青海草原毛蟲微衛星位點,獲得的青海草原毛蟲SSR位點為進一步研究其種群遺傳學和種群擴散及遷移奠定基礎。