南彥斌, 許嘉誠, 何啟玥, 潘學(xué)能, 周淵濤
(青海大學(xué)農(nóng)牧學(xué)院, 青海 西寧 810016)
微衛(wèi)星(Microsatellite)又名簡單重復(fù)序列(Single sequence reperts,SSR),它由1~6個核苷酸堿基串聯(lián)重復(fù)而形成,在真核生物基因組中分布比較廣泛[1]。微衛(wèi)星具有幾個非常顯著的特征,如可重復(fù)性高,共顯性,多態(tài)性高,數(shù)量豐富,對基因組有較好的覆蓋性等[2]。因此,其被廣泛應(yīng)用于群體遺傳結(jié)構(gòu)[3]、生物遺傳多樣性[4]、系統(tǒng)發(fā)生和遺傳圖譜繪制[5]等研究領(lǐng)域。雖然微衛(wèi)星標(biāo)記優(yōu)點突出但是無法掩蓋傳統(tǒng)的SSR引物開發(fā)方法成本極高且成功率低的缺點[6]。目前,高通量測序技術(shù)和生物信息學(xué)的的快速發(fā)展,使得SSR標(biāo)記開發(fā)更加高效便捷[7]。基于轉(zhuǎn)錄組數(shù)據(jù)已成功挖掘出多種昆蟲的SSR位點。譬如李敏等[8]從中華大仰蝽(Notonectachinensis) 37 801條unigenes中搜索到3124個SSR位點,利用半定量實驗篩選出16個SSR位點,他們認(rèn)為SSR新標(biāo)記的開發(fā)在中華大仰蝽的遺傳多樣性分析基因圖譜構(gòu)建中具有重要作用。王定鋒等[9]基于轉(zhuǎn)錄組數(shù)據(jù)挖掘灰茶尺蠖(Ectropisgrisescens)微衛(wèi)星標(biāo)記,發(fā)現(xiàn)灰茶尺蠖轉(zhuǎn)錄組中SSR位點數(shù)量大、出現(xiàn)頻率高、基元類型十分豐富,這些信息將為今后更好地利用SSR標(biāo)記研究灰茶尺蠖種群間的遺傳分化和制定灰茶尺蠖綜合防治措施在分子生物學(xué)方面提供借鑒。韓小紅[10]等通過對星天牛(Anoplophorachinensis)轉(zhuǎn)錄組SSR位點特征分析,發(fā)現(xiàn)其轉(zhuǎn)錄組中SSR的主要重復(fù)類型為單堿基重復(fù),其次是三堿基重復(fù)。SSR位點信息和位點特征是昆蟲分子標(biāo)記研究的重要內(nèi)容。然而截至目前,對于青海草原毛蟲微衛(wèi)星位點的研究還處于空白狀態(tài)。
草原毛蟲(Gynaephora)俗稱草原毒蛾,隸屬于鱗翅目毒蛾科。在世界范圍內(nèi),目前有記載的草原毛蟲屬昆蟲共有15種,而我國獨占8種且大部分分布在青藏高原上,對青藏高原高寒甸草產(chǎn)生了嚴(yán)重的危害[11]。該蟲通常一年發(fā)生一代,其1齡幼蟲在土塊、牛糞和草根下越冬,與國外的草原毛蟲生活史周期相比國內(nèi)的生活史周期較短[12-13]。由于青海草原毛蟲對生境的要求比較特殊,故而在國內(nèi)對它的研究較為有限,國外的相關(guān)研究也僅限于毒蛾科的其他蟲種[14-15]。近十年從分子生物學(xué)的角度出發(fā)關(guān)于青海草原毛蟲的研究,主要是從腸道菌、線粒體基因、序列變化和轉(zhuǎn)錄組基因表達(dá)等方面解析青海草原毛蟲適應(yīng)高海拔環(huán)境機制[15-18]。近年來有研究者從生物防治的角度開始探究青海草原毛蟲的防治。如楊忠岐等研究金小蜂對青海草原毛蟲的寄生作用[19],陳青紅等研究4種生物藥劑對青海草原毛蟲的防控效果[20],白海濤等利用昆蟲病原線蟲防治青海草原毛蟲21等。目前,隨著國家對生態(tài)環(huán)境的重視程度的增加,未來借助于分子生物學(xué)研究,推動青海草原毛蟲的生物防治是大勢所趨。本研究通過高通量測序技術(shù)獲得青海草原毛蟲的轉(zhuǎn)錄組數(shù)據(jù),結(jié)合生物信息學(xué)方法,對青海草原毛蟲的SSR位點進行挖掘并進行序列特征的統(tǒng)計與分析,設(shè)計和篩選能夠穩(wěn)定擴增的SSR引物,為青海草原毛蟲的分子鑒定與遺傳多樣性分析提供參考。
本試驗中的青海草原毛蟲樣品于2020年7月采自青海省海北州(36°59′ N,100°52′ E),設(shè)3個生物學(xué)重復(fù),每組3頭4齡幼蟲及雌雄成蟲,使用液氮快速冷凍后于—80℃超低溫冰箱中保存。取5只單頭4齡幼蟲用于SSR引物驗證,使用TIANamp Genomic DNA Kit基因組DNA試劑盒提取總DNA,經(jīng)過2%的瓊脂糖凝膠電泳和微量核酸分光光度計(北京凱奧科技發(fā)展有限公司)檢測,結(jié)果滿足進一步實驗要求,置于-20℃保存?zhèn)溆谩?/p>
委托北京諾禾致源科技股份有限公司通過Trizol法提取青海草原毛蟲樣品RNA,經(jīng)純度和濃度檢測合格后等物質(zhì)的量混合,采用Illumina HiSeqTM2000高通量測序平臺對青海草原毛蟲雌雄成蟲和4齡幼蟲進行全長轉(zhuǎn)錄組測序、文庫構(gòu)建和數(shù)據(jù)評估。在得到原始數(shù)據(jù)后通過SMRTlink軟件獲得未組裝的Subreads,進行校正,之后產(chǎn)生環(huán)形一致性序列(Circular-consensus sequence,CCS)。根據(jù)3′和5′引物及PloyA特征在CCS是否存在,從而確定全長序列和非全長序列,并分類找出全長非嵌合序列,通過同種類型聚類法(Iterative isoform-clustering,ICE)可將相似的全長非嵌合序列(Full length non-chimeric,FLNC)聚成一簇,得到Cluster,每個Cluster獲得一條一致性序列,結(jié)合非全長序列對其進行校對(Polishing),篩選得到高質(zhì)量序列用于后續(xù)分析[22]。
在Trinity v2.0.2拼接的基礎(chǔ)上,使用Corset v4.6軟件進行轉(zhuǎn)錄本聚類,去除冗余轉(zhuǎn)錄本,得到非冗余unigenes。對聚類得到的unigenes使用blast2 go v2.5軟件同GO,KOG和KEGG數(shù)據(jù)庫比對,進行基因功能注釋。
利用MISA軟件(Micro satellite identification tool,www.pgrc.ipk-gatresleben.de/misa)對青海草原毛蟲轉(zhuǎn)錄組中1 Kb以上的unigene進行SSR位點分析。對應(yīng)的各個unitsize的最少重復(fù)次數(shù)分別為10,6,5,5,5和5,重復(fù)單元的長度大小分別是1,2,3,4,5和6 bp。兩個相鄰的微衛(wèi)星重復(fù)單元間的間隔長度至少為100 bp。
在1.4節(jié)篩選到的SSR位點的基礎(chǔ)上,運用Primer3軟件[23]對引物進行批量設(shè)計。目標(biāo)擴增片段必須設(shè)置為包含SSR起始-3 bp,終止為+6 bp,擴增片段的長度為700~300 bp。引物的長度設(shè)為18~25 bp,最適長度為22 bp,引物最大能允許的不能識別的堿基數(shù)為1個,引物的最適退火溫度為60℃,退火溫度設(shè)置為55~65℃,上下游引物間的退火溫度差異不能超過3℃,引物末端穩(wěn)定性不能超過250[24]。
為了驗證引物的穩(wěn)定性,從批量設(shè)計的引物中隨機選出25對引物,交由生工生物工程(上海)股份有限公司。PCR反應(yīng)體系為25 μL:青海草原毛蟲幼蟲DNA模板1 μL,上下游引物(10 μmol·L-1)各加1 μL,2×EasyTaq?PCR SuperMix 12.5 μL,無核酸酶水9.5 μL。反應(yīng)程序:94℃預(yù)變性5 min;94℃變性30 s,57℃退火30 s,72℃延伸40 s,總共35個循環(huán);72℃延伸7 min,保存在4℃下。擴增后的產(chǎn)物用2%瓊脂糖凝膠電泳檢測(100V,30 min)。
采用Illumina HiSeqTM2000高通量測序平臺對青海草原毛蟲雌雄成蟲、4齡幼蟲轉(zhuǎn)錄組數(shù)據(jù)測序,獲得71 038 419條原始reads;測序獲得unigenes序列總長度為69 070 017 bp,堿基錯誤率為0.02%;Q20(質(zhì)量值≥20堿基所占百分比)含量平均約為98.18%,Q30(質(zhì)量值≥30堿基所占百分比)含量平均約為94.60%,測序質(zhì)量較好。采用Trinity軟件對clean reads數(shù)據(jù)進行拼接組裝,組裝后獲得63 335條all unigenes。其中,N90(拼接轉(zhuǎn)錄本按照長度從長到短排序,累加轉(zhuǎn)錄本的長度,到不小于總長90%的拼接轉(zhuǎn)錄本的長度)和N50(拼接轉(zhuǎn)錄本按照長度從長到短排序,累加轉(zhuǎn)錄本的長度,到不小于總長50%的拼接轉(zhuǎn)錄本的長度)的長度分別為440 bp和1 714 bp,平均長度為1 091 bp。所有63 335條unigenes序列用于后續(xù)的SSR搜索。
利用Blast2GO v2.5進行GO注釋,將注釋成功的基因按照GO三個大類(生物學(xué)過程,細(xì)胞組成,分子功能)的下一層級進行分類。結(jié)果表明(圖1),18 588條unigenes共得到78 461條功能注釋,分布在43個功能區(qū)。其中,獲得注釋序列數(shù)最多的是生物學(xué)過程分類,有39474條,分布在26個功能區(qū);其次是分子功能分類的注釋序列,有21 443條,分布在12個功能區(qū);注釋序列最少的是細(xì)胞組分類別,為17 544條,分布在5個功能區(qū)。在生物學(xué)過程類別中,細(xì)胞進程和代謝進程的注釋序列數(shù)占主導(dǎo)地位,分別為10 846和9 355條。在分子功能分類中,注釋到結(jié)合和催化活性的序列數(shù)最多,分別為9 932和7 295條。在細(xì)胞組分分類中,注釋到細(xì)胞結(jié)構(gòu)體和胞內(nèi)部分的序列占主導(dǎo)地位,分別為8 267和4 877條。有11個功能區(qū)的注釋序列超過2 000條,分別是細(xì)胞進程(10 846),結(jié)合(9 932),代謝過程(9 355),細(xì)胞結(jié)構(gòu)體(8 267),催化活性(7 295),生物調(diào)節(jié)(4 085),生物調(diào)控進程(3 803),含蛋白復(fù)合物(3 490),定位(3 002)和應(yīng)激響應(yīng)(2 607)。
圖1 青海草原毛蟲unigenes GO功能注釋
KEGG代謝通路富集分析結(jié)果表明(圖2),共有11 846條unigenes參與到細(xì)胞進程、環(huán)境信息處理、遺傳信息處理、新陳代謝和有機體系統(tǒng)這五大類生化代謝通路中。其中代謝系統(tǒng)途徑中基因最多(3 346條),它們主要參與碳水化合物代謝和脂質(zhì)代謝等過程,分別為566條和436條。在這34組代謝通路子類別中,信號轉(zhuǎn)導(dǎo)(Signal transduction)通路的基因最多,為1 302條。這些代謝通路相關(guān)的基因分布于已知的284個代謝通路,其中富集最多的10條通路分別是核糖體,內(nèi)質(zhì)網(wǎng)蛋白加工,嘌呤代謝,RNA轉(zhuǎn)運,剪接體,氧化磷酸化,碳代謝,嘧啶代謝,cAMP信號通路,真核生物的核糖體生物發(fā)生,真核生物中的核糖體生物合成(表1)。
表1 KEGG數(shù)據(jù)庫中單基因的十大代謝途徑富集
對有注釋的基因進行直系同源分類(圖3),9 521條unigenes共得到10 628條注釋,分為26類,其中,一般功能預(yù)測和翻譯后修飾的基因最多,分別為1 487條(13.99%) 和1 114條(10.48%);其余依次是T類信號轉(zhuǎn)導(dǎo)機制,有1 005條(9.46%);J類翻譯、核糖體結(jié)構(gòu)和生物合成,有864條(8.13%);K類轉(zhuǎn)錄,有644條(6.06%);U類細(xì)胞內(nèi)運輸、分泌和囊泡轉(zhuǎn)運,有639條(6.01%);S類功能未知,有630條(5.93%);A類核酸處理與修飾,有598條(5.63%)。其余的都低于500條,X類未命名蛋白質(zhì),有1條(0.01%)。
圖3 青海草原毛蟲unigenes KOG分類注釋
對青海草原毛蟲轉(zhuǎn)錄組中的63 335條unigenes的數(shù)據(jù)利用MISA軟件進行簡單重復(fù)序列(SSR)位點搜索??偣舱业?2 597個微衛(wèi)星位點,分布于9 851條unigenes中,其微衛(wèi)星出現(xiàn)頻率(含有SSR的unigenes數(shù)量與總unigenes數(shù)量之比)是15.55%,其中含有1個以上SSR位點的獨立基因序列有1 852條,以復(fù)合型形式存在的SSR有801個,SSR的分布頻率(SSR個數(shù)與總unigene數(shù)量之比)為19.89%(表2)。這些SSR基序包含1~5 bp的串聯(lián)重復(fù)序列。
表2 青海草原毛蟲轉(zhuǎn)錄組SSR分布及特征
從青海草原毛蟲轉(zhuǎn)錄組數(shù)據(jù)中可以得到,單核苷酸重復(fù)是最主要的重復(fù),占SSR總數(shù)的71.83%。二核苷酸重復(fù)居其次,占SSR總數(shù)的14.27%。三、四核苷酸重復(fù),分別占SSR總數(shù)的10.59%和2.99%。五核苷酸和六核苷酸重復(fù)相較于前幾種核苷酸的含量很少,分別占SSR總數(shù)的0.29%和0.0.3%。出現(xiàn)頻率最高的是10次重復(fù)基元,共5 265個SSR位點,占總SSR位點的41.77%(表3)。
表3 青海草原毛蟲轉(zhuǎn)錄組數(shù)據(jù)的SSR重復(fù)類型分布
對5~99次SSR基元重復(fù)次數(shù)進行分析結(jié)果表明(圖4),SSR位點重復(fù)在5~10次的有8 335個,占總SSR位點的66.17%;重復(fù)次數(shù)在11~20次的有3 896個,占總SSR位點的30.93%;重復(fù)次數(shù)在21~30次的有267個,占總SSR位點的2.12%;重復(fù)次數(shù)在31~40次的有53個,占總SSR位點的0.42%;重復(fù)次數(shù)在41~50次的有20個,占總SSR位點的0.16%;重復(fù)次數(shù)在51~99次的有26個,占總SSR位點的0.21%。對單核苷酸重復(fù),二核苷酸重復(fù)和三核苷酸重復(fù)的分析發(fā)現(xiàn),SSR的數(shù)量與重復(fù)次數(shù)成反比出現(xiàn),隨著重復(fù)次數(shù)的增加,SSR的數(shù)量反而在減少。在圖中未呈現(xiàn)的四、五、六核苷酸重復(fù)具有與單、二、三核苷酸重復(fù)相同的趨勢表現(xiàn)。
圖4 青海草原毛蟲轉(zhuǎn)錄組中SSR基元重復(fù)次數(shù)
青海草原毛蟲轉(zhuǎn)錄組SSR基元類型結(jié)果顯示(圖5和表4),有60個SSR重復(fù)基元類型,共獲得12 597個SSR位點。其中,2種類型的基元在單核苷酸重復(fù)中呈現(xiàn),(A/T)n是最主要的基元類型,形成8990個SSR位點,占同類基元SSR位點的99.36%,占總SSR位點的71.37%。4種類型的基元在二核苷酸重復(fù)中呈現(xiàn),(AC/GT)n為主要優(yōu)勢基元,形成763個SSR位點(42.46%,6.06%),其次為(AT/AT)n和(AG/CT)n基元,分別形成619個(34.45%,4.91%)和370個SSR位點(20.59%,2.94%)。10種類型的基元在三核苷酸重復(fù)中呈現(xiàn),(AAT/ATT)n是最主要的優(yōu)勢基元,有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種類型的基元在四核苷酸重復(fù)中呈現(xiàn),(AAAT/ATTT)n是最主要的優(yōu)勢基元,有170個SSR位點(45.09%,1.35%);其他優(yōu)勢基元分別為(ACAT/ATGT)n、(AGAT/ATCT)n,分別形成57個(15.12%,0.45%)和44個(11.67%,0.35%)SSR位點。在五核苷酸重復(fù)中有22個重復(fù)基元類型,最主要的優(yōu)勢基元是(AAAAT/ATTTT)n,(AACGT/ACGTT)n,在這兩個基元中都形成了4個(11.11%,0.03%)SSR位點。六核苷酸重復(fù)有5個重復(fù)基元類型,全都形成1個SSR位點(20%,0.01%)。
表4 青海草原毛蟲轉(zhuǎn)錄組不同重復(fù)基元SSR位點數(shù)及頻率
圖5 青海草原毛蟲轉(zhuǎn)錄組不同重復(fù)類型優(yōu)勢基元SSR位點數(shù)分布
本研究從被篩選出的微衛(wèi)星引物中隨機選取25對,對5個不同樣品(海北州同一地點的5個青海草原毛蟲個體重復(fù))的青海草原毛蟲DNA樣本進行PCR擴增,結(jié)果11對引物在每個樣品中擴增取得的片段與預(yù)期的片段幾乎一致。有效擴增率為44%(表5,圖6)。
表5 青海草原毛蟲部分EST-SSR篩選引物信息
圖6 青海草原毛蟲5個樣品25個SSR位點引物的PCR擴增產(chǎn)物電泳結(jié)果
伴隨著高通量測序技術(shù)的快速發(fā)展,對昆蟲轉(zhuǎn)錄組學(xué)的研究變得非常廣泛[25]。本研究測序獲得unigenes序列總長度為69 070 017 bp,平均長度為1 091 bp,N50長度為1 714 bp,N90長度為440 bp。由于N50的值較大,因此其組裝完整性較高,保證了轉(zhuǎn)錄組分析的準(zhǔn)確性。
本研究發(fā)現(xiàn),青海草原毛蟲轉(zhuǎn)錄組數(shù)據(jù)庫中存在大量的微衛(wèi)星,且種類豐富,核苷酸的重復(fù)類型出現(xiàn)了6種。單核苷酸重復(fù)中最主要的重復(fù)類型是A/T,這與對印度谷螟(Plodiainterpunctella)[2]、黏蟲(Mythimnaseparata)[26]、褐飛虱(Nilaparuatalugens)[27]SSR位點特征的研究結(jié)果基本相同。但本次試驗結(jié)果與對黃粉蟲(Tenebriomolitor)[28]、麥紅吸漿蟲(Sitodiplosismosellana)[29]、沙蔥葉甲(Galerucadaurica)[30]等昆蟲的SSR位點特征研究結(jié)果相比較存在差異。這些研究發(fā)現(xiàn)核苷酸的重復(fù)類型以單核苷酸為主,三核苷酸次之。本研究的結(jié)果是單核苷酸重復(fù)是最主要的重復(fù)類型,二核苷酸次之。該結(jié)果與意大利微蝗(Calliptamusitalicus)[31]、梨小食心蟲(Grapholithamolesta)[32]、沙棘木蠹蛾(Holcocerushippophaecolus)的研究結(jié)果類似。通常情況下,在ESTs里,單核苷酸為主要重復(fù)類型,三核苷酸重復(fù)居于其次。原因在于,與編碼區(qū)其他的重復(fù)基元類型相比,三核苷酸核心基元更加穩(wěn)定,基本不會出現(xiàn)編碼框滑動突變現(xiàn)象[33]。本研究中出現(xiàn)差異,可能是不同樣品之間本身存在差異,從而導(dǎo)致結(jié)果不同。
為驗證青海草原毛蟲SSR的穩(wěn)定性,本研究利用其轉(zhuǎn)錄組數(shù)據(jù)設(shè)計出2 714對引物,從中隨機選出25對引物進行擴增試驗,結(jié)果表明11對引物在每個單個樣品中都能成功擴增,并得到預(yù)期條帶。這些位點的多態(tài)性目前還不明確,需要進一步試驗評估,但它的數(shù)量足以開展青海草原毛蟲種群遺傳學(xué)相關(guān)研究。導(dǎo)致不同樣品的EST-SSR位點擴增失敗的原因有很多,如:由于基因組DNA中的內(nèi)含子在轉(zhuǎn)錄后沒表現(xiàn)出來,對擴增過程造成干擾從而導(dǎo)致預(yù)期的條帶與實際擴增的條帶相比,會出現(xiàn)擴增條帶遠(yuǎn)遠(yuǎn)大于預(yù)期的現(xiàn)象;在基因組中EST-SSR的表達(dá)豐富度太低,導(dǎo)致產(chǎn)物的濃度低,無法檢測[2]。
微衛(wèi)星的重要作用體現(xiàn)在基因調(diào)控、蛋白功能、種群遺傳結(jié)構(gòu)、種群擴散及遷移等方面[34]。引物的開發(fā)是微衛(wèi)星應(yīng)用的前提,傳統(tǒng)的方法開發(fā)成本太高,不合時宜。通過在轉(zhuǎn)錄組高通量數(shù)據(jù)中直接發(fā)掘微衛(wèi)星并設(shè)計引物的方法,大大減少了開發(fā)成本,很大程度上對微衛(wèi)星的開發(fā)進程起到推動作用[30]。本研究在青海草原毛蟲轉(zhuǎn)錄組數(shù)據(jù)庫中共篩出12 597個SSR位點,并對其分布情況及可用性進行了詳細(xì)敘述,這些SSR的發(fā)掘,對于之后研究青海草原毛蟲的遺傳多樣性和功能基因組學(xué)奠定了堅實的基礎(chǔ),也更進一步表明通過高通量測序技術(shù)發(fā)掘昆蟲的SSR引物是一種符合人類需求的研究方法。
本研究通過KOG,GO注釋和KEGG通路三大數(shù)據(jù)庫分析,發(fā)現(xiàn)unigenes注釋主要集中于一般功能預(yù)測、細(xì)胞進程和碳水化合物代謝。本研究共篩出12 597個SSR位點,發(fā)生頻率為19.89%,SSR位點數(shù)量豐富,(A/T)n,(AC/GT)n分別是單核苷酸重復(fù)和二核苷酸的優(yōu)勢基元,SSR數(shù)量隨重復(fù)次數(shù)增加而降低,重復(fù)次數(shù)類型隨基元序列長度增長而減少。從1 556條青海草原毛蟲unigenes中成功設(shè)計出2 714對青海草原毛蟲SSR引物,隨機挑出25對引物進行PCR驗證,有11對引物擴增出目的DNA片段。本研究基于轉(zhuǎn)錄組數(shù)據(jù)成功篩選出青海草原毛蟲微衛(wèi)星位點,獲得的青海草原毛蟲SSR位點為進一步研究其種群遺傳學(xué)和種群擴散及遷移奠定基礎(chǔ)。