樣本采集與處理條件對中華按蚊轉(zhuǎn)錄組測序結(jié)果的影響*

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

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

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

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

1 材料與方法

1.1 樣本采集

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

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

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

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

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

測序原始數(shù)據(jù)使用Trimmomatic v0.33軟件(Bolgeretal., 2014)過濾包含接頭、未知堿基N含量大于5%等低質(zhì)量序列。所有過濾后的清潔標(biāo)簽(Clean reads)用于之后輸入Trinity v2.11.0軟件(Grabherretal., 2011)合并組裝，然后使用Tgicl v2.1軟件(Perteaetal., 2003)對轉(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)行評估，評估采用昆蟲(insecta)共有SCO基因數(shù)據(jù)庫。使用TransDecoder v5.3.0軟件(Grabherretal., 2011)識別Unigene中的編碼區(qū)域(Coding domain sequence, CDS)，識別結(jié)合Blast 2.2.28軟件(Camachoetal., 2009)比對SwissProt數(shù)據(jù)庫(McMillanetal., 2008)和Hmmerscan軟件(Chaietal., 2019)搜索Pfam 數(shù)據(jù)庫中(Finnetal., 2014)蛋白同源性的結(jié)果，以獲得更完整編碼區(qū)域。

功能注釋通過將組裝得到的Unigene在常用功能數(shù)據(jù)庫中依據(jù)序列相似性比對獲得KEGG(Kanehisaetal., 2017)、GO(The gene ontology resource)、NR、NT(Bensonetal., 2018)、SwissProt(McMillanetal., 2008)、Pfam(Finnetal., 2014)和KOG(Mudado Mdeetal., 2006)，此過程中采用的軟件分別為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)將各個樣品與Unigene進(jìn)行比對，計算基因表達(dá)水平。R中deseq軟件包(Andersetal., 2010)被用于計算不同樣本組(重復(fù)樣品歸為一組)的差異表達(dá)量分析，使用每百萬堿基檢出片段數(shù)量(Fragments per Kilobase per million，F(xiàn)PKM)值作為衡量表達(dá)量差異的參數(shù)。差異倍數(shù)達(dá)到兩倍以上并且二輪校準(zhǔn)后的P-value(Q-value)≤ 0.001 被用來篩選顯著差異表達(dá)基因。差異基因的GO、KEGG富集分析使用R軟件中的phyper軟件包(Berlivetetal., 2007)，Q-value≤ 0.05的功能視為顯著富集。結(jié)果的可視化通過R軟件中的GG-plot2軟件包(Maag, 2018)完成。

2 結(jié)果

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

轉(zhuǎn)錄組組裝結(jié)果所得的Unigene基因數(shù)據(jù)集及其相關(guān)功能注釋信息是之后比較分析的基礎(chǔ)。本研究所選擇的12個中華按蚊樣本共測得75.58 Gb 數(shù)據(jù)，平均每個樣本測序6.30 Gb。最終去冗余從頭組裝結(jié)果包含48 775個Unigene。序列總大小為 68.1 Mb，平均長度為1 395 bp、N50值為2 349 bp，整體GC含量為47.75 %。進(jìn)一步從這些Unigene中檢測出28 265個編碼序列CDS，同時預(yù)測出5 263個編碼轉(zhuǎn)錄因子。綜合統(tǒng)計Unigene比對到不同基因功能數(shù)據(jù)庫的結(jié)果，顯示共有35 567(72.92%)個基因能夠獲得功能注釋信息。其中各數(shù)據(jù)庫分別注釋了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%)個基因。

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

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

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

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

研究首先比較了各組樣本相對于整體的組裝完整性。結(jié)果顯示，整體Unigene的SCO完整性極高，達(dá)到了98.7%。各組間比較結(jié)果顯示全部樣本的完整性均低于整體，即各組均有部分基因未能檢出。之后，組間比較也發(fā)現(xiàn)明顯差異，其中M1組與M4組的SCO完整性低于M2組與M3組，而M1組2號重復(fù)樣本與M4組3號重復(fù)樣本的完整SCO基因比例甚至低于60%，存在大量片段及未能檢出的缺失的SCO基因(圖3)。該結(jié)果說明至少從SCO抽樣水平，各組檢出的基因數(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個基因。而M1與M2組檢出基因數(shù)量更多(分別為45 887 與44 478個)，二者共有42 584個基因，占全部的87.31%(圖4)。進(jìn)一步，為了揭示M1與M2組所多檢出的基因的功能分布，本研究對在M3與M4組中均未能檢出而M1與M2組均檢出的11 358條基因的功能進(jìn)行了富集。結(jié)果顯示，GO與KEGG富集的通路均與基礎(chǔ)代謝機(jī)制相關(guān)(圖5～6)。如GO中最顯著富集基因為發(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á)水平

對于所有組均檢出的基因，研究比較了其表達(dá)量以分析不同變量在解析表達(dá)水平層面上的影響。結(jié)果顯示，任一兩組間均存在大量差異基因，發(fā)生顯著上調(diào)與下調(diào)的基因總數(shù)均超過10 000 條，超過全部共有基因(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)計圖

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

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

3 討論

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

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

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

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

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

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