耿合員 何 江 王 淵 孫立新 聶國(guó)嬡 胡雙雙 楊慶貴**

(1.江蘇省檢驗(yàn)檢疫科學(xué)技術(shù)研究院, 江蘇南京 210019; 2.江蘇國(guó)際旅行衛(wèi)生保健中心, 江蘇南京 210019; 3.陜西國(guó)際旅行衛(wèi)生保健中心, 陜西西安 710068)

中華按蚊Anophelessinensis是嚴(yán)重威脅人類健康的病媒生物，野外種群數(shù)量大，滋生環(huán)境復(fù)雜，適應(yīng)能力強(qiáng)，能夠通過(guò)叮咬宿主傳播包括瘧疾、馬來(lái)絲蟲病、流行性乙型腦炎等多種傳染病，因此長(zhǎng)期以來(lái)是預(yù)防醫(yī)學(xué)的重點(diǎn)監(jiān)測(cè)防控對(duì)象(Zhangetal., 2015; Fengetal., 2017; Fangetal., 2018)。隨著測(cè)序技術(shù)應(yīng)用的普及，以轉(zhuǎn)錄組測(cè)序?yàn)榇淼幕蜓芯恐饾u成為中華按蚊等害蟲的主要研究手段之一。目前常規(guī)意義上的轉(zhuǎn)錄組測(cè)序是指基于Illumina等短讀長(zhǎng)平臺(tái)對(duì)樣本mRNA的序列獲取，結(jié)果分析包含不同基因mRNA的數(shù)量、類型與表達(dá)量等信息?；谵D(zhuǎn)錄組數(shù)據(jù)庫(kù)的比較能夠解析目標(biāo)物種包括不同發(fā)育階段、生理反應(yīng)、脅迫應(yīng)答等條件差異下參與的基因及其可能的作用機(jī)制。對(duì)于中華按蚊基于不同時(shí)空與發(fā)育、脅迫條件下的比較轉(zhuǎn)錄組分析已在此物種吸血分子基礎(chǔ)、消化代謝與免疫以及抗藥性機(jī)理等方面取得了大量進(jìn)展，為深入理解其生理生態(tài)基因表達(dá)調(diào)控機(jī)制提供證據(jù)支持(Sunetal., 2017; Liuetal., 2018; Yanetal., 2018; Zhouetal., 2018)。但是，轉(zhuǎn)錄組測(cè)序技術(shù)的成熟也伴隨著新問(wèn)題的出現(xiàn)，其中樣本不同采集與保存方法所導(dǎo)致的差異凸顯，成為了影響測(cè)序與分析結(jié)果的重要因素(Shenetal., 2018; Sonetal., 2018)。以中華按蚊為例，其個(gè)體藏身隱蔽且十分脆弱，野外活體采集具備難度，部分樣本采集完拿到實(shí)驗(yàn)室時(shí)已死。而且，由于野外采集的實(shí)驗(yàn)條件通常有限，因此樣本在送達(dá)實(shí)驗(yàn)室安全保存前通常也會(huì)發(fā)生不同程度的RNA降解。此外，中華按蚊體內(nèi)常具有已吸取的宿主血液，從而對(duì)樣本的RNA提取造成干擾。這些干擾因素均會(huì)影響樣本構(gòu)建測(cè)序文庫(kù)中mRNA的豐度，降低后續(xù)分析的準(zhǔn)確性。

鑒于此，本研究擬對(duì)中華按蚊在吸血、存活狀態(tài)及投入RNA貯存液前的冷凍時(shí)間等變量下轉(zhuǎn)錄組測(cè)序結(jié)果差異進(jìn)行比較分析，分析這些變量下組裝的完整性，鑒定檢出基因的組成及表達(dá)量差異，以及這部分差異基因主要的功能。本研究結(jié)果將從方法學(xué)角度，為之后中華按蚊的轉(zhuǎn)錄組測(cè)序研究在樣本采集與保存方面的策略優(yōu)化，以及有效評(píng)估現(xiàn)有轉(zhuǎn)錄組數(shù)據(jù)有效性等方面，提供詳盡數(shù)據(jù)支持。

1 材料與方法

1.1 樣本采集

中華按蚊樣本采集自云南瑞麗一養(yǎng)牛場(chǎng)。依據(jù)本研究的目的，共采集4類樣本，樣本均為雌性，分別為未吸血活體、吸血活體、吸血活體在-20 ℃冷凍20 min、吸血死樣本在-20 ℃冷凍20 h。詳細(xì)的樣本分組如表1所示。

表1 中華按蚊轉(zhuǎn)錄組測(cè)序樣本信息

1.2 RNA樣本提取與轉(zhuǎn)錄組測(cè)序

使用天根試劑盒提取各組樣本的總RNA，操作方法參考試劑盒說(shuō)明書。每組選擇5個(gè)個(gè)體，混樣提取，提取出的RNA樣本使用Bioanalyzer平臺(tái)進(jìn)行檢測(cè)合格后，用于轉(zhuǎn)錄組測(cè)序。轉(zhuǎn)錄組文庫(kù)構(gòu)建與測(cè)序委托深圳華大基因公司完成，采用BGISEQ-500平臺(tái)，每組3個(gè)生物學(xué)重復(fù)，4組共計(jì)12個(gè)樣本，每個(gè)樣本測(cè)序至少6 Gb。

1.3 數(shù)據(jù)處理與分析

測(cè)序原始數(shù)據(jù)使用Trimmomatic v0.33軟件(Bolgeretal., 2014)過(guò)濾包含接頭、未知堿基N含量大于5%等低質(zhì)量序列。所有過(guò)濾后的清潔標(biāo)簽(Clean reads)用于之后輸入Trinity v2.11.0軟件(Grabherretal., 2011)合并組裝，然后使用Tgicl v2.1軟件(Perteaetal., 2003)對(duì)轉(zhuǎn)錄本進(jìn)行聚類去冗余得到最終轉(zhuǎn)錄組基因數(shù)據(jù)集Unigene。Unigene的質(zhì)量使用BUSCO v3(Benchmarking Universal Single-Copy Orthologs)軟件(Seppeyetal., 2019)中單拷貝直系同源基因(Single-Copy Orthologs, SCO)完整性比例進(jìn)行評(píng)估，評(píng)估采用昆蟲(insecta)共有SCO基因數(shù)據(jù)庫(kù)。使用TransDecoder v5.3.0軟件(Grabherretal., 2011)識(shí)別Unigene中的編碼區(qū)域(Coding domain sequence, CDS)，識(shí)別結(jié)合Blast 2.2.28軟件(Camachoetal., 2009)比對(duì)SwissProt數(shù)據(jù)庫(kù)(McMillanetal., 2008)和Hmmerscan軟件(Chaietal., 2019)搜索Pfam 數(shù)據(jù)庫(kù)中(Finnetal., 2014)蛋白同源性的結(jié)果，以獲得更完整編碼區(qū)域。

功能注釋通過(guò)將組裝得到的Unigene在常用功能數(shù)據(jù)庫(kù)中依據(jù)序列相似性比對(duì)獲得KEGG(Kanehisaetal., 2017)、GO(The gene ontology resource)、NR、NT(Bensonetal., 2018)、SwissProt(McMillanetal., 2008)、Pfam(Finnetal., 2014)和KOG(Mudado Mdeetal., 2006)，此過(guò)程中采用的軟件分別為KAAS(Moriyaetal., 2007)(KEGG)、Blast2GO basic(Conesaetal., 2005)(GO)、Blast v2.2.28(Camachoetal., 2009)(NR、NT、SwissProt、KOG)、Hmmscan v3.2.1(Chaietal., 2019)(Pfam)。

使用 RSEM v1軟件(Lietal., 2011)將各個(gè)樣品與Unigene進(jìn)行比對(duì)，計(jì)算基因表達(dá)水平。R中deseq軟件包(Andersetal., 2010)被用于計(jì)算不同樣本組(重復(fù)樣品歸為一組)的差異表達(dá)量分析，使用每百萬(wàn)堿基檢出片段數(shù)量(Fragments per Kilobase per million，F(xiàn)PKM)值作為衡量表達(dá)量差異的參數(shù)。差異倍數(shù)達(dá)到兩倍以上并且二輪校準(zhǔn)后的P-value(Q-value)≤ 0.001 被用來(lái)篩選顯著差異表達(dá)基因。差異基因的GO、KEGG富集分析使用R軟件中的phyper軟件包(Berlivetetal., 2007)，Q-value≤ 0.05的功能視為顯著富集。結(jié)果的可視化通過(guò)R軟件中的GG-plot2軟件包(Maag, 2018)完成。

2 結(jié)果

2.1 轉(zhuǎn)錄組組裝結(jié)果

轉(zhuǎn)錄組組裝結(jié)果所得的Unigene基因數(shù)據(jù)集及其相關(guān)功能注釋信息是之后比較分析的基礎(chǔ)。本研究所選擇的12個(gè)中華按蚊樣本共測(cè)得75.58 Gb 數(shù)據(jù)，平均每個(gè)樣本測(cè)序6.30 Gb。最終去冗余從頭組裝結(jié)果包含48 775個(gè)Unigene。序列總大小為 68.1 Mb，平均長(zhǎng)度為1 395 bp、N50值為2 349 bp，整體GC含量為47.75 %。進(jìn)一步從這些Unigene中檢測(cè)出28 265個(gè)編碼序列CDS，同時(shí)預(yù)測(cè)出5 263個(gè)編碼轉(zhuǎn)錄因子。綜合統(tǒng)計(jì)Unigene比對(duì)到不同基因功能數(shù)據(jù)庫(kù)的結(jié)果，顯示共有35 567(72.92%)個(gè)基因能夠獲得功能注釋信息。其中各數(shù)據(jù)庫(kù)分別注釋了33 330(NR：68.33%)、19 966(NT：40.93%)、25 033(SwissProt：51.32%)、24 512(KOG：50.26%)、26 293(KEGG：53.91%)、26 704(GO：54.75%)以及25 756(Pfam：52.81%)個(gè)基因。

基于統(tǒng)一的Unigene數(shù)據(jù)集，研究獲取了各樣本中所包含的基因及其表達(dá)量。表達(dá)量定量結(jié)果顯示，各組間全部基因表達(dá)量分布范圍接近，log10(FPKM+1)的值主體介于0到1之間，離散點(diǎn)較少(圖1)?；虮磉_(dá)密度統(tǒng)計(jì)顯示各樣本表達(dá)分布高度重疊，只存在一個(gè)單一的峰(圖2)。這些結(jié)果說(shuō)明構(gòu)建的Unigene數(shù)據(jù)集在解析各樣本時(shí)的結(jié)果可靠，可用于后續(xù)比較分析。

圖1 表達(dá)量箱線圖

圖2 表達(dá)量密度圖

2.2 中華按蚊轉(zhuǎn)錄組組間組裝質(zhì)量比較

研究首先比較了各組樣本相對(duì)于整體的組裝完整性。結(jié)果顯示，整體Unigene的SCO完整性極高，達(dá)到了98.7%。各組間比較結(jié)果顯示全部樣本的完整性均低于整體，即各組均有部分基因未能檢出。之后，組間比較也發(fā)現(xiàn)明顯差異，其中M1組與M4組的SCO完整性低于M2組與M3組，而M1組2號(hào)重復(fù)樣本與M4組3號(hào)重復(fù)樣本的完整SCO基因比例甚至低于60%，存在大量片段及未能檢出的缺失的SCO基因(圖3)。該結(jié)果說(shuō)明至少?gòu)腟CO抽樣水平，各組檢出的基因數(shù)量及質(zhì)量有區(qū)別。

圖3 組裝結(jié)果SCO基因完整性示意圖

2.3 組間差異基因鑒定

之后進(jìn)一步從全轉(zhuǎn)錄組完整Unigene水平鑒定并比較了各組檢出基因的差異。結(jié)果顯示，四組共有的Unigene為 28 183條，僅占全部的57.78%(圖4)。從韋恩圖上可知，M3與M4組中檢出基因較少是導(dǎo)致共有基因占比較低的原因，二者分別只檢出33 267與33 058個(gè)基因。而M1與M2組檢出基因數(shù)量更多(分別為45 887 與44 478個(gè))，二者共有42 584個(gè)基因，占全部的87.31%(圖4)。進(jìn)一步，為了揭示M1與M2組所多檢出的基因的功能分布，本研究對(duì)在M3與M4組中均未能檢出而M1與M2組均檢出的11 358條基因的功能進(jìn)行了富集。結(jié)果顯示，GO與KEGG富集的通路均與基礎(chǔ)代謝機(jī)制相關(guān)(圖5～6)。如GO中最顯著富集基因?yàn)榘l(fā)生在核內(nèi)(Nucleus)的核酸結(jié)合(Nucleic acid binding)與鋅離子結(jié)合(Zinc ion binding)，為基礎(chǔ)代謝中最常應(yīng)用的分子功能(圖5)。而KEGG中最顯著富集的基因主要參與了核糖體生物合成(Ribosome biogenesis in eukaryotes)、谷胱甘肽生物合成(Glycosaminoglycan biosynthesis)、基礎(chǔ)轉(zhuǎn)錄因子(Basal transcription factors)等基本物質(zhì)代謝途徑(圖6)。

圖4 轉(zhuǎn)錄組組間基因韋恩圖

圖5 GO類群功能富集圖

圖6 KEGG通路富集圖

2.4 組間共有基因差異表達(dá)水平

對(duì)于所有組均檢出的基因，研究比較了其表達(dá)量以分析不同變量在解析表達(dá)水平層面上的影響。結(jié)果顯示，任一兩組間均存在大量差異基因，發(fā)生顯著上調(diào)與下調(diào)的基因總數(shù)均超過(guò)10 000 條，超過(guò)全部共有基因(28 183條)的30%。具體而言，研究發(fā)現(xiàn)M3與M4組間的差異表達(dá)基因數(shù)量最少，且上下調(diào)接近；M1與M2組相比上調(diào)基因多于下調(diào)，而其余組比較均為下調(diào)基因多于上調(diào)基因；M1與M3組發(fā)生顯著調(diào)控的基因數(shù)量最多，且上調(diào)與下調(diào)基因數(shù)量均為各比較組中最高(圖7)。

圖7 差異表達(dá)基因統(tǒng)計(jì)圖

通過(guò)對(duì)各組差異表達(dá)基因分別進(jìn)行GO功能富集以比較其組成與功能上的差異，發(fā)現(xiàn)這些基因幾乎不存在組間功能重疊，表現(xiàn)出強(qiáng)組間特異性特征(圖8)。對(duì)于M2組與M3組而言，其差異表達(dá)基因的作用位置甚至都與其他組不同，該組基因主要作用于胞外區(qū)域(extracellular region)，而其他組差異表達(dá)基因則主要作用于核內(nèi)。該組中與能量代謝相關(guān)的通路(ATP synthesis)也同樣被富集出來(lái)。M1組與M2組中與細(xì)胞代謝活性相關(guān)的通路發(fā)生顯著富集，如DNA復(fù)制(DNA replication)與RNA結(jié)合(RNA binding)等。M3與M4組比較，沒(méi)有獲得有實(shí)際生物學(xué)意義的GO分類。M1與M4組相比，富集了多個(gè)發(fā)揮免疫與解毒相關(guān)的分子功能，包括內(nèi)肽酶活性(Endopeptidase activity)與去氧化活性(Oxidoreductase activity)等。這些功能均反映著各比較組中所對(duì)應(yīng)的變量對(duì)基因表達(dá)量的影響，提示不同的保存與處理方法對(duì)于樣本RNA質(zhì)量的干擾明顯。

圖8 差異表達(dá)基因GO類群富集圖

3 討論

樣本的RNA是轉(zhuǎn)錄組測(cè)序的物質(zhì)基礎(chǔ)，其質(zhì)量決定了結(jié)果的可靠性。在本研究中，我們以中華按蚊這一需監(jiān)控的且具有復(fù)雜樣本采集與保存過(guò)程的病媒生物為材料在4種不同條件變量(非吸血活體、吸血活體、吸血活體短期冷凍、吸血死體長(zhǎng)時(shí)間冷凍)下的轉(zhuǎn)錄組測(cè)序結(jié)果差異進(jìn)行了分析與比較。本研究依次從轉(zhuǎn)錄組組裝水平、檢出核心直系同源基因比例、整體檢出基因數(shù)量差異及功能富集、共有基因的表達(dá)差異及功能富集等4個(gè)層面進(jìn)行了詳盡的量化評(píng)估。

經(jīng)結(jié)果整合，研究發(fā)現(xiàn)吸血活體為最優(yōu)樣本，其3個(gè)重復(fù)樣本的核心基因完整性均超過(guò)90%，檢出基因總數(shù)占完整Unigene數(shù)據(jù)庫(kù)的比例也達(dá)到了44 478個(gè)，占全部的91.19%。這與研究預(yù)期的非吸血活體樣本為最優(yōu)存在一定偏差，即吸血樣本中存有的宿主血液中的核酸對(duì)轉(zhuǎn)錄組測(cè)序樣本的干擾影響十分微弱。進(jìn)一步推測(cè)其背后的原因，研究首先排除了多檢出基因?yàn)樗拗饕氲目赡苄裕次铙w與非吸血樣本相比較其檢出的基因具有高度重疊，之后基于兩組間差異表達(dá)基因的功能富集結(jié)果主要為能量代謝相關(guān)，推測(cè)吸血樣本具有更活躍的代謝水平可能導(dǎo)致其基因表達(dá)豐度更高，從而更容易被轉(zhuǎn)錄組測(cè)序平臺(tái)捕獲。之前認(rèn)為的宿主外源污染，其可能在樣本捕獲時(shí)便已降解，或在之后組裝前的數(shù)據(jù)過(guò)濾階段被清除，對(duì)測(cè)序結(jié)果不構(gòu)成實(shí)際影響。

與吸血組相比，非吸血活體樣本檢出基因數(shù)量最多，達(dá)到45 887個(gè)，占總數(shù)的94.01%。但是，該樣本的核心直系同源基因完整性卻極低，平均僅為63.3%。從核心基因的構(gòu)成看，推測(cè)該樣本可能存在組裝問(wèn)題，即測(cè)序?qū)虻母采w不完整，呈現(xiàn)出碎片化特征。導(dǎo)致該情況出現(xiàn)的原因從樣本生理狀態(tài)的角度較難找到合理解釋，因?yàn)槿绻莻€(gè)體虛弱或降解的原因，其差異表達(dá)基因與檢出基因數(shù)量應(yīng)該均明顯少于目前鑒定的狀況。于是，研究只能暫且推測(cè)該樣本是在建庫(kù)過(guò)程中mRNA捕獲時(shí)便出現(xiàn)了缺失，具體根源還有待于之后通過(guò)其他測(cè)序技術(shù)方法進(jìn)行探究。

對(duì)于樣本的捕獲狀態(tài)，研究發(fā)現(xiàn)采集時(shí)立即處死對(duì)于轉(zhuǎn)錄組測(cè)序的樣本影響較小。與活體樣本相比，死亡個(gè)體樣本檢出基因最少，且具有大量差異表達(dá)基因，即存在明顯降解。然而，出乎意料的是，就核心基因完整性而言，死亡樣本并未發(fā)生斷崖式下降，其與非吸血活體樣本接近，存在大量片段化基因。除此之外，在檢出的基因中，其與活體短期冷凍20 min樣本間具有廣泛重疊，而與其他活體組鑒定出的差異表達(dá)基因中，活體樣本中發(fā)生下調(diào)的基因數(shù)量均多于上調(diào)，即死體檢出基因表達(dá)表現(xiàn)出受刺激應(yīng)答反應(yīng)的特征。這些結(jié)果與通常認(rèn)為的死亡個(gè)體后mRNA會(huì)快速完全降解的認(rèn)知不符，推測(cè)是由于本研究所選擇的個(gè)體是采集時(shí)才處死的，采集之前此個(gè)體的躲避行為可能導(dǎo)致了其表現(xiàn)出刺激應(yīng)答特征。該推測(cè)也一定程度上得到了差異基因功能富集的支持，本研究中的死體組與其他吸血組相比光傳導(dǎo)(Phototransduction)均被顯著富集，與非吸血組相比，則有多個(gè)解毒與消化相關(guān)的基因發(fā)生表達(dá)改變，均符合各組樣本間的生理差異。

此外，研究還發(fā)現(xiàn)冷凍處理對(duì)于轉(zhuǎn)錄組測(cè)序樣本的影響極為顯著，其嚴(yán)重性遠(yuǎn)超本研究中所討論的其他變量。在研究所設(shè)計(jì)的兩組經(jīng)過(guò)冷凍處理的樣本中均出現(xiàn)了極為明顯的降解，表現(xiàn)為檢出的基因數(shù)量較未經(jīng)冷凍組減少了均超過(guò)10 000 條，占全部Unigene的總數(shù)超過(guò)20%。而通過(guò)對(duì)這部分未檢出的基因(11 358條)進(jìn)行功能富集也證實(shí)了廣泛降解的發(fā)生，即這部分基因功能參與了包括核酸結(jié)合及生物合成等多種基礎(chǔ)物質(zhì)代謝機(jī)制。之后，基于冷凍20 min與冷凍20 h的樣本基因檢出總數(shù)相近這一結(jié)果，研究認(rèn)為冷凍處理導(dǎo)致的樣本降解主要發(fā)生在處理前期，且很可能在20 min內(nèi)，而隨著冷凍時(shí)間的增加，降解過(guò)程逐漸變得緩慢。

綜上所述，本研究通過(guò)對(duì)中華按蚊非吸血活體、吸血活體、吸血活體短期冷凍、吸血死體長(zhǎng)時(shí)間冷凍下的轉(zhuǎn)錄組測(cè)序結(jié)果差異進(jìn)行分析與比較，發(fā)現(xiàn)吸血活體為最優(yōu)樣本。各變量間吸血可能引入的污染對(duì)樣本的影響最小，樣本捕獲導(dǎo)致的死亡與冷凍均會(huì)導(dǎo)致樣本降解，其中冷凍時(shí)間對(duì)樣本質(zhì)量的影響極大，可在短時(shí)間內(nèi)造成樣本嚴(yán)重降解。