大噪鹛多組織轉(zhuǎn)錄組測(cè)序組裝和基因注釋

李仁明秀，劉乙，蘭月，耿煬，楊楠，岳碧松*

（1.四川大學(xué)生命科學(xué)學(xué)院,生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,成都 610065;2.四川省生態(tài)環(huán)境科學(xué)研究院,成都 610041;3.西南民族大學(xué)青藏高原研究院,成都 610041）

大噪鹛隸屬于雀形目Passeriformes噪鹛科Leiothrichdae噪鹛屬，為該科體型最大的物種之一（Cheng，1987；Rasmussen&Anderton，2005），主要分布在四川、重慶、湖北、甘肅南部、西藏東南部以及云南西北部，我國(guó)特有種（鄭光美，2017），國(guó)家二級(jí)重點(diǎn)保護(hù)野生動(dòng)物（國(guó)家林業(yè)和草原局，農(nóng)業(yè)農(nóng)村部，2021），其垂直分布范圍跨度較大（海拔2 135～4 115 m），甚至在海拔4 200 m以上地區(qū)也能生存（張強(qiáng)等，2010），棲息環(huán)境主要包括亞高山、高山深林灌叢及其周圍的林緣地帶（del Hoyo，2007）。

隨著高通量測(cè)序技術(shù)的飛速發(fā)展，利用轉(zhuǎn)錄組測(cè)序分析技術(shù)研究基因的表達(dá)調(diào)控的差異越來(lái)越廣泛（Bai，2016；Zhang，2016；Wang，2017；Zhao，2018）。目前通過(guò)高通量測(cè)序分析，對(duì)于高海拔低氧環(huán)境的綠尾紅雉（Qu，2013）、藏豬（Ai，2014）、褐 背 擬 地 鴉（Wang，2015）、藏雞（Cui，2019）等已有相關(guān)適應(yīng)性的研究。大噪鹛分子遺傳信息匱乏，目前僅有關(guān)于甘肅蓮花山大噪鹛的繁殖習(xí)性研究（Wang，2010）。本研究通過(guò)二代測(cè)序技術(shù)對(duì)采自四川省甘孜藏族自治州帕姆嶺和察青松多白唇鹿國(guó)家級(jí)自然保護(hù)區(qū)的8只大噪鹛的39個(gè)組織樣本進(jìn)行轉(zhuǎn)錄組測(cè)序，通過(guò)組裝、注釋，獲得了第一個(gè)組裝質(zhì)量較高的大噪鹛轉(zhuǎn)錄組，為進(jìn)一步研究其遺傳學(xué)和瀕危機(jī)制提供基礎(chǔ)數(shù)據(jù)，通過(guò)不同組織間的差異表達(dá)分析大噪鹛高海拔生存適應(yīng)機(jī)制。

1 材料與方法

1.1 材料

大噪鹛樣品共8只（均為多年鳥類監(jiān)測(cè)過(guò)程中，撞網(wǎng)受傷而無(wú)法救治的個(gè)體）：5只采自四川省甘孜藏族自治州雅江縣八角樓鄉(xiāng)帕姆嶺（101°18′E，30°10′N，海拔4 147 m），編號(hào)分別為PA、PB、PC、PD、PE，取心臟、肝臟、胃、肺臟、腎臟和肌肉6種組織，編號(hào)分別為1～6；3只采自四川省甘孜藏族自治州白玉縣察青松多白唇鹿國(guó)家級(jí)自然保護(hù)區(qū)（99°20′E，30°59′N，海拔3 800 m），編號(hào)分別為CA、CB、CC，取心臟、肺臟、肌肉3種組織，編號(hào)分別為1～3；共計(jì)39個(gè)組織樣本，液氮保存。由于察青松多采集的心臟與肺臟組織樣本未明確標(biāo)記，為避免混淆，組裝過(guò)程采用39個(gè)組織樣本以獲得更完整的轉(zhuǎn)錄本，后續(xù)分析過(guò)程中僅采用帕姆嶺地區(qū)的組織樣本。

1.2 轉(zhuǎn)錄組測(cè)序

使 用Invitrogen的TRIzol試 劑 盒（Promega，USA）進(jìn)行總RNA提取，操作流程按照說(shuō)明書進(jìn)行，提取后交付北京諾禾致源科技股份有限公司建庫(kù)，并在Illumina Hiseq2000高通量測(cè)序平臺(tái)進(jìn)行雙末端測(cè)序，測(cè)序片段長(zhǎng)約200 bp。

1.3 轉(zhuǎn)錄組組裝

在Trinity 2.4.0（Grabherr，2011）中 以把39個(gè)組織合并組裝成共同的無(wú)參轉(zhuǎn)錄本（Geng，2017），只保留片段長(zhǎng)度在300 bp以上的contig序列（Haas，2013）。組裝參數(shù)：Trinitymax_memory 200G-CPU 80-min_contig_length 300。組裝完成后再通過(guò)EvidentialGene的tr2aacds（http：//arthropods.eugenes.org/EvidentialGene/about/EvidentialGene_trassembly_pipe.html）去 除 片 段 化的、低編碼潛力的以及序列高度相似的轉(zhuǎn)錄本。組裝完成后，選取最長(zhǎng)的轉(zhuǎn)錄本作為Unigene。

1.4 功能注釋

使用BLAST（Altschul，2012）將組裝完成的轉(zhuǎn)錄本與Swiss-prot、Pfam、EggNOG數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，E值設(shè)置為1E-5。將Swiss-port數(shù)據(jù)庫(kù)的注釋結(jié)果導(dǎo)入Blast2GO（Conesa，2005）進(jìn)行GO注釋，通過(guò)WEGO 2.0（Ye，2018）進(jìn)行可視化展示。以真核生物為背景基因集，使用KEGG Automatic Annotation Server（KAAS）（Yuki，2007）對(duì)轉(zhuǎn)錄本進(jìn)行KEGG通路分析。

1.5 基因差異表達(dá)分析

通過(guò)Salmon（Patro，2017）計(jì)算轉(zhuǎn)錄本表達(dá)量得到每千個(gè)堿基的轉(zhuǎn)錄每百萬(wàn)映射讀取的轉(zhuǎn)錄本數(shù)（transcripts per kilobase of exon model per million mapped reads，TPM），對(duì)帕姆嶺的30個(gè)樣本進(jìn)行重復(fù)性對(duì)比和主成分分析，通過(guò)Trinity的EdgeR包（Smyth，2010）進(jìn)行基因差異表達(dá)分析。差異表達(dá)基因的篩選標(biāo)準(zhǔn)為|log（fold change）|≥2且＜0.05。使用ClusterProfile（Yu，2012）對(duì)差異表達(dá)基因進(jìn)行GO功能注釋和KEGG通路富集分析。

2 結(jié)果

2.1 測(cè)序結(jié)果

共得到原始reads 914 244 993條（約274 Gb），過(guò)濾后得到的高質(zhì)量reads共896 246 793條（約269 Gb），有效reads的比例高于96.3%，質(zhì)量在Q20以上的reads超過(guò)95%，質(zhì)量在Q30以上的reads超過(guò)88%（表1）。證明本次大噪鹛轉(zhuǎn)錄組樣品測(cè)序質(zhì)量與過(guò)濾結(jié)果較好，可用于后續(xù)分析。

表1 各樣品測(cè)序結(jié)果匯總Table 1 Summary of sequencing results of each sample

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

通過(guò)Trinity拼接組裝共獲得848 451條轉(zhuǎn)錄本（表2），拼接的總長(zhǎng)度為1 165 234 924 bp，平均長(zhǎng)度為1 373.37 bp，N50為4 474 bp。通過(guò)tr2aacds去除冗余片段后，共獲得308 343條轉(zhuǎn)錄本，總長(zhǎng)度為406 642 215 bp，平均長(zhǎng)度為1 318.80 bp，N50為2 748 bp。

表2 大噪鹛的轉(zhuǎn)錄組組裝Table 2 Transcriptome assembly of Garrulax maximus

2.3 轉(zhuǎn)錄本注釋

將組裝好的轉(zhuǎn)錄本與5個(gè)數(shù)據(jù)庫(kù)進(jìn)行同源性比對(duì)和功能注釋，共獲得有注釋信息的83 379條轉(zhuǎn)錄本：Swiss-port數(shù)據(jù)庫(kù)75 138條、KEGG數(shù)據(jù)庫(kù)59 118條、GO數(shù) 據(jù) 庫(kù)66 742條、Pfam數(shù) 據(jù) 庫(kù)39 738條、EggNOG數(shù)據(jù)庫(kù)55 604條。同時(shí)在5個(gè)數(shù)據(jù)庫(kù)中有注釋信息的轉(zhuǎn)錄本共31 019條（圖1）。

圖1 大噪鹛轉(zhuǎn)錄本功能注釋Fig.1 Function annotation of Garrulax maximus transcripts

按照細(xì)胞過(guò)程和信號(hào)傳遞、信息儲(chǔ)存與處理、新陳代謝以及無(wú)特征基因?qū)⒆⑨尩降霓D(zhuǎn)錄本分為4大類25小類：無(wú)顯著特征的轉(zhuǎn)錄本15 830條（46.9%），其次為信號(hào)傳導(dǎo)機(jī)制（8 373條）、翻譯后修飾轉(zhuǎn)運(yùn)（5 911條）和細(xì)胞內(nèi)運(yùn)輸（3 590條）（圖2）。

圖2 大噪鹛轉(zhuǎn)錄本的EggNOG功能分類Fig.2 EggNOG analysis of Garrulax maximus transcripts

將轉(zhuǎn)錄本分為生物學(xué)過(guò)程、細(xì)胞組分和分子功能，其中，細(xì)胞組分中注釋到的轉(zhuǎn)錄本最多，共60 651條，占90.8%；生物學(xué)過(guò)程中注釋到的轉(zhuǎn)錄本最少，只有1 380條，占2.1%（圖3）。

圖3 大噪鹛轉(zhuǎn)錄本的GO功能分類Fig.3 GO analysis of Garrulax maximus transcripts

注釋到轉(zhuǎn)錄本最多的前5條通路分別為：碳水化合物代謝（ko01200）、核糖體（ko03010）、氨基酸生物合成（ko01230）、嘌 呤 代 謝（ko00230）和 癌 癥 通 路（ko05200），分別注釋到1 345條、1 268條、1 103條、1 045條和1 004條。

2.4 樣本重復(fù)性分析

帕姆嶺樣本中，PC4、PC6和PE6為離群樣本，剔除后構(gòu)建相關(guān)性矩陣，可分為6組，依次為心臟、肌肉、胃、肝臟、腎臟和肺臟，各組的樣本重復(fù)性均較好。

主成分分析結(jié)果顯示，27個(gè)組織樣本可被清楚分為5個(gè)聚類簇，其中，心臟、肌肉和胃有一定的重疊，彼此顯示出較近的聚類關(guān)系；其余3種組織分別聚類在不同方位（圖4）。

圖4 帕姆嶺地區(qū)大噪鹛6種組織樣本主成分分析結(jié)果Fig.4 Principle component analysis of 6 tissue samples from Garrulax maximus in Pamuling area

2.5 差異表達(dá)基因分析

共得到36 059個(gè)差異表達(dá)基因（表3）：肺臟與肝臟的差異表達(dá)基因數(shù)量最多（12 442個(gè)），肌肉與心臟的最少（4 005個(gè)）。

表3 大噪鹛6種組織間差異表達(dá)基因數(shù)Table 3 Number of differentially expressed genes among 6 tissues of Garrulax maximus

心臟的差異表達(dá)基因共富集在228條GO條目上，主要是與心臟形態(tài)結(jié)構(gòu)發(fā)育相關(guān)的條目，如心房肌細(xì)胞發(fā)育（GO：0055014）、心肌組織形態(tài)發(fā)生（GO：0055008）、胚胎心管形態(tài)發(fā)生（GO：0003143）、房間隔形態(tài)發(fā)生（GO：0060413）及心肌纖維發(fā)育（GO：0048739）等；還顯著富集到心臟收縮調(diào)節(jié)的相關(guān)條目，如心肌收縮（GO：0060048）、肌肉收縮（GO：0006936）、肌肉收縮調(diào)節(jié)（GO：0006937）及通過(guò)調(diào)節(jié)游離鈣離子的釋放來(lái)調(diào)節(jié)心肌收縮（GO：0010881）等；此外，還富集到呼吸鏈上能量傳遞的許多通路，如呼吸鏈復(fù)合物Ⅳ（GO：0045277）、氧化磷酸化（GO：0006119）、線粒體呼吸鏈復(fù)合體Ⅰ（GO：0005747）等。心臟一共富集到了18條KEGG通路，除了心臟收縮功能相關(guān)的通路，還富集到了與體液滲透壓調(diào)節(jié)相關(guān)的腎素分泌（ko04924）、醛固酮合成與分泌（ko04925），并且心鈉肽差異表達(dá)基因（natriuretic peptide precursor A，NPPA）富集在心臟的多條通路上。

肺臟一共富集到了43條GO條目，其中有21條富集結(jié)果與免疫反應(yīng)相關(guān)，如：?jiǎn)魏思?xì)胞趨化性（GO：0002548）、抗菌肽介導(dǎo)的抗菌體液免疫應(yīng)答（GO：0061844）、淋 巴 細(xì) 胞 趨 化 性（GO：0048247）、趨化因子活性（GO：0008009）、炎癥反應(yīng)（GO：0006954）、炎 癥 反 應(yīng) 的 正 調(diào) 控（GO：0050729）、趨化因子介導(dǎo)的信號(hào)通路（GO：0070098）、細(xì)胞因子介導(dǎo)的信號(hào)通路（GO：0019221）、整合素介導(dǎo)的細(xì)胞粘附負(fù)調(diào)控（GO：0033629）、巨噬細(xì)胞趨化性（GO：0048246）、中性粒細(xì)胞趨化性（GO：0030593）、中性粒細(xì)胞活化（GO：0042119）、脂多糖結(jié)合（GO：0001530）、細(xì)胞對(duì)脂多糖的反應(yīng)（GO：0071222）和細(xì)胞對(duì)白細(xì)胞介素-1的反應(yīng)（GO：0071347）等。此外，還富集到與氧氣運(yùn)輸、血管生成、造血相關(guān)的通路。肺臟富集到了21條KEGG通路，其中大部分也是免疫與疾病相關(guān)的通路，包括Toll樣受體信號(hào)通路（ko04620）、趨化 因 子 信 號(hào) 通 路（ko04062）、IL-17信 號(hào) 通 路（ko04657）、JAK-STAT信號(hào)通路（ko04630）等。差異表達(dá)基因IL8、NFKB1、IL3RB等被富集在KEGG通路中。

肝臟一共富集到了298條GO條目和44條KEGG通路。GO富集結(jié)果分為以下三大類，一是維生素A代謝相關(guān)通路，如對(duì)維生素A的反應(yīng)（GO：0033189）、視 黃 醇 代 謝 過(guò) 程（GO：0042572）和維甲酸代謝過(guò)程（GO：0042573）等；二是脂類代謝的相關(guān)通路，如脂肪酸結(jié)合（GO：0005504）、高 密 度 脂 蛋 白 顆 粒（GO：0034364）、膽固醇代謝過(guò)程（GO：0008203）等；三是發(fā)揮解毒功能的相關(guān)通路，如有毒物質(zhì)結(jié)合（GO：0015643）、藥物結(jié)合（GO：0008144）、細(xì)胞葡萄糖醛酸化（GO：0052695）、葡萄糖醛酸轉(zhuǎn)移酶活性的負(fù)調(diào)節(jié)（GO：1904224）等。KEGG富集結(jié)果主要也是與物質(zhì)合成代謝相關(guān)的通路，以及解毒功能相關(guān)的細(xì)胞色素P450對(duì)外源物質(zhì)的代謝（ko00980）、藥物代謝-細(xì)胞色素P450（ko00982）。

3 討論

本研究利用Illumina測(cè)序技術(shù)首次對(duì)8只大噪鹛的肝臟、心臟、肌肉、胃、肺臟以及腎臟等6個(gè)組織進(jìn)行了轉(zhuǎn)錄組測(cè)序、組裝以及功能注釋。在差異表達(dá)基因的功能富集結(jié)果中，肺臟組織中富集到大量與免疫相關(guān)的Toll樣受體，Toll樣受體可以識(shí)別革蘭氏陽(yáng)性菌細(xì)胞壁的肽多糖、革蘭氏陰性菌細(xì)胞壁的脂多糖以及病毒的遺傳物質(zhì)雙鏈RNA等相對(duì)保守的結(jié)構(gòu)序列（Aderem & Ulevitch，2000），屬于機(jī)體的固有免疫（Akira & Takeda，2004），可以快速激活細(xì)胞內(nèi)信號(hào)級(jí)聯(lián)反應(yīng)過(guò)程中的一系列催化蛋白及信號(hào)因子，使機(jī)體快速產(chǎn)生促炎性細(xì)胞因子、趨化因子等，從而調(diào)控免疫細(xì)胞的功能并清除病原體。由此說(shuō)明了長(zhǎng)期生活在高海拔地區(qū)的大噪鹛可能通過(guò)增強(qiáng)先天免疫以抵御呼吸過(guò)程中外界病原體的侵襲。

長(zhǎng)期生活在高海拔地區(qū)，體內(nèi)血液的分布也會(huì)發(fā)生變化，心血管系統(tǒng)的血流分布會(huì)相對(duì)增多，而腎臟的血流會(huì)相對(duì)減少，刺激腎素-血管緊張素-醛固酮系統(tǒng)（Hurtado-Arestegui，2018），導(dǎo)致血壓升高。大噪鹛心臟組織中發(fā)現(xiàn)了腎素分泌、醛固酮合成與分泌通路，并發(fā)現(xiàn)NPPA富集在多條通路上。NPPA基因編碼心鈉肽（atrial natriuretic peptide，ANP）的前體，在高血容量以及高血壓情況下，心肌細(xì)胞釋放ANP到血液循環(huán)中，ANP的主要作用是促進(jìn)腎臟利鈉利尿，放松血管平滑肌，從而調(diào)節(jié)血容量和血壓（Song，2015）。此外，ANP還能抑制醛固酮的產(chǎn)生（Maack，1984）。大噪鹛通過(guò)腎素、醛固酮相關(guān)通路以及ANP協(xié)同調(diào)節(jié)血壓，減少低氧刺激導(dǎo)致的功能紊亂的影響，這可能也是大噪鹛在高海拔生存時(shí)的機(jī)體優(yōu)化策略之一。

在肝臟的差異表達(dá)分析中，篩選到了大量上調(diào)的谷胱甘肽轉(zhuǎn)移酶（GST）基因。GST是谷胱甘肽過(guò)氧化物酶家族的成員，主要在清除自由基系列反應(yīng)的前期催化氧自由基與谷胱甘肽結(jié)合（Hayes，2005）。氧自由基不及時(shí)清除會(huì)引起細(xì)胞破壞、蛋白質(zhì)破壞、DNA破壞，最終導(dǎo)致氧化應(yīng)激。長(zhǎng)期的高原生活還可能因?yàn)槁匀毖醵鹑毖跣愿渭?xì)胞損傷，其臨床病理表現(xiàn)為肝細(xì)胞腫脹和脂肪變性（萬(wàn)平新等，2016）。缺氧情況下機(jī)體會(huì)產(chǎn)生更多的活性氧（Dawson，1993）。因此，雖然氧化應(yīng)激反應(yīng)在低氧環(huán)境下難以避免，大噪鹛肝臟組織通過(guò)上調(diào)谷胱甘肽過(guò)氧化物酶家族的相關(guān)基因，可以盡可能降低氧化應(yīng)激反應(yīng)所帶來(lái)的傷害。

綜上所述，本研究基于大噪鹛心臟、肝臟、胃、肺臟、腎臟、肌肉6種組織轉(zhuǎn)錄組表達(dá)量的差異，對(duì)各個(gè)組織的差異表達(dá)基因進(jìn)行了富集分析，結(jié)果發(fā)現(xiàn)心臟、肝臟、肺臟3種組織中有許多組織特異性的差異表達(dá)基因，并且表達(dá)量差異顯著，表明這些組織在大噪鹛低氧習(xí)服的過(guò)程中發(fā)揮著重要作用。但受限于保護(hù)物種，本次研究的樣品數(shù)量較少，后續(xù)還可以進(jìn)行多地非損傷性采樣，為大噪鹛的保護(hù)提供更多的線索與數(shù)據(jù)。