基于高通量測序的玉米象轉(zhuǎn)錄組分析

李媛媛，趙金紅

(皖南醫(yī)學(xué)院醫(yī)學(xué)寄生蟲學(xué)教研室，安徽蕪湖 241002)

玉米象Sitophiluszeamais屬鞘翅目Coleoptera象甲科Curculionidae，別名米牛、鐵嘴，在全世界范圍內(nèi)均有分布，是儲糧的世界性頭號害蟲。玉米象可孳生于多種谷物及加工品、干果、豆類、藥材、油料等，尤其是陰暗潮濕的倉庫為最為嚴(yán)重(鄭旭等，2014；姜洪等，2018)，不僅可引起儲糧中谷物數(shù)量減少，而且還可改變谷物及其產(chǎn)品的質(zhì)量，導(dǎo)致種子活力的退化(Haqetal.，2000；王鵬杰等，2020)。目前對于玉米象的研究主要集中在生物、生態(tài)學(xué)特性以及防治研究方面，對其使用化學(xué)合成熏蒸劑進(jìn)行防治，但長期使用化學(xué)制劑可導(dǎo)致玉米象產(chǎn)生嚴(yán)重抗藥性(呂建華等，2010；侯昌亮等，2016)。目前，由于缺乏在基因組和轉(zhuǎn)錄組方面研究，有效的分子標(biāo)記尚未被發(fā)現(xiàn)，導(dǎo)致玉米象的抗藥性和種群遺傳結(jié)構(gòu)的研究尚未開展。因此，本文通過獲得大量轉(zhuǎn)錄組信息，以期為功能基因的挖掘提供信息資源，為玉米象種群遺傳結(jié)構(gòu)、抗藥性監(jiān)測以及防治奠定理論基礎(chǔ)。

轉(zhuǎn)錄組廣義上是指生物體細(xì)胞或組織在特定狀態(tài)下所轉(zhuǎn)錄出來的所有RNA的總和(張春蘭等，2012)。轉(zhuǎn)錄組學(xué)是從整體轉(zhuǎn)錄水平系統(tǒng)研究基因轉(zhuǎn)錄圖譜，并揭示復(fù)雜生物學(xué)通路和性狀調(diào)控網(wǎng)絡(luò)分子機(jī)制的學(xué)科(崔凱等，2019)。隨著測序技術(shù)的發(fā)展，核酸的檢測和定量更加便捷和準(zhǔn)確，越來越多的研究學(xué)者將高通量測序技術(shù)應(yīng)用到轉(zhuǎn)錄組研究中。對于危害農(nóng)作物的昆蟲轉(zhuǎn)錄組研究工作也有相關(guān)報(bào)道，如大墊尖翅蝗Epacromiuscoerulipes(金永玲等，2015)、東方粘蟲Mythimnaseparata(胡艷華等，2015；李微等，2017)、溝眶象Eucryptorrhynchuschinensis(武政梅，2016)、疣蝗Trilophidiaannulata(邱忠營等，2016)、印度谷螟Plodiainterpunctella(李慧等，2019)、麥紅吸漿蟲Sitodiplosismosellana(蔣月麗等，2020)、紅棕象甲Rhynchophorusferrugineus(Yangetal.，2020)。本研究采用高通量測序技術(shù)對玉米象進(jìn)行轉(zhuǎn)錄組測序、拼接與組裝，對所獲得的unigenes進(jìn)行基因功能注釋、代謝途徑及SSR檢測等分析，為玉米象進(jìn)一步的分子標(biāo)記開發(fā)和基因功能研究提供分子基礎(chǔ)。

1 材料和方法

1.1 試蟲

供試的玉米象成蟲采集于安徽省蕪湖市糧庫，按照文獻(xiàn)(鄧樹華等，2011)培養(yǎng)條件進(jìn)行培養(yǎng)。恒溫培養(yǎng)箱：溫度27℃±2℃，相對濕度75%±5%，光周期L ∶D=16 h ∶8 h，小麥飼養(yǎng)7 d后，篩去成蟲，待新一代成蟲羽化大量出現(xiàn)一周左右，選取2頭羽化后7～14 d的健壯活潑玉米象成蟲作為研究對象，用無菌水清洗蟲體，液氮速凍后保存于-80℃冰箱備用。

1.2 方法

1.2.1總RNA提取

采用Total RNA Extractor試劑盒提取液氮冷凍保存的玉米象總RNA，將提取的RNA進(jìn)行電泳檢測，Nanodrp檢測RNA的純度(OD260/OD280)，利用Qubit對RNA濃度進(jìn)行精確定量，RNA的完整性用Agilent 2100進(jìn)行檢測。將滿足測序要求的RNA樣本進(jìn)行轉(zhuǎn)錄組測序。

1.2.2文庫構(gòu)建和測序

樣品檢測合格后，先采用Oligo(dT)方法分離并純化出mRNA，隨后將mRNA打斷為短片段，用隨機(jī)引物合成一鏈cDNA(互補(bǔ)DNA)和二鏈cDNA，對雙鏈cDNA進(jìn)行純化、末端修飾、片段大小的選擇，最后進(jìn)行PCR富集得到最終的cDNA文庫。使用Agilent 2100對文庫進(jìn)行準(zhǔn)確定量，以保證文庫的質(zhì)量。庫檢合格后進(jìn)行Illumina Hiseq測序。

1.2.3數(shù)據(jù)組裝與功能注釋

采用Illumina Hiseq測序得到原始圖像數(shù)據(jù)，經(jīng)堿基識別轉(zhuǎn)化為序列數(shù)據(jù)，分別去除帶接頭的序列、無法確定的堿基數(shù)量(N%)>10%的序列和低質(zhì)量(質(zhì)量值≤Q20的堿基數(shù)占50%以上)的序列。利用Trinity軟件對樣品數(shù)據(jù)進(jìn)行組裝拼接后，分別與7個數(shù)據(jù)庫(Nr、GO、Swissprot、KOG、Pfam、KEGG、TrEMBL)中的蛋白質(zhì)序列進(jìn)行BLAST比對，從而對其功能進(jìn)行注釋。

1.2.4簡單重復(fù)序列(SSR)及單核苷酸多態(tài)性(SNP)分析

利用MISA軟件對玉米象中篩選得到的1 kb以上的Unigene進(jìn)行SSR位點(diǎn)分析。通過Bcftools軟件對單核苷酸多態(tài)性(SNP)進(jìn)行鑒定。

2 結(jié)果與分析

2.1 玉米象轉(zhuǎn)錄組測序質(zhì)控與組裝

對獲得的原始數(shù)據(jù)進(jìn)行統(tǒng)計(jì)，共產(chǎn)生44 277 838條原始reads，總堿基數(shù)為6.64 G。再消除帶接頭和低質(zhì)量的reads后，獲得pair-end Reads總數(shù)為43 331 830，總堿基數(shù)為6.34 G，GC含量為38.66%，Q20堿基比例平均超過99.27%。經(jīng)過組裝，共獲得64 358條unigenes，總長度39 481 752 bp，最短201 bp，最長29 046 bp，平均長度613 bp，N50長度為871 bp。(圖1)

圖1 組裝序列長度分布圖Fig.1 Size distribution of transcripts and unigenes

2.2 功能注釋

分別在7個數(shù)據(jù)庫中對組裝得到的64 358條unigenes進(jìn)行注釋(表1)。不同數(shù)據(jù)庫中注釋成功的unigenes基因數(shù)及其所占比例有所差別。有915條unigenes同時(shí)在7個數(shù)據(jù)庫中注釋成功，占總數(shù)的1.42%；在這7個數(shù)據(jù)庫中至少有一個數(shù)據(jù)庫注釋成功的有26 580條，占總數(shù)的41.3%。在Nr數(shù)據(jù)庫當(dāng)中，有25 271條unigenes能夠?qū)ひ挼浇菩蛄校s占unigenes總數(shù)的39.27%；在GO數(shù)據(jù)庫當(dāng)中，有20 747條unigenes獲取了注釋，約為總數(shù)的32.24%；在Pfam和KEGG數(shù)據(jù)庫匯中獲得注釋的unigenes數(shù)量都10 000條以下，分別為9 471(占總體數(shù)的14.78%)和3 370(占總體數(shù)的5.24%)。

2.2.1Unigenes的GO分類注釋

按照GO功能分類方式，將GO注釋的20 747條unigenes分為生物過程(biological process)、細(xì)胞組分(cellular component)和分子功能(molecular function)3大類(圖2)。

所有大類可細(xì)分為68個二級分類。對unigenes在二級分類中的分布情況進(jìn)行統(tǒng)計(jì)分析，“生物學(xué)過程”包含26個不同的亞類，它也是3大類中最多的亞類，其中細(xì)胞過程(cellular process)和生物調(diào)節(jié)(biological regulation)所占比例較高，分別為13 426條和7 814條，分別占GO注釋信息總數(shù)的67.71%和37.66%；在“細(xì)胞組分”類別中有22個亞類，其中細(xì)胞(cell)和細(xì)胞部分(cell part)類型的unigenes較高，分別為13 710(66.08%)和13 664(65.86%)；在“分子功能”類別中有20個亞類，其中連接(binding)和催化活性(catalytic activity)的unigenes較高，分別為12 704(61.23%)和9 356(45.1%)。

2.2.2Unigenes的KOG分類注釋

在KOG數(shù)據(jù)庫中比對玉米象轉(zhuǎn)錄組unigenes，結(jié)果顯示共12 060條unigenes獲得注釋信息，根據(jù)其功能大致可分為25類(圖3)。其中，信號傳導(dǎo)機(jī)制(Signal transduction mechanisms)、一般功能預(yù)測基因(General function prediction only)、翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(Posttranslational modification, protein turnover， chaperones)獲得注釋最多，分別為1 962條、1 648條、1 135條，而細(xì)胞活性(Cell motility)獲得注釋最少，僅20條。

圖2 玉米象GO功能分類圖Fig.2 Statistical diagram of GO function of Sitophilus zeamais注：1,抗氧化功能；2,連接功能；3,催化功能；4,趨化活性；5,化學(xué)抗拮；6,電子轉(zhuǎn)移活性；7,金屬伴侶活性；8,分子功能調(diào)節(jié)；9,分子傳導(dǎo)活性；10,形態(tài)發(fā)生活性；11,營養(yǎng)庫活性；12,蛋白標(biāo)志物；13,信號轉(zhuǎn)導(dǎo)活性；14,結(jié)構(gòu)分子活性；15,轉(zhuǎn)錄因子活性，蛋白質(zhì)結(jié)合；16,翻譯調(diào)節(jié)功能；17,載體活性；18,通道調(diào)節(jié)活性；19,酶調(diào)節(jié)功能；20,受體調(diào)節(jié)功能；21,生物習(xí)性；22,生物黏附；23,生物相；24,生物調(diào)節(jié)；25,細(xì)胞聚集；26,細(xì)胞殺傷；27,組織或生物源的細(xì)胞組分；28,細(xì)胞過程；29,解毒；30,發(fā)育過程；31,生長；32,免疫反應(yīng)過程；33,定位；34,移動；35,新陳代謝過程；36,多有機(jī)體過程；37,多細(xì)胞生物過程；38,生物負(fù)調(diào)控過程；39,生物正調(diào)控過程；40,生物過程調(diào)控；41,生殖；42,生殖過程；43,刺激反應(yīng)；44,節(jié)律進(jìn)程；45,信號；46,建立定位；47,細(xì)胞；48,細(xì)胞結(jié)合；49,細(xì)胞成分；50,細(xì)胞外基質(zhì)；51,細(xì)胞外基質(zhì)成分；52,細(xì)胞外區(qū)域；53,細(xì)胞外區(qū)域成分；54,蛋白質(zhì)復(fù)合物；55,細(xì)胞膜；56,細(xì)胞膜成分；57,膜包圍的內(nèi)腔；58,內(nèi)核；59,細(xì)胞器；60,細(xì)胞器成分；61,其他有機(jī)體；62,其他有機(jī)體成分；63,超分子纖維；64,共質(zhì)體；65,突觸；66,突觸成分；67,病毒；68,病毒部分。

圖3 玉米象KOG功能統(tǒng)計(jì)圖Fig.3 Statistical diagram of KOG function of Sitophilus zeamais

2.2.3Unigenes的KEGG分類注釋

使用KEGG數(shù)據(jù)庫對玉米象的unigenes可能參與或涉及的代謝途徑進(jìn)行了統(tǒng)計(jì)分析，在組裝得到的unigenes中，有3 370條得到注釋。這些通路信息分可為5大類，分別為細(xì)胞進(jìn)程(Cellular Processes)、環(huán)境信息處理(Environmental Information Processing)、遺傳信息處理(Genetic Information Processing)、新陳代謝(Metabolism)和有機(jī)系統(tǒng)(Organismal Systems)，這5大類可以進(jìn)一步分為33個亞類274個功能通路(圖4)。其中信號轉(zhuǎn)導(dǎo)(Signal transduction)獲得注釋信息最多，有724條，參與較多的幾條信號轉(zhuǎn)導(dǎo)通路分別是PI3K-Akt信號通路(102條)、MAPK信號通路(89條)、mTOR信號通路(89條)、cMAP信號通路(79條)，這些代謝通路與環(huán)境信息大類中信號轉(zhuǎn)導(dǎo)有關(guān)；翻譯(Translation)獲得的注釋信息有435條，位居第2位；膜運(yùn)輸(Membrane transport)和其它此生代謝產(chǎn)物的生物合成(Biosynthesis of other secondary metabolites)的注釋最少，僅有24條和26條。

圖4 玉米象KEGG通路分類統(tǒng)計(jì)圖Fig.4 Statistical diagram of KEGG Pathway of Sitophilus zeamais

2.3 玉米象SNP特征分析

對玉米象序列進(jìn)行檢測，共發(fā)現(xiàn)146 081個單核苷酸多態(tài)性(SNP，Single Nucleotide Polymorphsims)位點(diǎn)。對這些SNP位點(diǎn)進(jìn)行類型統(tǒng)計(jì)，結(jié)果顯示轉(zhuǎn)換突變類型的SNP有102 850個(占70.41%)，顛換突變類型的SNP有43 231個(占29.59%)。在轉(zhuǎn)換突變類型中，由胞嘧啶轉(zhuǎn)換為胸腺嘧啶的突變最多(25 808個)，其次為鳥嘌呤轉(zhuǎn)換為腺嘌呤的突變(25 786個)；在顛換突變類型中，由腺嘌呤顛換為胸腺嘧啶的突變最多(7 152個)，而鳥嘌呤顛換為胞嘧啶的突變最少(4 084個)(圖5)。

圖5 玉米象SNP突變類型統(tǒng)計(jì)圖Fig.5 Statistical diagram of SNP mutation of Sitophilus zeamais

2.4 玉米象微衛(wèi)星位點(diǎn)(SSR)分析

使用軟件MISA對玉米象轉(zhuǎn)錄組中的64 358條unigenes進(jìn)行搜索，共發(fā)現(xiàn)5 002個SSR位點(diǎn)，出現(xiàn)頻率為7.772%(檢測出的SSR數(shù)量與總unigenes數(shù)目的比值)，分布于4 606條unigenes中，發(fā)生頻率為7.157%(含有SSR的unigenes的數(shù)目與總unigenes數(shù)目的比值)。玉米象SSR重復(fù)類型豐富，一至六核苷酸重復(fù)均有出現(xiàn)，但是數(shù)量分布差異較大。單核苷酸重復(fù)所占比例最大，為84.09%；其次為三核苷酸和二核苷酸重復(fù)，分別所占比例為9.42%和5.12%；四核苷酸重復(fù)、五核苷酸重復(fù)和六核苷酸重復(fù)分別占1.10%、0.26%和0.01%(表2)。其中，發(fā)生從分布情況來看，玉米象轉(zhuǎn)錄組序列中平均每11.805 kb(序列總長度與SSR總數(shù)目的比值)出現(xiàn)一個SSR，表明該物種轉(zhuǎn)錄組SSR數(shù)量較為豐富。

表2 玉米象轉(zhuǎn)錄組中SSR不同重復(fù)類型數(shù)量和分布

3 結(jié)論與討論

轉(zhuǎn)錄組學(xué)能從整體水平上研究細(xì)胞中基因轉(zhuǎn)錄的情況及其轉(zhuǎn)錄調(diào)控規(guī)律(Shendure，2008；Grabherr，2011)，通過高通量測序技術(shù)和生物信息分析手段對玉米象的RNA進(jìn)行測序和分析，獲得的轉(zhuǎn)錄組數(shù)據(jù)可為玉米象分子生物學(xué)研究提供基因組數(shù)據(jù)來源。

本研究在無參考基因組的情況下，采用Illumina HiSeqTM測序技術(shù)對玉米象進(jìn)行轉(zhuǎn)錄組測序，分析玉米象的轉(zhuǎn)錄組特征，共獲得6.64 G原始數(shù)據(jù)，組裝后得到64 358條unigenes，平均長度分別613 bp，N50長度為871 bp。結(jié)合生物信息學(xué)分析方法，對玉米象unigenes與7個公共數(shù)據(jù)庫(Nr、GO、Swissprot、KOG、Pfam、KEGG、TrEMBL)進(jìn)行比對分析，分別有25 271、20 747、16 477、12 060、9 471、3 370、25 500條unigenes獲得注釋，在7個數(shù)據(jù)庫中同時(shí)得到注釋的unigenes有26 580條，占全部unigenes總數(shù)的41.3%，但仍有37 778條(58.7%)未得到準(zhǔn)確的定位，出現(xiàn)這一現(xiàn)象的原因，可能與目前測序技術(shù)的局限性組裝得到的unigenes片段長度太短、數(shù)據(jù)庫基因注釋信息不足以及該物種存在新基因等因素有關(guān)。

其中，與Nr數(shù)據(jù)庫進(jìn)行比對時(shí)，注釋的25 271條unigenes與中歐山松大小蠹Dendroctonusponderosae匹配率最高，為51.40%，這個結(jié)果和之前學(xué)者對溝眶象E.chinensis(武政梅，2016)的轉(zhuǎn)錄組進(jìn)行研究時(shí)，得出的結(jié)論相同，其對獲得的65 186條unigenes進(jìn)行功能注釋，有21 742條注釋到Nr數(shù)據(jù)庫，與中歐山松大小蠹的匹配率最高達(dá)到49.41%。

在GO注釋的unigenes中，與細(xì)胞過程和生物調(diào)節(jié)、細(xì)胞組件及連接催化活性相關(guān)的基因最多，這些可能與玉米象細(xì)胞的增殖和分化，自身的生長發(fā)育及生命活動過程有著密切聯(lián)系。通過比對KOG數(shù)據(jù)庫，獲得12 060個unigenes注釋信息，對玉米象的unigenes的功能分布狀況有了初步的了解。

KEGG通路注釋分類中，共有3 370條玉米象unigenes參與到了274個代謝通路中，其中信號轉(zhuǎn)導(dǎo)獲得注釋最多，參與PI3K-Ak(102條)、MAPK(89條)、mTOR(89條)信號通路的unigenes數(shù)量最多，這些信號通路在細(xì)胞生長、分化和凋亡方面都起到重要的作用，目前上述信號通路作為腫瘤領(lǐng)域的研究熱點(diǎn)，被證明在多種人類腫瘤中異常激活，通路中的某些組分的突變可引起細(xì)胞轉(zhuǎn)化、調(diào)控腫瘤細(xì)胞的存活、增殖及遷移，同時(shí)抑制自噬(舒婷等，2016；柳望舒等，2016)。在現(xiàn)今大多數(shù)昆蟲抗藥性研究中，功能基因組技術(shù)已被廣泛用于候選基因調(diào)控抗藥性機(jī)制的研究(陳皓等，2020)；隨著昆蟲基因組學(xué)、蛋白質(zhì)組學(xué)、遺傳學(xué)與分子生物學(xué)的應(yīng)用，昆蟲抗藥性機(jī)制的研究已取得突破性進(jìn)展(吳有剛等，2019)。通過對玉米象轉(zhuǎn)錄組的注釋分析，可以從代謝調(diào)控中找到相關(guān)功能基因，并了解基因在代謝通路中的調(diào)控位置，對后續(xù)玉米象的抗藥性、防治研究方面提供了基礎(chǔ)數(shù)據(jù)。

SSR分子標(biāo)記具有操作簡便、多態(tài)性好、共顯性以及數(shù)量豐富等特點(diǎn)(何海等，2015)。目前高通量測序可以從上萬條基因序列中獲得大量的微衛(wèi)星。本研究對拼接得到的64 358條unigenes進(jìn)行分析，發(fā)現(xiàn)5 002個SSR位點(diǎn)，統(tǒng)計(jì)發(fā)現(xiàn)單核苷酸和三核苷酸重復(fù)所占比例最大，分別為84.09%和9.42%。在之前的研究中，如溫帶臭蟲Cimeslectularius(李敏等，2019)的SSR、東方粘蟲(李微等，2017)的SSR、星天牛Anoplophorachinensis(韓小紅等，2019)的SSR、大猿葉甲Colaphellusbowringi(沙君雪等，2018)以單核苷酸為主要類型，其次為三核苷酸重復(fù)，與本研究相符；疣蝗(邱忠營等，2016)的SSR主要以二核苷酸為主要重復(fù)類型；與玉米象同屬于象甲科的溝眶象(武政梅等，2016)的SSR以三核苷酸為主。不同種屬的昆蟲SSR的主要重復(fù)類別不同，可以通過對SSR的密度、分布特點(diǎn)對不同的種進(jìn)行區(qū)分。

目前對于象甲科昆蟲分子生物學(xué)的研究大多集中在線粒體全基因組序列的主要結(jié)構(gòu)特征特點(diǎn)分析上，轉(zhuǎn)錄組的相關(guān)報(bào)道并不多見。近期，有學(xué)者對和玉米象同屬于隱頦象亞科Dryophthorinae的紅棕象甲的轉(zhuǎn)錄組進(jìn)行了研究，獲得了蛹、幼蟲和成蟲3個時(shí)期大量的轉(zhuǎn)錄組數(shù)據(jù)，對昆蟲發(fā)育相關(guān)基因轉(zhuǎn)錄組數(shù)據(jù)的分析將有助于害蟲防治。本研究獲得的數(shù)據(jù)可為玉米象后續(xù)的多態(tài)性監(jiān)測、種群遺傳結(jié)構(gòu)提供基礎(chǔ)資料，對玉米象的生物防治和其抗藥性的研究提供參考依據(jù)。