蔣月麗,李 彤,武予清,苗 進(jìn),鞏中軍,段 云,劉啟航
(河南省農(nóng)業(yè)科學(xué)院植物保護(hù)研究所,河南省農(nóng)作物病蟲害防治重點(diǎn)實(shí)驗(yàn)室,農(nóng)業(yè)部華北南部有害生物治理重點(diǎn)實(shí)驗(yàn)室,鄭州 450002)
麥紅吸漿蟲Sitodiplosismosellana(Gehin)屬雙翅目癭蚊科Cecidomyiidae,是一種世界性的小麥害蟲,主要分布在美國、加拿大、歐洲、俄羅斯和中國(Lambetal., 2000; Berzonskyetal., 2003),以幼蟲潛伏在穎殼內(nèi)吸食正在灌漿的汁液,造成麥粒癟瘡、空殼或霉?fàn)€而減產(chǎn),具有很大的危害性,一般減產(chǎn)10%~20%,重者減產(chǎn)30%~50%,嚴(yán)重的乃至顆粒無收,全球每年因麥紅吸漿蟲的危害而造成的經(jīng)濟(jì)損失高達(dá)數(shù)億美元(Doane and Olfert, 2008; Oakley, 2008)。由于麥紅吸漿蟲具有蟲體小,在土壤中滯育時(shí)間長,在田間呈島嶼型分布和發(fā)生危害具有隱蔽性、間歇性和周期性的特點(diǎn)(Wise and Lamb, 2004; Liatukasetal., 2009; 郁振興等,2011),長期以來一直是威脅小麥生產(chǎn)的重要害蟲。目前,關(guān)于吸漿蟲的研究大多集中在生態(tài)和生物化學(xué)(Hinks and Doane, 1988; Hu and Zhang, 1995; Wise and Lamb, 2004; Wu and Yuan, 2004; Chenetal., 2009),但是由于基因組和轉(zhuǎn)錄組信息的缺乏,很少有關(guān)于其功能基因的研究,本文通過獲得大量轉(zhuǎn)錄組信息,可以為其功能基因的挖掘提供巨大的信息資源。
轉(zhuǎn)錄組(transcriptome)是指細(xì)胞在特定狀態(tài)下全部表達(dá)的RNA的總和,反映相同基因在不同條件下表達(dá)水平的差異,并能揭示不同基因的相互作用及各自功能(Velculescu, 1997)。轉(zhuǎn)錄組學(xué)是一門基于轉(zhuǎn)錄組在整體水平之上研究生物基因轉(zhuǎn)錄情況的學(xué)科。近年來,轉(zhuǎn)錄組學(xué)技術(shù)在揭示細(xì)胞生理活動(dòng)規(guī)律和代謝機(jī)理的研究中廣泛應(yīng)用(Lockhart & Winzeler, 2008)。轉(zhuǎn)錄組測序(RNA sequencing)分為三代。第一代測序技術(shù)以基因表達(dá)系列分析技術(shù)(Serial Analysis of GeneExpression,SAGE)和大規(guī)模平行測序技術(shù)(Massively Paralel Signature Sequencing,MPSS)為代表,以雙脫氧終止法(Sangeretal., 1997)和化學(xué)降解法(Whitfeldetal., 1954; Maxametal., 1977)為基礎(chǔ)。第二代測序技術(shù)—高通量測序技術(shù)(High-throughput sequencing)主要包括Roche 454焦磷酸測序技術(shù)、Illumina/Solexa測序技術(shù)、Solexa測序技術(shù)和SOLiD測序技術(shù)等,其原理為提取樣品總RNA后,用高通量測序平臺進(jìn)行樣品文庫的深度測序,經(jīng)圖像識別,獲得DNA序列(張春蘭,2012)。第三代測序技術(shù)包括Helicos單分子測序儀、Pacific Bioscience的SMRT技術(shù)以及Oxford Nanopore Technologies公司的納米孔單分子測序技術(shù)。其中,Illumina測序技術(shù)成本低、周期短、通量高,具有明顯的優(yōu)勢。轉(zhuǎn)錄組測序可用于研究基因功能、可變剪接、新轉(zhuǎn)錄本預(yù)測和結(jié)構(gòu)性變異(Alagnaetal., 2009; Barakatetal., 2009; Dassanayakeetal., 2010; Lietal., 2010)。相對于傳統(tǒng)的芯片雜交平臺,轉(zhuǎn)錄組測序可對任意物種的整體轉(zhuǎn)錄活動(dòng)進(jìn)行檢測,提供更精確的數(shù)字化信號,更高的檢測通量及更廣泛的檢測范圍。因此,對于缺乏基因組信息的物種而言,采用轉(zhuǎn)錄組測序技術(shù)可獲得大量的轉(zhuǎn)錄本信息,從中發(fā)掘重要功能基因的重要研究手段(Franssenetal., 2011; 楊楠等,2012;Lietal., 2013)。
本研究中將Illumina HiSeq 2000高通量測序技術(shù)應(yīng)用到麥紅吸漿蟲轉(zhuǎn)錄組研究,將測序得到的大量數(shù)據(jù)進(jìn)行拼接與組裝,結(jié)合生物信息學(xué)方法對所獲得的Unigene進(jìn)行基因功能注釋、功能分類、代謝途徑分析等,研究結(jié)果和分析為進(jìn)一步的分子標(biāo)記開發(fā)和基因功能研究奠定基礎(chǔ),也為麥紅吸漿蟲的基因組學(xué)研究提供豐富的資源。
麥紅吸漿蟲土樣,與2014年在河南省新鄉(xiāng)市輝縣采集,將所采集的吸漿蟲土樣放入春化室進(jìn)行春化處理,與次年春季移入氣候室,待成蟲羽化當(dāng)天收集,選擇健壯活蟲經(jīng)處理后,立即放入液氮保存。
將液氮冷凍保存的試蟲送至北京諾禾致源生物信息科技有限公司。使用試劑盒提取RNA。采用瓊脂糖凝膠電泳分析RNA降解程度以及是否有污染,Nanodrop檢測RNA的純度(OD260/OD280比值),Qubit對RNA濃度進(jìn)行精確定量,Agilent 2100精確檢測RNA的完整性。將滿足測序要求的RNA樣本進(jìn)行轉(zhuǎn)錄組測序及原始數(shù)據(jù)分析。
1.2.1建庫流程
樣品檢測合格后,用帶有Oligo (dT)的磁珠富集真核生物mRNA(若為原核生物,則通過試劑盒去除rRNA來富集mRNA)。隨后加入fragmentation buffer將mRNA打斷成短片段,以mRNA為模板,用六堿基隨機(jī)引物(random hexamers)合成一鏈cDNA,然后加入緩沖液、dNTPs(dNTP中的dTTP用dUTP取代)和DNA polymerase I和RNase H合成二鏈cDNA,再用AMPure XP beads純化雙鏈cDNA,之后用USER酶降解含有U的cDNA第二鏈。純化的雙鏈cDNA先進(jìn)行末端修復(fù)、加A尾并連接測序接頭,再用AMPure XP beads進(jìn)行片段大小選擇。最后進(jìn)行PCR擴(kuò)增,并用AMPure XP beads純化PCR產(chǎn)物,得到最終的文庫。文庫構(gòu)建完成后,先使用Qubit 2.0進(jìn)行初步定量,稀釋文庫至1.5 ng/μL,隨后使用Agilent 2100對文庫的insert size進(jìn)行檢測,insert size符合預(yù)期后,使用Q-PCR方法對文庫的有效濃度進(jìn)行準(zhǔn)確定量(文庫有效濃度>2 nM),以保證文庫質(zhì)量。庫檢合格后,把不同文庫按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling后進(jìn)行Illumina HiSeq測序。
1.2.2轉(zhuǎn)錄組信息分析
1.2.2.1 測序數(shù)據(jù)的預(yù)處理
測序得到的原始測序序列(Sequenced Reads)或者raw reads,里面含有帶接頭的、低質(zhì)量的reads,為了保證信息分析質(zhì)量,必須對raw reads過濾,得到clean reads,后續(xù)分析都基于clean reads。數(shù)據(jù)處理的步驟如下:(1)去除帶接頭(adapter)的reads;(2)去除N(N表示無法確定堿基信息)的比例大于10%的reads;(3)去除低質(zhì)量reads(質(zhì)量值Qphred <= 20的堿基數(shù)占整個(gè)reads的50%以上的reads)。
1.2.2.2 Trinity拼接
利用Trinity軟件對樣品數(shù)據(jù)進(jìn)行組裝,通過序列之間的overlap信息組裝得到重疊群(Contig),然后在局部進(jìn)行組裝得到轉(zhuǎn)錄本(Transcripts),最后從局部中挑選最主要的轉(zhuǎn)錄本作為單基因簇(Unigene)(Grabherretal., 2011)。
1.2.2.3 轉(zhuǎn)錄本的功能注釋
拼接獲得轉(zhuǎn)錄本后,與7個(gè)公共數(shù)據(jù)庫(Nr,Nt,Pfam,KOG/COG,Swiss-Prot,KEGG和GO)。中的蛋白質(zhì)序列進(jìn)行比對,取閾值e ≤10-10,通過序列相似性進(jìn)行功能注釋。序列相似性比較使用BLAST(Altschuletal.)算法,然后用WEGO軟件對所有的Unigene進(jìn)行GO功能分類統(tǒng)計(jì)(Conesaetal., 2005; Yeetal., 2006)。然后對Unigene分別進(jìn)行蛋白質(zhì)直系同源數(shù)據(jù)庫(euKaryotic Ortholog Groups,KOG)功能分類和東京基因與基金組百科全書(Kyoto encyclopedia of genes and genomes,KEGG)代謝途徑分析(Romanetal., 2003; Minoruetal., 2004)。
1.2.2.4 CDS預(yù)測
CDS預(yù)測分為兩步:1. 將Unigenes按照NR蛋白庫,Swissprot蛋白庫的優(yōu)先級順序進(jìn)行比對,若比對上,則從比對結(jié)果中提取出轉(zhuǎn)錄本的ORF編碼框信息,并按照標(biāo)準(zhǔn)密碼子表將編碼區(qū)序列翻譯成氨基酸序列(按照5′>3′的順序);2. 對于沒有比對上NR和Swissprot蛋白庫的序列,或者比對上未預(yù)測出結(jié)果的序列,則采用estscan (3.0.3)軟件預(yù)測其ORF,從而得到這部分基因編碼的核酸序列和氨基酸序列。
1.2.2.5 SSR分析
利用MISA軟件對麥紅吸漿蟲轉(zhuǎn)錄組中篩選得到1 kb以上的Unigene進(jìn)行簡單重復(fù)序列(Simple sequence repeats,SSR)位點(diǎn)分析,搜索標(biāo)準(zhǔn)為:單、二、三、四、五、六核苷酸基序(motif)至少重復(fù)次數(shù)分別為10、8、5、4、3、3,對查找的SSR類型進(jìn)行特征分析。
轉(zhuǎn)錄組樣本數(shù)據(jù)量為27.88 G,總的CG含量平均值為40.26%,數(shù)據(jù)中96.70%的序列的質(zhì)量參數(shù)Q≥20,93.53%的Q≥30。經(jīng)過分析,共獲得59 257個(gè)Unigenes,總長度49 861 164 bp,最短201 bp,最長29 282 bp,平均長度841 bp(圖1)。
2.2.1Unigene在各數(shù)據(jù)庫中的比對
將獲得的59 257個(gè)Unigenes分別在7個(gè)公共數(shù)據(jù)庫中進(jìn)行比對,取閾值e≤10-10,分別獲得18 596、2 868、9 110、14 680、18 912、19 584和11 279個(gè)結(jié)果,共計(jì)95 029個(gè)結(jié)果。
2.2.2Unigene的GO分類注釋
共有19 584個(gè)Unigenes在GO數(shù)據(jù)庫中3大類50個(gè)功能組中找到對應(yīng)(圖2)。其中生物學(xué)進(jìn)程組中較大的是細(xì)胞進(jìn)程(Cellular process)及新陳代謝進(jìn)程(Metabolic process);細(xì)胞成分組中較大的細(xì)胞(Cell)和細(xì)胞成分(Cell part);分子功能組中較大的是結(jié)合活性(binging)和催化活性(Catalytic activity)。
圖1 拼接Unigenes與轉(zhuǎn)錄本長度分布圖Fig.1 length distribution of splice Unigenes and transcript
圖2 Unigenes的GO分類圖Fig.2 GO functional categories of unigenes
2.2.3Unigene的KOG分類注釋
蛋白質(zhì)直系同源數(shù)據(jù)庫(euKaryotic Ortholog Groups, KOG)是對基因產(chǎn)物進(jìn)行直系同源分類的數(shù)據(jù)庫。將麥紅吸漿蟲Unigenes與KOG數(shù)據(jù)庫進(jìn)行比對,預(yù)測Unigenes功能并進(jìn)行分類統(tǒng)計(jì)。研究結(jié)果表明,共有11 279個(gè)麥紅吸漿蟲Unigenes被注釋,根據(jù)其功能大致可分為26類,對每類的Unigenes進(jìn)行了統(tǒng)計(jì)分析(圖3 A-Z)。從圖中可以看出,Unigenes涉及的KOG功能類別比較全面,涉及了大多數(shù)的生命活動(dòng)。其中,一般功能預(yù)測類基因最多(2 226個(gè));其次是信號傳導(dǎo)機(jī)制類基因(1 467個(gè))、翻譯后修飾,蛋白質(zhì)折疊和分子伴侶(1 265個(gè))、轉(zhuǎn)錄(734個(gè));未知蛋白質(zhì)(1個(gè))最少,其他類別的基因表達(dá)豐度都各不相同。
圖3 麥紅吸漿蟲Unigenes的KOG分類圖Fig.3 KOG functional categories of Sitodiplosis mosellana Unigenes
2.2.3Unigenes的KEGG分類注釋
KEGG(Kyoto Encyclopedia of Genes and Genomes)系統(tǒng)分析基因產(chǎn)物和化合物在細(xì)胞中的代謝途徑以及這些基因產(chǎn)物的功能的數(shù)據(jù)庫。它整合了基因組、化學(xué)分子和生化系統(tǒng)等方面的數(shù)據(jù),包括代謝通路(KEGG PATHWAY)、藥物(KEGG DRUG)、疾病(KEGG DISEASE)、功能模型(KEGG MODULE)、基因序列(KEGG GENES)及基因組(KEGG GENOME)等等。結(jié)合KEGG數(shù)據(jù)庫,對麥紅吸漿蟲的Unigene可能參與或涉及的代謝途徑進(jìn)行了統(tǒng)計(jì)分析。結(jié)果顯示,共有9 110個(gè)麥紅吸漿蟲Unigenes被注釋,分別歸屬于細(xì)胞進(jìn)程、環(huán)境信息進(jìn)程、遺傳信息進(jìn)程、新陳代謝和有機(jī)體系統(tǒng)五大類代謝途徑,主要包括細(xì)胞生長與死亡、細(xì)胞運(yùn)動(dòng)、信號轉(zhuǎn)導(dǎo)、能量代謝等32類代謝途徑。其中占比例最高的是信號轉(zhuǎn)導(dǎo)有1 078個(gè)(占11.83%),其次是轉(zhuǎn)錄有730個(gè)(占8.01%)和折疊、排序與退化有661個(gè)(占7.26%),比例最小的是生物合成與次生代謝僅有52個(gè)(占0.57%)(圖4)。
CDS預(yù)測結(jié)果顯示,共30 088條序列可被編碼,占全部基因的50.78%(30088/59257)。圖5表示所預(yù)測CDS的長度統(tǒng)計(jì)。
SSR,簡單重復(fù)序列(simple sequence repeats),又稱為短串聯(lián)重復(fù)序列和微衛(wèi)星標(biāo)記。是一類由幾個(gè)核苷酸(1~5個(gè))為重復(fù)單位組成的長達(dá)幾十個(gè)核苷酸的重復(fù)序列,長度較短,且廣泛均勻分布于真核生物基因組中(Chambersetal., 2000)。由于重復(fù)單位的次數(shù)不同或重復(fù)程度的不完全相同,造成了SSR長度的高度變異性,其中最常見的雙核苷酸重復(fù)類型,如(CA) n(Jérmeetal., 2010)。本研究分析中采用MISA 1.0,默認(rèn)參數(shù);對應(yīng)的各個(gè)unit size的最少重復(fù)次數(shù)分別為(1~10,2~6,3~5,4~5,5~5,6~5)對Unigene進(jìn)行SSR檢測,共檢測到36 323個(gè)微衛(wèi)星序列,發(fā)生率為61.30%(36 323/59 257)。其中單核苷酸重復(fù)、二核苷酸重復(fù)、三核苷酸重復(fù)、四核苷酸重復(fù)、五核苷酸重復(fù)、六核苷酸重復(fù)微衛(wèi)星序列分別為21 508、9 217、5 290、261、22、25條,依次占微衛(wèi)星序列的59.21%(21 508/36 323)、25.38%、14.56%、0.72%、0.06%、0.07%,其中單核苷酸重復(fù)占比最高,其次是二核苷酸和三核苷酸重復(fù)。并且對不同SSR類型在基因轉(zhuǎn)錄本的密度分布進(jìn)行統(tǒng)計(jì)(圖6)。
圖4 麥紅吸漿蟲Unigenes的KEGG分類圖Fig.4 KEGG functional categories of Sitodiplosis mosellana Unigenes
圖5 麥紅吸漿蟲Unigenes的CDS的長度分布統(tǒng)計(jì)圖Fig.5 Length distribution of CDS of Sitodiplosis mosellana unigenes
圖6 麥紅吸漿蟲Unigene的不同SSR類型在基因轉(zhuǎn)錄本的密度分布Fig.6 Density distribution of Unigene SSR in different types gene transcripts of Sitodiplosis mosellana
轉(zhuǎn)錄組學(xué)(Transcriptomics)是功能基因組學(xué)研究的一個(gè)重要內(nèi)容,它是從整體水平上研究細(xì)胞中基因轉(zhuǎn)錄的情況及其轉(zhuǎn)錄調(diào)控規(guī)律?;诟咄繙y序技術(shù)的轉(zhuǎn)錄組測序(RNA-seq)通過對組織中的RNA(包括mRNA和非編碼RNA)進(jìn)行測序,能夠全面快速地獲得某一物種特殊組織或器官在某一特定狀態(tài)下的幾乎所有轉(zhuǎn)錄本信息,具有高準(zhǔn)確性、高通量、高靈敏度和低運(yùn)行成本等突出優(yōu)勢,已經(jīng)廣泛應(yīng)用于各種生物轉(zhuǎn)錄組的研究(Shendure, 2008; Grabherretal., 2011; 張春蘭等,2012)。為獲得麥紅吸漿蟲的轉(zhuǎn)錄組信息,本研究應(yīng)用Illumina高通量測序技術(shù)對麥紅吸漿蟲進(jìn)行了轉(zhuǎn)錄組測序,通過序列的預(yù)處理Trinity拼接、 轉(zhuǎn)錄本功能注釋等步驟,得到了數(shù)據(jù)量為27.88 G的轉(zhuǎn)錄組信息。大量的麥紅吸漿蟲轉(zhuǎn)錄本信息可為其今后的分子標(biāo)記和功能基因挖掘提供有價(jià)值的信息。
通過高通量測序,經(jīng)拼接組裝,共獲得77 500個(gè)轉(zhuǎn)錄本,進(jìn)一步分析獲得59 257個(gè)Unigenes,總長度49 861 164 bp,最短201 bp,最長29 282 bp,平均長度841 bp,N50值為1 802 bp(N50是指將拼接轉(zhuǎn)錄本按照長度從長到短排序,累加轉(zhuǎn)錄本的長度,到不小于總長50%的拼接轉(zhuǎn)錄本的長度),N50值越大,反映組裝得到的長片段越多,組裝效果就越好。測序數(shù)據(jù)產(chǎn)量和數(shù)據(jù)組裝質(zhì)量是轉(zhuǎn)錄組測序完成情況的重要指標(biāo)。以上研究結(jié)果表明,此次序列組裝的質(zhì)量和長度可以滿足轉(zhuǎn)錄組分析的基本要求,且新一代高通量測序技術(shù)是批量麥紅吸漿蟲功能基因的更為有效手段,進(jìn)一步說明Illumina HiSeq 2000是高通量轉(zhuǎn)錄組測序的可靠平臺。
結(jié)合生物信息學(xué)分析方法,對麥紅吸漿蟲Unigene與7個(gè)公共數(shù)據(jù)庫(Nr,Nt,Pfam,KOG/KOG,Swiss-Prot,KEGG和GO)進(jìn)行比對,以強(qiáng)大的生物信息學(xué)平臺作支撐,根據(jù)“基因結(jié)構(gòu)相似,功能同源”的原理,對基因的功能進(jìn)行注釋。分別獲得18 596、2 868、9 110、14 680、18 912、19 584和11 279個(gè)結(jié)果,共計(jì)95 029個(gè)結(jié)果。表明在對麥紅吸漿蟲基因組不清楚的情況下,高通量測序技術(shù)是批量發(fā)現(xiàn)麥紅吸漿蟲功能基因的有效且簡易手段。
本研究中利用Nr和SwissProt蛋白數(shù)據(jù)庫對所獲得的59 257個(gè)Unigenes進(jìn)行BLASTX比對分析,共有30 088條序列可被編碼,占全部基因的50.78%(30 088/59 257)。29 169個(gè)Unigenes與其它物種蛋白序列無匹配,占總體的49.22%,此部分Unigenes主要包括以下3種類型:(a)Unigenes序列片段長度過短,不能獲得同源性比對結(jié)果;(b)基因注釋信息的暫時(shí)缺乏。昆蟲基因組學(xué)及轉(zhuǎn)錄組學(xué)研究剛剛起步,生物信息數(shù)據(jù)庫仍在不斷更新和完善中,基因功能注釋信息不全會造成部分序列暫時(shí)無法獲得對應(yīng)的功能注釋信息;(c)昆蟲特有的新基因,可能存在一些昆蟲特有的新基因,這或許也是導(dǎo)致其同源序列較難發(fā)現(xiàn)的原因之一。隨著今后研究的深入,進(jìn)一步將麥紅吸漿蟲與其它昆蟲進(jìn)行比較分析,為研究昆蟲的分子生物學(xué)提供更為重要的信息。SSR分子標(biāo)記具有操作簡便、重復(fù)性好、多態(tài)性豐富、遺傳信息量大和共顯性遺傳等優(yōu)點(diǎn),已在遺傳多樣性分析、遺傳圖譜構(gòu)建、功能基因發(fā)掘、分子標(biāo)記輔助育種等研究中得到了廣泛應(yīng)用(宋建等,2006;Yangetal., 2009; 劉峰等,2012)。采取實(shí)驗(yàn)室手段開發(fā)SSR引物費(fèi)時(shí)、費(fèi)力、成本高且試驗(yàn)復(fù)雜,基于轉(zhuǎn)錄組數(shù)據(jù)庫信息進(jìn)行SSR分子標(biāo)記開發(fā)將是一種既經(jīng)濟(jì)又有效的方法。本研究通過查找發(fā)現(xiàn)了36 323個(gè)SSR位點(diǎn),SSR不但出現(xiàn)頻率高(發(fā)生率為61.30%),而且類型豐富。目前,已經(jīng)有很多研究者通過轉(zhuǎn)錄組數(shù)據(jù)庫發(fā)掘分子標(biāo)記(Zhaietal., 2013; Lietal., 2016),解決了傳統(tǒng)開發(fā)分子標(biāo)記的方法耗時(shí)長、成本高的問題。本研究將為麥紅吸漿蟲后續(xù)分子鑒定、多態(tài)性檢測和種群遺傳結(jié)構(gòu)的研究提供一定的參考價(jià)值。
本研究采用Illumina HiSeq 2000高通量測序技術(shù)建立了麥紅吸漿蟲轉(zhuǎn)錄組數(shù)據(jù)庫,獲得了大量的轉(zhuǎn)錄本信息,并對表達(dá)基因進(jìn)行了序列組裝、功能注釋及代謝途徑等分析,為今后更深入研究麥紅吸漿蟲功能基因組、基因克隆及麥紅吸漿蟲抗蟲品種的研究提供了大量的信息資源,該轉(zhuǎn)錄組數(shù)據(jù)也還可以作為今后基因組的參考序列,為麥紅吸漿蟲的分子生物學(xué)研究提供寶貴的基因組數(shù)據(jù)來源。