銹赤扁谷盜的全長轉(zhuǎn)錄組測序分析

胡懷月 陳二虎 韋磊 王康旭 李孟凡 唐培安

摘要：銹赤扁谷盜是一種世界性儲糧害蟲，但對其分子生物學(xué)方面的基礎(chǔ)研究較少。研究基于PacBio IsoSeq平臺對銹赤扁谷盜進行全長轉(zhuǎn)錄組測序，并對其數(shù)據(jù)進行生物學(xué)分析。測序數(shù)據(jù)經(jīng)過校正后，銹赤扁谷盜全長轉(zhuǎn)錄組平均長度為1 845.06 bp，N50長度為2 130 bp；經(jīng)去冗余后，獲得轉(zhuǎn)錄本32 855條。在Nr、KEGG、KOG、Swiss-Prot四大數(shù)據(jù)庫中，分別有29 404、27 794、20 286、21 197條轉(zhuǎn)錄本被注釋；其中，19 346條轉(zhuǎn)錄本在4個數(shù)據(jù)庫中均有注釋，29 416條轉(zhuǎn)錄本至少在一個數(shù)據(jù)庫中有注釋，占總數(shù)的89.53%。此外，經(jīng)鑒定，獲得1 139個轉(zhuǎn)錄因子（TFs）、2 636條長鏈非編碼RNA（LncRNA）和4 197個SSR，預(yù)測后獲得32 800條CDS。

關(guān)鍵詞：銹赤扁谷盜；全長轉(zhuǎn)錄組；基因組注釋

中圖分類號：S379.5 文獻標志碼：A DOI：10.16465/j.gste.cn431252ts.20230326

基金項目：國家重點研發(fā)計劃專項（2021YFD2100604-01）；江蘇省重點研發(fā)計劃項目（BE2022377）；國家自然科學(xué)基金項目（32272388）；江蘇省高校優(yōu)勢學(xué)科建設(shè)工程資助（PAPD）；江蘇省研究生科研與實踐創(chuàng)新計劃（KYCX21_1528）。

Sequencing analysis of the full-length transcriptome data of Cryptolestes ferrugineus

Hu Huaiyue1， Chen Erhu1， Wei Lei2， Wang Kangxu1， Li Mengfan1， Tang Peian1

（ 1. College of Food Science and Engineering， Nanjing University of Finance and Economics/ Collaborative Innovation Center for Modern Grain Circulation and Safety， Nanjing， Jiangsu 210023； 2. Jiangsu Wanjiafu Rice Industry Co.， Taizhou， Jiangsu 225700 ）

Abstract： Cryptolestes ferrugineus is a worldwide grain storage pest， but there are few basic studies on its molecular biology. In this study， we sequenced the full-length transcriptome of C. ferrugineus based on PacBio Iso-Seq platform and analyzed the data biologically. After the sequencing data were corrected， the average length of the full-length transcriptome was 1 845.06 bp， and the N50 length was 2 130 bp； after redundancy removal， 32 855 transcripts were obtained. Among the four major databases， Nr， KEGG， KOG and Swiss-Prot， 29 404， 27 794， 20 286 and 21 197 transcripts were annotated respectively； among them， 19 346 transcripts were annotated in all four databases， and 29 416 transcripts were annotated in at least one database， accounting for 89.53% of the total. In addition， 1 139 transcription factors （TFs）， 2 636 long-stranded non-coding RNAs（LncRNAs） and 4 197 SSR were identified， and 32 800 CDS were obtained after prediction.

Key words： Cryptolestes ferrugineus， full-length transcriptome， genome annotation

銹赤扁谷盜Cryptolestes ferrugineus（Stephens）屬于鞘翅目扁谷盜科，主要分布在溫帶和熱帶地區(qū)。在全球范圍內(nèi)，銹赤扁谷盜是發(fā)生最嚴重、最廣泛的儲糧害蟲之一，一旦發(fā)生便會造成嚴重的不可逆損失[1-2]。銹赤扁谷盜在我國主要分布于糧食生態(tài)區(qū)，在廣東、云南和海南等地都有發(fā)生，常見于米廠、面粉廠、酒廠和一些糧倉內(nèi)[3]。銹赤扁谷盜在糧堆的分布也受溫度和濕度的影響[4]，在糧堆發(fā)熱的情況下，更容易造成銹赤扁谷盜的繁殖，嚴重時可出現(xiàn)糧堆結(jié)露，進而使糧食發(fā)生霉變，給儲糧安全帶來極大的威脅[5]。

盡管二代測序高通量測序技術(shù)對轉(zhuǎn)錄組學(xué)的相關(guān)研究起到了極大的推動作用，但還是存在測序讀長較短[6]、重復(fù)區(qū)域拼接效果差、無法得到較完整的轉(zhuǎn)錄本[7]等局限性。近年來，與第二代測序相比，三代轉(zhuǎn)錄組測序可快速全面地獲得某一物種的全長轉(zhuǎn)錄本信息，逐漸被應(yīng)用于昆蟲學(xué)研究中，例如，按蚊[8]、家蠶[9]、小菜蛾[10]等，主要集中在生長發(fā)育抗藥性和遺傳變異等方面。在銹赤扁谷盜的研究中，關(guān)于三代測序的報道相對較少，因此，本文對銹赤扁谷盜進行全長轉(zhuǎn)錄組測序并通過基因功能注釋、CDS預(yù)測、TFs分析等，進一步從分子水平對銹赤扁谷盜的研究奠定基礎(chǔ)，也為探究新的靶標基因、利用分子設(shè)計實現(xiàn)害蟲防治提供理論依據(jù)并豐富鞘翅目昆蟲的全長轉(zhuǎn)錄組數(shù)據(jù)庫。

1 材料與方法

1.1 供試昆蟲

本試驗銹赤扁谷盜種群采自上海福星面粉廠，于南京財經(jīng)大學(xué)糧食儲運國家工程實驗室培養(yǎng)。在實驗室內(nèi)以人工飼料（全麥粉∶燕麥∶小麥碎∶酵母粉=5∶3∶3∶1）、恒溫29～31 ℃、相對濕度為70%～80% RH、無光照條件下飼養(yǎng)。在培養(yǎng)過程中，銹赤扁谷盜以3 d產(chǎn)的卵為同一齡期，每隔3 d將瓶中銹赤扁谷盜成蟲全部挑出放入相同飼料的新瓶子中，重復(fù)此操作4～5次之后即可獲得齡期一致的試蟲。

1.2 方法

1.2.1 總RNA提取及文庫構(gòu)建

對銹赤扁谷盜幼蟲、蛹、成蟲3個階段分別進行取樣，每個階段3個重復(fù)，對混合樣品進行RNA提取和富集，使用NanoDrop分光光度計檢測樣品RNA的濃度；使用Agilent2100檢測樣品RNA的完整度?？俁NA（含有poly（A））檢測合格后，利用Clontech SMARTer PCR cDNA Synthesis Kit反轉(zhuǎn)錄成第一鏈cDNA，PCR擴增合成雙鏈cDNA并對PCR產(chǎn)物進行純化，對cDNA進行DNA損傷修復(fù)、末端修復(fù)、連接Adapter和文庫質(zhì)量評估工作后形成完整的SMRT bell文庫，再對文庫進行測序。

1.2.2 全長轉(zhuǎn)錄組數(shù)據(jù)分析和校正

使用SMRT Link V8.0原始數(shù)據(jù)進行分析。首先從下機數(shù)據(jù)中提取高質(zhì)量的CCS（circular consensus sequence）序列，移除引物和barcode、poly（A）和連環(huán)結(jié)構(gòu)，得到全長非嵌合序列（FLNC），再對相似的FLNC reads進行聚類，合并成一個完整的轉(zhuǎn)錄本序列。使用LoRDEC軟件進行轉(zhuǎn)錄本校正，LoRDEC采用混合策略，需要使用兩組數(shù)據(jù)：參考reads（二代短reads）、PacBio長reads。通過讀取短reads，建立DBG圖，針對長reads中的錯誤區(qū)域?qū)ふ倚Ｕ蛄?。之后，通過Nr、KEGG、KOG、Swiss-Prot、Pfam和GO數(shù)據(jù)庫分別對銹赤扁谷盜全長轉(zhuǎn)錄組進行功能注釋，并進行轉(zhuǎn)錄因子（TFs）鑒定、長鏈非編碼RNA（lncRNAs）預(yù)測和SSR分析等。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)分析

測序結(jié)果表明，銹赤扁谷盜全長轉(zhuǎn)錄組數(shù)據(jù)庫包含414 552條CCSs，平均長度為1 868 bp。此外，結(jié)果顯示含有12 009 936條子序列，子序列平均長度為1 845 bp，N50長度為2 094 bp。在聚類和校正之后，對序列進行去冗余，獲得32 855條轉(zhuǎn)錄本，平均長度為1 845.06 bp（如圖1），N50長度為2 130 bp，如表1。

2.2 功能注釋

2.2.1 基本注釋

對得到的32 855條轉(zhuǎn)錄本進行Nr、KEGG、KOG、Swiss-Prot四大數(shù)據(jù)庫注釋，如圖2所示，其中Nr數(shù)據(jù)庫中有29 404條轉(zhuǎn)錄本被注釋，KEGG數(shù)據(jù)庫中有27 794條轉(zhuǎn)錄本被注釋，KOG數(shù)據(jù)庫中有20 286條轉(zhuǎn)錄本被注釋，SwissProt數(shù)據(jù)庫中有21 197條轉(zhuǎn)錄本被注釋，其中至少被一個數(shù)據(jù)庫注釋的轉(zhuǎn)錄本數(shù)目為29 416條（89.53%），19 346條轉(zhuǎn)錄本在四大數(shù)據(jù)庫中都有注釋。

2.2.2 E值分布

將所有轉(zhuǎn)錄本在 Nr、Swiss-Prot、KEGG 和KOG 四大數(shù)據(jù)庫中的最佳比對結(jié)果的E值進行統(tǒng)計，將其分為5個范圍。由E值可以看出，轉(zhuǎn)錄與數(shù)據(jù)庫中匹配序列為同源序列的假陽性概率。如圖3所示，Nr數(shù)據(jù)庫中，在0≤E值≤1E-150范圍內(nèi)有16 004條轉(zhuǎn)錄本，占54.43%；SwissProt數(shù)據(jù)庫中，在0≤E值≤1E-150范圍內(nèi)有4 843條轉(zhuǎn)錄本，占22.85%；KEGG數(shù)據(jù)庫中，在0≤E值≤1E-150范圍內(nèi)有11 321條轉(zhuǎn)錄本，占40.73%；KOG 數(shù)庫中，在0≤E值≤1E-150范圍內(nèi)有5 454條轉(zhuǎn)錄本，占26.89%。

2.2.3 物種分布

利用Blastx將轉(zhuǎn)錄本序列與Nr數(shù)據(jù)庫進行比對后，取每個轉(zhuǎn)錄本在Nr數(shù)據(jù)庫中比對結(jié)果最好（E值最低）的那一條序列作為對應(yīng)同源序列確定同源序列所屬物種，統(tǒng)計比對到各個物種的同源序列數(shù)量，同源性較高的前3位分別是光肩星天牛、赤擬谷盜和蜂箱小甲蟲，如表2所示。

2.2.4 GO分類

在GO數(shù)據(jù)庫中獲得注釋轉(zhuǎn)錄本173 206條，如圖4所示，生物學(xué)過程（biological process）包括25個功能組，其中細胞過程（cellular process）包含轉(zhuǎn)錄本最多，有13 585條，其次是單一生物過程（single-organism process）12 740條，最少的是生物階段（biological phase），只有11條；細胞組分（cellular component）包括21個功能組，其中細胞（cell）和細胞組分（cell part）包含轉(zhuǎn)錄本最多，都有11 053條，最少的是擬核（nucleoid），只有1條；分子功能（molecular function）包括12個分子功能，其中結(jié)合活性（binding）最多，有10 176條，其次是催化活性（catalytic activity）9 217條，電子載體活性（electron carrier activity）最少，只有14條。

2.3 高級注釋及基因結(jié)構(gòu)分析

2.3.1 TFs鑒定

將預(yù)測的蛋白序列同相應(yīng)的TF數(shù)據(jù)庫（Animal TFdb）進行 hmmscan 比對。結(jié)果顯示，一共預(yù)測到1 139個轉(zhuǎn)錄因子，屬于55個TF家族，對轉(zhuǎn)錄本數(shù)目最多的前10個TF家族進行繪圖，如圖5所示。其中Zf-C2H2家族有355個、ZBTB家族有74個、TF_bzip家族有72個，HMG家族有71個。

2.3.2 LncRNA預(yù)測

使用CNCI軟件和CPC軟件進行編碼能力預(yù)測，取兩個軟件都預(yù)測為“非編碼”的結(jié)果作為最終的lncRNA結(jié)果。如圖6所示，CNCI預(yù)測lncRNA 2 871個，CPC預(yù)測lncRNA 2 846個，兩個軟件同時預(yù)測lncRNA 2 636個。

2.3.3 SSR分析

利用MISA軟件對轉(zhuǎn)錄組的所有轉(zhuǎn)錄本進行搜索，尋找轉(zhuǎn)錄本中的SSR，分析結(jié)果顯示，不同串聯(lián)重復(fù)單元類型的SSR在總SSR中所占比例不同，其中AAT/ATT的頻率最高，所占比例為31.3%。如圖7所示，在SSR核苷酸基序類型中，三核苷酸重復(fù)最多（tri nucleotides，69.19%），依次是四核苷酸重復(fù)（tetra nucleotides，5.23%）、雙核苷酸重復(fù)（di nucleotide motifs，13.27%）、五核苷酸重復(fù)（penta nucleotides，2.14%）和六核苷酸重復(fù)（hexa-nucleotide motifs，0.17%）。

3 討論與結(jié)論

二代測序雖在轉(zhuǎn)錄組學(xué)相關(guān)研究中應(yīng)用廣泛，但存在測序讀長短，對高重復(fù)區(qū)域無法較好地拼接等缺點。相比于二代轉(zhuǎn)錄組測序定量的特點，全長轉(zhuǎn)錄組還可對轉(zhuǎn)錄本進行定性分析，即分析轉(zhuǎn)錄本的結(jié)構(gòu)。第三代測序技術(shù)具有超長讀長、無需組裝即可直接獲得RNA的全長序列，還可鑒定基因的可變剪切、可變多聚核苷酸化APA、非編碼RNA等，讓測序的轉(zhuǎn)錄組數(shù)據(jù)更加豐富可靠。

本研究利用第三代測序技術(shù)對銹赤扁谷盜不同齡期混合樣品進行測序，校正后去冗余共獲得32 855條轉(zhuǎn)錄本，平均長度為1 845.06 bp，N50長度為2 130 bp。通過在Nr、KEGG、KOG、Swiss-Prot四大數(shù)據(jù)庫中注釋，至少被一個數(shù)據(jù)庫注釋的轉(zhuǎn)錄本數(shù)目為29 416條，占總轉(zhuǎn)錄本的89.53%，被注釋到Nr數(shù)據(jù)庫的最多，有29 404條。而對蓮草直胸跳甲卵、幼蟲、蛹、成蟲4個發(fā)育階段的樣本的全長轉(zhuǎn)錄組測序，共獲得28 982條高質(zhì)量的轉(zhuǎn)錄本，并預(yù)測得到4 198個lncRNA，24 040個開放閱讀框ORF[11]；對家蠶絲腺的測序獲得11 697條轉(zhuǎn)錄本[12]。lncRNA具有調(diào)控轉(zhuǎn)錄、剪切或相鄰基因的表達等功能。例如在果蠅中，lncRNA可以調(diào)控雄性果蠅的性別決定[13]、睡眠[14]、運動行為[15]等生物過程；在家蠶絲腺中發(fā)現(xiàn)了一些與蠶絲基因表達相關(guān)的lncRNA[16]。近年來，大量的lncRNA已從家蠶[16]、黑腹果蠅[17]、小菜蛾、岡比亞按蚊、中華蜜蜂[18]等昆蟲中鑒定出來。本文共獲得2 636個lncRNA，為后續(xù)在基因表達和調(diào)控提供參考。在銹赤扁谷盜中共獲得1 139個轉(zhuǎn)錄因子，Zf-C2H2家族有355個、ZBTB家族有74個、TF_bzip家族有72個，HMG家族有71個。不同的TFs可能參與不同的代謝過程，也可能有多種不同的功能[19]。在很多研究中發(fā)現(xiàn)，轉(zhuǎn)錄因子與解毒酶基因的表達相關(guān)。例如，在赤擬谷盜中，轉(zhuǎn)錄因子CncC和Maf通過調(diào)控CYP6BQ基因上調(diào)增強其對溴氰菊酯的代謝[20]；轉(zhuǎn)錄因子FTZ-F1可以調(diào)控小菜蛾CYP6BG1基因的表達，從而使小菜蛾對氯蟲苯酰胺的代謝增強[21]。所以，對銹赤扁谷盜的TFs分析會對其免疫和抗藥性方面的分析提供依據(jù)。通過對SSR開發(fā)標記發(fā)掘其功能，也為銹赤扁谷盜遺傳多樣性分析提供更多數(shù)據(jù)支持。

綜上，本文基于全長轉(zhuǎn)錄組測序技術(shù)獲得了銹赤扁谷盜的全長轉(zhuǎn)錄組，并注釋了銹赤扁谷盜全長轉(zhuǎn)錄組相關(guān)基因功能信息，也為進一步研究銹赤扁谷盜的防治篩選靶標，并為其遺傳機制的研究奠定基礎(chǔ)。

參 考 文 獻

[1]PHILLIPS T W， THRONE J E. Biorational approaches to managing stored-product insects[J]. Annual review of entomology， 2010， 55： 375-397.

[2] LOSEY S M， DAGLISH G J， PHILLIPS T W. Orientation of rusty grain beetles， Cryptolestes ferrugineus （Coleoptera：Laemophloeidae）， to semiochemicals in field and laboratory experiments[J]. Journal of Stored Products Research， 2019， 84： 101513.

[3] 張瑞杰，王殿軒，白春啟，等.中國12省80地市儲糧中扁谷盜屬昆蟲分布調(diào)查[J].河南工業(yè)大學(xué)學(xué)報（自然科學(xué)版），2018， 39（5）：112-118.

[4] JIAN F J， LARSON R， JAYAS D S， et al. Three-dimensional spatial distribution of adults of Cryptolestes ferrugineus （Coleoptera： Laemophloeidae） in stored wheat under different temperatures， moisture contents， and adult densities[J]. Journal of Stored Products Research， 2011， 47（4）： 293-305.

[5] PARDE S R， JAYAS D S， WHITE N D G. Movement of Cryptolestes ferrugineus （Coleoptera： Cucujidae） in grain columns containing pockets of high moisture content wheat and carbon dioxide gradients[J]. Journal of Stored Products Research， 2004， 40（3）： 299-316.

[6] KOREN S， SCHATZ M C， WALENZ B P， et al. Hybrid error correction and de novo assembly of single-molecule sequencing reads[J]. Nature Biotechnology， 2012， 30（7）： 693-700.

[7]田李，張穎，趙云峰.新一代測序技術(shù)的發(fā)展和應(yīng)用[J].生物技術(shù)通報，2015，31（11）：1-8.

[8] JIANG X， HALL A B， BIEDLER J K， et al. Single molecule RNA sequencing uncovers trans‐splicing and improves annotations in Anopheles stephensi[J]. Insect Molecular Biology， 2017， 26（3）： 298-307.

[9] KAWAMOTO M， JOURAKU A， TOYODA A， et al. High-quality genome assembly of the silkworm， Bombyx mori [J]. Insect Biochemistry and Molecular Biology， 2019， 107： 53-62.

[10] ZHAO Q， ZHONG W M， HE W Y， et al. Genome-wide profiling of the alternative splicing provides insights into development in Plutella xylostella [J]. BMC Genomics， 2019， 20（1）： 1-14.

[11] JIA D， WANG Y X， LIU Y H， et al. SMRT sequencing of fulllength transcriptome of flea beetle Agasicles hygrophila （Selman and Vogt）[J]. Scientific Reports， 2018， 8（1）： 1-8.

[12] CHEN T， SUN Q W， MA Y， et al. A transcriptome atlas of silkworm silk glands revealed by PacBio single-molecule long-read sequencing[J]. Molecular Genetics and Genomics， 2020， 295（5）： 1227-1237.

[13] MULVEY B B， OLCESE U， CABRERA J R， et al. An interactive network of long non-coding RNAs facilitates the Drosophila sex determination decision[J]. Biochimica et Biophysica Acta （BBA）-Gene Regulatory Mechanisms， 2014， 1839（9）： 773-784.

[14] SOSHNEV A A， ISHIMOTO H， MCALLISTER B F， et al. A conserved long noncoding RNA affects sleep behavior in Drosophila[J]. Genetics， 2011， 189（2）： 455-468.

[15] LI M X， WEN S Y， GUO X Q， et al. The novel long non-coding RNA CRG regulates Drosophila locomotor behavior[J]. Nucleic Acids Research， 2012， 40（22）： 11714-11727.

[16] WU Y Q， CHENG T C， LIU C， et al. Systematic identification and characterization of long non-coding RNAs in the silkworm， Bombyx mori[J]. PloS one， 2016， 11（1）： e0147147.

[17] YOUNG R S， MARQUES A C， TIBBIT C， et al. Identification and properties of 1，119 candidate lincRNA loci in the Drosophila melanogaster genome[J]. Genome Biology and Evolution， 2012， 4（4）： 427-442.

[18] JAYAKODI M， JUNG J W， PARK D， et al. Genome-wide characterization of long intergenic non-coding RNAs （lincRNAs） provides new insight into viral diseases in honey bees Apis cerana and Apis mellifera[J]. BMC Genomics， 2015， 16（1）： 1-12.

[19] CHEN K， RAJEWSKY N. The evolution of gene regulation by transcription factors and microRNAs[J]. Nature Reviews Genetics， 2007， 8（2）： 93-103.

[20] KALSI M， PALLI S R. Transcription factors， CncC and Maf， regulate expression of CYP6BQ genes responsible for deltamethrin resistance in Tribolium castaneum[J]. Insect Biochemistry and Molecular Biology， 2015， 65： 47-56.

[21] LI X X， SHAN C Y， LI F， et al. Transcription factor FTZ‐F1 and cis‐acting elements mediate expression of CYP6BG1 conferring resistance to chlorantraniliprole in Plutella xylostella[J]. Pest Management Science， 2019， 75（4）： 1172-1180.