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

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

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

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

1 材料與方法

1.1 材料

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

1.2 轉(zhuǎn)錄組測序

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

1.3 轉(zhuǎn)錄組組裝

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

1.4 功能注釋

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

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

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

2 結(jié)果

2.1 測序結(jié)果

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

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

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

通過Trinity拼接組裝共獲得848 451條轉(zhuǎn)錄本（表2），拼接的總長度為1 165 234 924 bp，平均長度為1 373.37 bp，N50為4 474 bp。通過tr2aacds去除冗余片段后，共獲得308 343條轉(zhuǎn)錄本，總長度為406 642 215 bp，平均長度為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個數(shù)據(jù)庫進(jìn)行同源性比對和功能注釋，共獲得有注釋信息的83 379條轉(zhuǎn)錄本：Swiss-port數(shù)據(jù)庫75 138條、KEGG數(shù)據(jù)庫59 118條、GO數(shù) 據(jù) 庫66 742條、Pfam數(shù) 據(jù) 庫39 738條、EggNOG數(shù)據(jù)庫55 604條。同時在5個數(shù)據(jù)庫中有注釋信息的轉(zhuǎn)錄本共31 019條（圖1）。

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

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

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

將轉(zhuǎn)錄本分為生物學(xué)過程、細(xì)胞組分和分子功能，其中，細(xì)胞組分中注釋到的轉(zhuǎn)錄本最多，共60 651條，占90.8%；生物學(xué)過程中注釋到的轉(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個組織樣本可被清楚分為5個聚類簇，其中，心臟、肌肉和胃有一定的重疊，彼此顯示出較近的聚類關(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個差異表達(dá)基因（表3）：肺臟與肝臟的差異表達(dá)基因數(shù)量最多（12 442個），肌肉與心臟的最少（4 005個）。

表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）及通過調(diào)節(jié)游離鈣離子的釋放來調(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)，如：單核細(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)的信號通路（GO：0070098）、細(xì)胞因子介導(dǎo)的信號通路（GO：0019221）、整合素介導(dǎo)的細(xì)胞粘附負(fù)調(diào)控（GO：0033629）、巨噬細(xì)胞趨化性（GO：0048246）、中性粒細(xì)胞趨化性（GO：0030593）、中性粒細(xì)胞活化（GO：0042119）、脂多糖結(jié)合（GO：0001530）、細(xì)胞對脂多糖的反應(yīng)（GO：0071222）和細(xì)胞對白細(xì)胞介素-1的反應(yīng)（GO：0071347）等。此外，還富集到與氧氣運輸、血管生成、造血相關(guān)的通路。肺臟富集到了21條KEGG通路，其中大部分也是免疫與疾病相關(guān)的通路，包括Toll樣受體信號通路（ko04620）、趨化 因 子 信 號 通 路（ko04062）、IL-17信 號 通 路（ko04657）、JAK-STAT信號通路（ko04630）等。差異表達(dá)基因IL8、NFKB1、IL3RB等被富集在KEGG通路中。

肝臟一共富集到了298條GO條目和44條KEGG通路。GO富集結(jié)果分為以下三大類，一是維生素A代謝相關(guān)通路，如對維生素A的反應(yīng)（GO：0033189）、視 黃 醇 代 謝 過 程（GO：0042572）和維甲酸代謝過程（GO：0042573）等；二是脂類代謝的相關(guān)通路，如脂肪酸結(jié)合（GO：0005504）、高 密 度 脂 蛋 白 顆 粒（GO：0034364）、膽固醇代謝過程（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對外源物質(zhì)的代謝（ko00980）、藥物代謝-細(xì)胞色素P450（ko00982）。

3 討論

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

長期生活在高海拔地區(qū)，體內(nèi)血液的分布也會發(fā)生變化，心血管系統(tǒng)的血流分布會相對增多，而腎臟的血流會相對減少，刺激腎素-血管緊張素-醛固酮系統(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ān)通路以及ANP協(xié)同調(diào)節(jié)血壓，減少低氧刺激導(dǎo)致的功能紊亂的影響，這可能也是大噪鹛在高海拔生存時的機體優(yōu)化策略之一。

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

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