顏氏大疣蛛毒腺轉錄組學研究

顧小亮,殷文豪,酉勇明,吳亞應,丁麗君,王 穎,巫秀美,張成桂,楊自忠

(大理大學昆蟲生物醫藥研究院昆蟲生物醫藥研發重點實驗室,藥用特種昆蟲開發國家地方聯合工程研究中心,大理大學藥用昆蟲及蛛形類生物資源數字化開發創新團隊,云南大理 671000)

轉錄組學(Transcriptomics)研究是功能基因組學研究的一個重要內容,它能從整體水平上研究細胞中基因轉錄調控的規律及其轉錄情況,是揭示復雜生物學通路和性狀調控網絡分子機制的學科(崔凱等, 2019)。基于高通量測序技術的轉錄組學主要原理是對真核生物在某個時期內特定組織或細胞轉錄出來的RNA(包括mRNA和非編碼RNA)進行測序,獲得幾乎所有轉錄本信息用于下游文庫構建、上機測序及分析轉錄組信息。具有數據量大、準確性高、效率高和低運行成本等突出優勢,已在各種生物物種轉錄組的研究中有著較為廣泛的應用(李媛媛等, 2021)。

作為一類常見動物,蜘蛛在生物防治、醫藥和材料學等領域有著很好的應用前景,但目前對蜘蛛的轉錄組學研究相對較少,已報道轉錄組學研究的物種包括家福塞勒捕鳥蛛Selenocosmiajiafu(卓城潔等, 2018)、擬環紋豹蛛Pardosapseudoannulata(劉靖等, 2021)、大腹園蛛Araneusventricosus(王康康等, 2020)、擬水狼蛛Piratasubpiraticus(卓俊哲等, 2019)、巴西流浪蜘蛛Phoneutrianigriventer、溝紋硬皮地蛛Calommatasignata等。動物毒素是藥物開發和臨床應用等多領域重要的資源,目前已用于對芋螺、水母、蛇和蝎類等的毒素研究(曹園等, 2018)。

顏氏大疣蛛Macrotheleyani隸屬于蛛形綱Arachinida蜘蛛目Araneae后紡亞目Opisthothelae原蛛下目Mygalomorphae大疣蛛科Macrothelidae,結漏斗狀網、多年生,為夜行性動物,白天隱蔽于土縫或石縫中,晚上爬出洞口捕食落在其網上的昆蟲或動物。主要分布于我國云南省的西北部和南部地區,該物種(Xuetal., 2002)于 2002年首次在云南省高黎貢山地區發現。目前全球已記錄描述的大疣蛛科蜘蛛共計1屬55種,其中分布于我國的大疣蛛有二十余種,占了世界大疣蛛科蜘蛛物種多樣性的一半左右(陳會明等, 2019)。大疣蛛個體相對較大、毒性較強,是農林生態系統中害蟲的重要天敵,也是進行新藥挖掘和毒液開發的重要對象,但其基因信息尚不完善,有待探索。本團隊前期研究發現顏氏大疣蛛毒液中含有豐富的多肽和蛋白類活性物質,藥用價值較高(王夢如等, 2019),具有很好的開發和應用前景。但毒液作為一種寶貴的天然產物,它的來源單一性、性質不穩定性以及成分復雜性等問題一直是影響其開發應用的關鍵因素。隨著近些年發展基因組學、蛋白質、多肽組學和轉錄組學技術越來越成熟,分子生物學方法廣泛應用于研究生物代謝機理規律或細胞生理活動等(賈新平等, 2014)。本研究對顏氏大疣蛛的毒腺轉錄組進行了測序生物信息學分析,開展基因表達譜和功能預測,并對雌雄個體之間的差異表達基因進行篩選,從轉錄組層面了解顏氏大疣蛛毒液作用的分子機制及毒液構成的分子基礎,也為其它有毒節肢動物或昆蟲相關研究提供參考。

1 材料與方法

1.1 樣品準備

顏氏大疣蛛雌雄成體采自于云南高黎貢山福貢縣,采集回實驗室后用透明且打孔的塑料盒于通風的室內進行飼養,飼養條件控制為:溫度25℃,相對濕度60%±5%,避光,每周定時以黃粉蟲Tenebriomolitor喂食并提供水源。選擇個體大小相近的健康雌雄成體各3頭作為實驗樣品,雌雄樣品組分別編號為DYZC和DYZX。所有樣品均為人工取毒后在無菌無酶環境下解剖出毒腺并立即放入液氮凍存。

1.2 方法

1.2.1總RNA提取

將液氮冷凍保存的毒腺組織樣品送至元莘生物醫藥(上海)科技有限公司進行總RNA提取。操作步驟包括:將毒腺組織置于研缽中加適量液氮研磨粉碎后采用Total RNA Extractor試劑盒(貨號:B511311-0500,生工生物工程(上海)股份有限公司)提取毒腺組織樣品總RNA,利用Nanodrop 2000(賽默飛世爾科技,美國)對所提取RNA的濃度和純度進行檢測,單次建庫要求RNA總量不少于5 μg,提取濃度不少于200 ng/μL,OD260/OD280值檢測結果介于1.8～2.2之間;采用2.0%瓊脂糖凝膠電泳回收目的條帶對提取RNA的完整性進行檢測,Agilent 2100(安捷倫科技,美國)測定RIN值。將滿足測序要求的RNA樣品采用雙末端測序(Pair-End)法進行高通量測序。

1.2.2cDNA文庫構建與測序

真核生物mRNA3′末端具有ployA尾的結構,樣品經檢測合格后,先采用帶有Oligo(dT)的磁珠與ploy A進行A-T堿基配對結合,可以富集分離出顏氏大疣蛛毒腺RNA中的mRNA,然后加入fragmentation buffer將mRNA打斷為200 bp左右的短片段序列,以mRNA作為模板,用隨機引物在逆轉錄酶的作用下合成單鏈cDNA(互補DNA),隨后進行雙鏈cDNA合成,再用AMPure XP beads對雙鏈cDNA進行純化,對純化的雙鏈進行末端修飾、片段大小的選擇,最后進行PCR富集,得到顏氏大疣蛛最終的cDNA文庫。使用Agilent 2100對文庫的有效濃度進行準確定量,以保證文庫的質量。庫檢合格后采用Illumina Hiseq測序平臺測序。

1.2.3轉錄組拼接和功能注釋

Illumina Hiseq-測序得到原始圖像數據經堿基識別得到的原始序列數據(raw reads)中包含有低質量序列片段、帶接頭序列、長度過短序列及N率較高序列,為了提高信息分析的最終質量,必須對原始序列進行過濾,得到高質量的clean reads,后續分析都基于此高質量序列數據。具體處理步驟包括:(1)去除低質量reads(測序質量值小于Q 20的堿基數占整個序列長度50%以上的reads);(2)去除N(即無法確定堿基信息)的比例大于10%的序列;(3)去除reads中的adapter序列;(4)剪切前去除5'端含有非AGCT的堿基;(5)除去adapter及質量修剪后的長度小于25 bp的一些小片段。獲得高質量RNA-seq測序數據后,使用Trinity(v 2.11.0)軟件將所有測序讀斷通過從頭組裝生成重疊群(contig)和單一序列(singleton),分別與數據庫(NR、String、Pfam、KOG/COG、GO、Swissprot、KEGG)中的蛋白質參考序列進行比對,取閾值e≤10-10,通過序列相似性對其功能進行注釋。序列相似性比較使用BLAST算法,再用blast2go (http://www. blast2go.com/b2ghome)軟件對所有的unigenes進行GO功能分類統計。然后對unigenes分別進行KOG功能分類和KEGG代謝途徑的分析。

1.2.4差異表達基因的篩選及分析

為了獲得雌性和雄性顏氏大疣蛛表達量不同的基因,使用edgeR v3.24 (http://www.bioconductor. org/packages/release/bioc/html/edgeR.html) R包分析DYZC和DYZX樣品組之間所有基因的差異表達。在進行差異表達檢測過程中,由于轉錄組測序分析會出現假陽性,需要對原有假設檢驗得到的顯著性P值(P-value)采用多重檢驗校正,并最終以log2 (Fold Change)>1,且FDR<0.05指標作為差異基因篩選標準。采用go_naot_exp軟件(v 1.4.0,元莘自定義腳本)和kegg_naot_exp軟件(v 1.4.0,元莘自定義腳本)基于超幾何分布原理對差異基因集進行顯著分析,并將所得差異基因進行GO功能注釋分析和KEGG通路注釋分析,采用go_enrichment(v2.1.1)和kegg_enrichment (v 2.1.1)軟件使用Fisher精確檢驗進行計算,對差異比較組合的所有差異基因集、上調基因集、下調基因集進行富集分析。

2 結果與分析

2.1 轉錄組測序質控與組裝

對獲得的原始數據進行統計,共獲得數據量37.8 G,共計443 073條transcripts,總GC含量平均值為38.17%,數據中96.27%序列的質量參數Q≥20%,90.79%序列的質量參數Q≥30% (NCBI SRA數據庫登錄號:PRJNA761501)。經過分析,最終得到301 024條unigenes,總長度170 512 372 bp,最短201 bp,最長36 105 bp,平均長度566 bp,GC含量平均值為38.17%,N50長度為705;說明測序效果較好,可以進行后續分析(圖1)。

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

2.2 Unigenes功能注釋分類

在公共數據庫中對組裝得到的77 275條unigenes序列進行注釋(表1),不同數據庫中注釋成功的unigenes基因數各不相同。在NR數據庫當中,有31 409條unigenes能夠尋覓到近似序列,約占unigenes總數的10.43%;在Swissprot數據庫中有11 817條unigenes獲得注釋,占3.93%;在GO數據庫當中,有11 480條unigenes獲取了注釋,約為總數的3.81%;在String數據庫中有6 549條 unigenes獲得注釋,占2.18%;在Pfam數據庫中有8 216條unigenes獲得注釋,占2.73%;在KEGG數據庫中獲得注釋的unigenes數量為7 804條,占總體數的2.59%。

表1 顏氏大疣蛛毒腺轉錄組功能注釋Table 1 Functional annotation of transcriptome in Macrothele yani

2.2.1顏氏大疣蛛轉錄組NR數據庫注釋

通過與NCBI-NR數據庫比對,獲得顏氏大疣蛛共31 409條unigenes得到注釋(圖2)。其中,相似性最大的物種為穹蛛Stegodyphusmimosarum,相似序列占總序列數的23.6%,其次是溫室希蛛Parasteatodatepidariorum占總序列數的10%,雕紋似刺尾蝎Centruroidessculpturatus占總序列數的6.7%,白蟻Cryptotermessecundus占總序列數的5.8%等,注釋到的其他物種和未知物種序列占31.9%,其中也包括章魚Octopusbimaculoides、廈門文昌魚Branchiostomabelcheri、德國小蠊Blattellagermanica、褐飛虱Nilaparvatalugens等。

圖2 NR數據庫物種注釋統計圖Fig.2 Statistical chart of species annotation in NR database

2.2.2Unigenes的KOG分類注釋

在KOG數據庫中比對顏氏大疣蛛轉錄組unigenes,結果顯示共8 799條unigenes序列獲得注釋信息,根據其功能大致可分為25類(圖3)。其中,信號傳導機制(Signal transduction mechanisms)、一般功能預測基因(General function prediction only)、翻譯后修飾、蛋白質折疊和分子伴侶(Posttranslational modification, Protein turnover, Chaperones)獲得注釋最多,分別為1 374條、998條、736條,而細胞活性(Cell motility)獲得注釋最少,僅12條。

圖3 顏氏大疣蛛KOG功能統計圖Fig.3 Statistical diagram of KOG function of Macrothele yani

2.2.3Unigenes的KEGG分類注釋

使用KEGG數據庫對顏氏大疣蛛的unigenes可能參與或涉及的代謝途徑進行了統計分析,在組裝得到的轉錄本中,有7 804條得到注釋。這些通路信息分可為5大類,分別為代謝(A, Metabolism)、遺傳信息處理(B, Genetic information processing)、環境信息處理(C, Environmental information processing)、細胞過程(D, Cellular process)、有機系統(E, Organismal systems)。這5大類可以進一步分為45個亞類367個功能通路(圖4)。其中氨基酸代謝(Amino acid metabolism)獲得注釋信息最多,有2 666條。參與較多的幾條信號轉導通路分別是代謝(Metabolism)通路中的氨基酸代謝(Amino acid metabolism)、全局和概覽圖(Global and overview maps),遺傳信息處理(Genetic information processing)中的翻譯(Translation),細胞過程(Cellular process)中的信號轉導(Signal transduction),均大于1 000條。

圖4 顏氏大疣蛛KEGG通路分類注釋結果Fig.4 Statistical diagram of KEGG Pathway of Macrothele yani

2.2.4Unigenes的GO分類注釋

按照GO功能分類方式,將注釋到的11 480條功能注釋序列分為分子功能(Molecular function)、細胞組分(Cellular process)和生物過程(Biological process)三大功能類別(圖5)。所有的這三大類可細分為62個二級分類,并對unigenes在62個二級分類中的分布情況進行統計分析發現:“生物學過程”包含26個不同的亞類,它也是3大類中最多的亞類,其中注釋條目較多的是細胞進程(Cellular process)和新陳代謝進程(Metabolic process),分別為8 372條和6 668條;“細胞成分:組中共分為19個亞類,序列條目較多的是細胞成分(Cell part)和細胞(Cell),分別為9 456條和9 490條;“分子功能”組中包含17個亞類,其中結合活性(Binding)和催化活性(Catalytic activity)的序列數占比最多,分別為8 172條和5 033條。注釋的序列超過5 000的功能條目共有7個,分別是代謝過程(Metabolic process)、細胞器(Organelle)、細胞進程(Cellular process)、細胞(Cell)、細胞組分(Cell part)、催化活性(Catalytic activity)和結合活性(Binding);這一結果從宏觀上展示了顏氏大疣蛛表達基因的功能分布特征,顯示了顏氏大疣蛛毒腺分泌毒液過程中基因表達譜的總體情況。

圖5 顏氏大疣蛛毒腺轉錄組unigenes的GO功能注釋結果Fig.5 Statistical diagram of GO function of Macrothele yani transcriptomic unigenes

2.3 顏氏大疣蛛雌雄基因差異表達分析

分別對DYZC組和DYZX組的轉錄組測序結果進行顯著差異表達基因分析,發現DYZF組和DYZM組有共318個較為顯著的差異表達基因,42 088個非顯著性差異基因,其中雌性毒腺組與雄性組比較表現出205個上調基因,113個下調基因,通過plot_volcano_exp(v 1.1.0)軟件制作火山圖,將非顯著性差異基因和差異表達基因可視化,推斷雌雄毒腺差異基因的整體分布情況(圖6)。從差異表達基因數據來看,雌雄個體顯著性差異表達基因數量與非顯著性差異基因數量相當,表明顏氏大疣蛛雌雄個體毒液成分或毒性大小等相差不大。

圖6 顏氏大疣蛛雌性與雄性間差異表達基因的比較Fig.6 Comparison of differential expressed genes between male and female Macrothele yani注:圖中每個點代表一個特定的基因或轉錄本,紅色圓點表示顯著上調的基因,綠色圓點表示顯著下調的基因,黑色圓點為非顯著差異基因。在左邊的點為表達差異下調的基因,右邊的點為表達差異上調的基因,越靠左邊和上邊的點表達差異越顯著。Note: Each dot in the figure represented a specific gene or transcript. Red dots represented genes that were significantly up-regulated, green dots represented genes that were significantly down-regulated, and black dots represented genes that were nonsignificant differential genes. The dots on the left were differentially down-regulated genes, and the dots on the right were differentially up-regulated genes. The more the dots on the left and top were, the more significant the differences in expression were.

2.4 差異表達基因GO和KEGG顯著性富集分析

對篩選的差異表達基因進行GO功能注釋分類顯示:共有440條unigenes被GO注釋,并被顯著富集到分子功能(Molecular function)、細胞進程(Cellular process)和生物過程(Biological process)三大類的38個亞類中(圖7)。其中生物調節(Biological regulation)、細胞成分組織或合成(Cellular component organization or biogenesis)、細胞過程(Cellular process)、代謝過程(Metabolic process)、應激反應(Response to stimulus)是生物過程中富集差異基因數量較多的亞類;細胞(Cell)、細胞成分(Cell part)、細胞器(Organelle)、細胞器成分(Organelle part)是細胞組分中富集差異基因數量較多的亞類;結合(Binding)和催化活性(Catalytic activity)是分子功能中富集差異基因數量最多的亞類。其中富集差異基因數目最多的是分子功能(Molecular function)大類中表達上調的結合(Binding)亞類的30條unigenes。表達下調最多的是細胞(Cell)和細胞成分(Cell part)均為12條unigenes。

圖7 顏氏大疣蛛差異表達基因GO功能注釋分類結果Fig.7 Differentially expressed genes of Macrothele yani in GO functional classification注:圖中下方橫坐標表示注釋到的GO term亞類的基因個數,上方橫坐標表示注釋到的GO term的序列數占所有注釋基因序列總數的比例(功能基因亞類和GO term大類是多對多的關系,即一個亞類基因序列可包含多個GO term的功能注釋,某一個GO term大類也會對應到多個不同的基因序列,并不是一對一的關系);縱坐標表示GO注釋的每一個詳細分類,包括左側三個大類和中間的38個亞類;紅色條形柱代表顯著上調基因,藍色條形柱代表顯著下調基因。Note: In the figure, the lower horizontal coordinate represented the number of genes in the annotated GO term subclass, while the upper horizontal coordinate represented the proportion of the annotated GO term sequence number to the total number of all annotated gene sequences (the relationship between functional gene subclass and GO term class was many-to-many, that was, one subclass gene sequence could contain multiple functional annotations of GO term. A GO term category would also correspond to multiple different gene sequences, rather than a one-to-one relationship). The vertical axis represented each detailed classification of the GO annotation, including the three categories on the left and 38 subcategories in the middle. The red bars represented significantly up-regulated genes and the blue bars represented significantly down-regulated genes.

2.5 功能蛋白鑒定

毒腺cDNA文庫構建后,通過分子生物學手段結合基因注釋信息將得到的蛋白多肽序列與Uniprot數據庫中的功能蛋白進行BLAST比對,進一步篩選毒腺中轉錄表達的毒素蛋白分子序列信息即搜索功能蛋白序列,對其進行功能預測,共計獲得16類68條表達毒素蛋白的相關基因序列(圖8)。其中轉錄表達的毒素蛋白包括乙酰膽堿脂酶(Acetylcholinesterase)、蝦紅素金屬蛋白酶(Astacin metalloprotease)富含半胱氨酸的分泌蛋白家族(CRISP(Allergen)、皮膚壞死毒素(Dermonecrotic toxin)、透明質酸酶(Hyaluronidase)、胰島素生長因子相關蛋白(Insult growth factor related proteia)、Kunitz型抑制劑(Kunitz-type)、腦啡肽酶(Neprilysin)、腫瘤蛋白同源物(controlled tumor protein homolog)等,它們在雌雄顏氏大疣蛛毒腺轉錄組編碼的蛋白序列中中均有發現。

圖8 顏氏大疣蛛毒腺轉錄本編碼的毒素蛋白Fig.8 Toxin proteins encoded by transcript of the venom glands of Macrothele yani

3 結論與討論

全世界目前已知蜘蛛物種近五萬種,是地球上生物多樣性最豐富的動物類群之一,以昆蟲或其他節肢動物為食。蜘蛛毒液的研究是繼蛇毒、蝎毒研究后的又一大熱點,但目前僅有少數物種有毒素研究,且研究深度和廣度都有待提升(丁麗君等, 2021)。現已證明蜘蛛毒液中的成分對神經遞質、離子通道及一些受體有明顯的激動或抑制作用,國內外不少學者對蜘蛛毒液的成分、體內外活性和作用機制進行研究,并取得了一些成果。例如:對棒毛絡新婦Trichonephilaclavata毒液的分離純化及各組分對昆蟲的生物活性研究(Jinetal., 2017);對塔蘭圖拉狼蛛Lycosatarantula毒液的蛋白質組學和轉錄組學研究以來分析蜘蛛毒液中的毒素蛋白(Kouaetal., 2020)等;但毒液本身因不穩定性、復雜性等生物學性質影響了其廣泛應用,因此對其進行基因組學研究,從生物信息學方面較易獲取顏氏大疣蛛毒液中有價值的數據,利用毒腺進行cDNA文庫的構建結合基因工程表達的方法可以有效解決毒液的來源問題(卓城潔等, 2018)。

由于缺乏蜘蛛的基因組、蛋白質序列等,本研究僅在數據庫中匹配到一部分已知的基因序列,其中從毒腺轉錄組數據中比對到68條毒素蛋白相關序列,包括大多數的酶和蛋白質等。其中,酶是一種催化生物轉化的蛋白質或RNA分子,一直是動物毒液的重要組成部分。如乙酰膽堿酯酶(Acetylcholinesterase)與受體結合并干擾乙酰膽堿的分解,導致乙酰膽堿在神經突觸中沉積并導致神經傳遞中斷(Thapaetal., 2017);透明質酸酶(Hyaluronidase)可以調節毒液擴散機制或防御機制,通過分解間質空間中的多糖來促進物種在細胞外基質中的擴散(Helm Iietal., 2019)。皮膚壞死毒素(Dermonecrotic toxin)屬于磷脂酶-D毒素家族一員,是蜘蛛毒素中研究最多的一類,這類毒素可誘發炎癥反應、溶血、血小板減少、皮膚壞死和腎功能衰竭等(Chaves-Moreiraetal., 2017)。在同為穴居型的溝紋硬皮地蛛Calommatasignata中也鑒定到5條kunitz型毒素和2條透明質酸酶毒素序列,其它家族的毒素未鑒定到(Hanetal., 2021);在Pamphobeteusverdolaga蜘蛛毒腺轉錄組中同樣鑒定到了kunitz型毒素和透明質酸酶(Estrada-Gómezetal., 2021),但磷脂酶等在顏氏大疣蛛中未鑒定到;在海南捕鳥蛛Haplopelmahainanum毒腺轉錄本數據中鑒定到了201種潛在的毒素蛋白序列,包括knottins、Kunitz型毒素、酶和其他蛋白質(Chengetal., 2016)。在巴西流浪蜘蛛Phoneutrianigriventer中除了富含半胱氨酸的肽毒素外,還鑒定了其他假定的毒液成分,如蛋白酶抑制劑、防御素和絲氨酸蛋白酶(Paivaetal., 2019)。一些毒素蛋白或酶在螞蟻、蝎子等節肢動物(Santosetal., 2017; Díazetal., 2019)毒液中也有報道。

物種的雌雄分化是進化的產物,大多數種類在形態上都表現出性別二態性,即同一物種除生殖器結構不同外,還表現出體型大小或顏色等的一些明顯差異(王孟卿等, 2005),但在轉錄組學方面還未充分了解其差異。為獲得雌雄顏氏大疣蛛毒腺的轉錄組信息,我們應用高通量測序技術對雌雄顏氏大疣蛛毒腺進行了轉錄組基因測序,通過序列的預處理、拼接、轉錄本功能注釋等步驟,結合生物信息學分析方法,對顏氏大疣蛛毒腺unigene與NR、String、Swissprot、Pfam、GO、KEGG數據庫進行比對(e≤10-10),以強大的生物信息學平臺作支撐,根據“基因結構相似,功能同源”的原理,對基因的功能進行注釋。分別獲得31 409、6 549、11 817、8 216、11 480、7 804 條unigenes注釋,共計301 024個結果,并分析了顏氏大疣蛛毒腺的轉錄組特征。表明在對顏氏大疣蛛基因組不清楚的情況下,高通量測序技術是批量發現物種功能基因的有效且簡易手段(蔣月麗等, 2020),更多具有研究價值的潛在信息將被展示出來。

本研究中利用高通量測序技術初步建立了顏氏大疣蛛毒腺轉錄組數據庫,獲得了大量有關毒腺轉錄組學的序列信息,對表達的基因進行序列拼接及組裝,與公共數據庫的功能蛋白序列比對,進行相應的功能注釋及代謝途徑分析,并比較了雌雄個體間的差異表達基因。該轉錄組數據為今后更深入研究顏氏大疣蛛毒液的功能基因組和基因克隆等提供了大量的信息資源,也可為其今后的分子標記和功能基因挖掘、農業新型殺蟲劑及新藥研究開發提供有價值的信息,為其它有毒節肢動物的基因組提供參考。