李媛媛，趙金紅

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

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

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

1 材料和方法

1.1 試蟲

供試的玉米象成蟲采集于安徽省蕪湖市糧庫，按照文獻(鄧樹華等，2011)培養(yǎng)條件進行培養(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進行電泳檢測，Nanodrp檢測RNA的純度(OD260/OD280)，利用Qubit對RNA濃度進行精確定量，RNA的完整性用Agilent 2100進行檢測。將滿足測序要求的RNA樣本進行轉錄組測序。

1.2.2文庫構建和測序

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

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

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

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

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

2 結果與分析

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

對獲得的原始數(shù)據(jù)進行統(tǒng)計，共產(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進行注釋(表1)。不同數(shù)據(jù)庫中注釋成功的unigenes基因數(shù)及其所占比例有所差別。有915條unigenes同時在7個數(shù)據(jù)庫中注釋成功，占總數(shù)的1.42%；在這7個數(shù)據(jù)庫中至少有一個數(shù)據(jù)庫注釋成功的有26 580條，占總數(shù)的41.3%。在Nr數(shù)據(jù)庫當中，有25 271條unigenes能夠?qū)ひ挼浇菩蛄?，約占unigenes總數(shù)的39.27%；在GO數(shù)據(jù)庫當中，有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)、細胞組分(cellular component)和分子功能(molecular function)3大類(圖2)。

所有大類可細分為68個二級分類。對unigenes在二級分類中的分布情況進行統(tǒng)計分析，“生物學過程”包含26個不同的亞類，它也是3大類中最多的亞類，其中細胞過程(cellular process)和生物調(diào)節(jié)(biological regulation)所占比例較高，分別為13 426條和7 814條，分別占GO注釋信息總數(shù)的67.71%和37.66%；在“細胞組分”類別中有22個亞類，其中細胞(cell)和細胞部分(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ù)庫中比對玉米象轉錄組unigenes，結果顯示共12 060條unigenes獲得注釋信息，根據(jù)其功能大致可分為25類(圖3)。其中，信號傳導機制(Signal transduction mechanisms)、一般功能預測基因(General function prediction only)、翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(Posttranslational modification, protein turnover， chaperones)獲得注釋最多，分別為1 962條、1 648條、1 135條，而細胞活性(Cell motility)獲得注釋最少，僅20條。

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

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

2.2.3Unigenes的KEGG分類注釋

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

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

2.3 玉米象SNP特征分析

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

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

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

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

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

3 結論與討論

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

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

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

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

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

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

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