段曉晨，程起群

(1.中國水產(chǎn)科學(xué)研究院東海水產(chǎn)研究所，上海 200090；2.上海海洋大學(xué)水產(chǎn)與生命學(xué)院，上海 201306)

綠鰭馬面鲀(Thamnaconusseptentrionalis)，隸屬于鲀形目(Tetraodontiformes)，單角鲀科(Monacanthidae),馬面鲀屬，又稱馬面魚、剝皮魚等，是我國沿海較為常見的經(jīng)濟魚類[1]。經(jīng)過近40年的高強度開發(fā)，綠鰭馬面鲀資源量迅速衰減，目前已不能充分滿足市場需求[2]。以“養(yǎng)”代“捕”，是緩解野生資源壓力、保障綠鰭馬面鲀市場供給的有效途徑。2011年，綠鰭馬面鲀的人工繁育技術(shù)取得突破[3]，之后開始應(yīng)用于工廠化養(yǎng)殖。

密度過高是水產(chǎn)養(yǎng)殖活動中的重要限制因子[4]，密度脅迫會使種群內(nèi)對空間、食物和氧氣等資源的競爭加劇，基因表達和適應(yīng)能力受到影響[5]。如SUN等[6]的研究表明，擁擠可以顯著誘導(dǎo)魚的免疫反應(yīng)，而長期處于密度脅迫下會損害魚類的免疫能力；QIANG等[7]的研究顯示，高密度條件會改變肌肉脂肪酸的組成，刺激脂質(zhì)代謝酶顯著上調(diào)。LIU等[8]發(fā)現(xiàn)，高密度養(yǎng)殖導(dǎo)致多種消化酶的活性顯著降低，進而影響生物的免疫能力。目前，關(guān)于綠鰭馬面鲀腦組織對密度脅迫的分子響應(yīng)機制尚不明確。

轉(zhuǎn)錄組學(xué)是揭示生物體內(nèi)分子機制的重要手段[9]，應(yīng)用轉(zhuǎn)錄組測序技術(shù)可以對水產(chǎn)動物進行高通量測序，通過測序結(jié)果進一步實現(xiàn)基因表達情況分析和基因功能注釋，明確特定條件下相關(guān)基因的表達情況[10]，目前轉(zhuǎn)錄組分析被普遍應(yīng)用到水產(chǎn)動物的生物學(xué)研究中。NORIKO等[11]發(fā)現(xiàn)，背景顏色不僅影響平均紅細(xì)胞血紅蛋白量(mean corpuscular hemoglobin，MCH)水平，而且影響大腦中的促性腺激素釋放激素II(chicken-type gonadotropin-releasing hormone II，cGnRH-II)水平，并表明cGnRH-II可能在調(diào)節(jié)比目魚大腦中的 MCH 神經(jīng)功能、食物攝入方面發(fā)揮作用。周江等[12]研究發(fā)現(xiàn)，星斑川鰈(Platichthysstellatus)腦組織中除促卵泡激素(follicle-stimulating hormone, FSH)外，促性腺激素(gonadotropins hormone, GnH)、促性腺激素釋放激素(gonadotropin-releasing hormone, GnRH)以及催乳素釋放肽受體(prolactin-releasing peptide-receptor, PrRPR)等均在卵巢退化期時期個體的腦組織中上調(diào)表達。劉怡南等[13]發(fā)現(xiàn)，igdcc3、prlh、mch、pitx2等基因在七彩神仙魚(Symphysodonharaldi)的腦組織中呈性別差異表達，可能在該魚腦組織類固醇激素形成及配子發(fā)生的調(diào)節(jié)過程中發(fā)揮重要作用。

本研究利用 Illumina 平臺對綠鰭馬面鲀的腦組織進行轉(zhuǎn)錄組測序，再對測序所得數(shù)據(jù)進行質(zhì)控、比對、差異性分析，然后再進行GO(Gene Ontology)富集分析、KEGG(Kyoto Encyclopedia of Genes and Genomes)富集分析、趨勢分析，以期找到綠鰭馬面鲀腦組織響應(yīng)密度脅迫的差異表達基因和關(guān)鍵通路,為了解綠鰭馬面鲀響應(yīng)密度脅迫的分子機制及綠鰭馬面鲀的集約化養(yǎng)殖提供參考依據(jù)。

1 材料與方法

1.1 實驗材料

實驗用綠鰭馬面鲀?nèi)∽灾袊a(chǎn)科學(xué)研究院東海水產(chǎn)研究所贛榆研究中心，系人工自繁。選取體質(zhì)健壯、規(guī)格較為一致且平均體質(zhì)量為(9.5± 0.5)g的個體進行實驗。

實驗前，將綠鰭馬面鲀放在水容量40 m3的室內(nèi)水泥池中暫養(yǎng)一周，密度為50尾·m-3，將此池樣本作為對照組；實驗時，從大池將魚分至水容量1 m3的聚乙烯微流水養(yǎng)殖桶中，設(shè)置2個密度試驗組，即100尾·m-3的中密度組和500尾·m-3的高密度組，每個密度設(shè)置3個重復(fù)。

實驗過程中，水溫保持(17±1)℃，pH為8.20±0.90，采用微流水飽和充氣養(yǎng)殖，溶氧保持在(6.7±0.4)mg·L-1，換水率120%。每日于8∶00、12∶00、15∶00分別投餌1次，總投餌量為魚平均體質(zhì)量的3%。實驗周期50 d。

1.2 組織樣本采集

首先，取對照組大池中的綠鰭馬面鲀腦組織；然后，在實驗開始后的第25天和第50天時，分別在兩個密度試驗組的養(yǎng)殖桶中取3尾魚，用丁香酚對實驗魚進行麻醉、致暈、致死，然后將魚撈起放置于干凈的解剖盤中，用酒精擦拭魚體表面，隨后用解剖剪剪開魚腦上方的骨骼，用鑷子取出腦組織。腦組織取出后，快速切成厚度<0.5 cm的任意片狀，迅速浸入10倍體積 RNAwait 液(Solarbio)中，然后按照RNAwait液的使用說明，將其4 ℃暫存一晝夜，再放置于-20 ℃冰箱冷凍保存，帶回實驗室進行后期實驗。為便于后續(xù)的敘述，各組代號的含義如下：1)CB指對照組的腦組織；2)D1B指第25天時100尾·m-3的腦組織；3)D2B指第25天時500尾·m-3的腦組織；4)d1B指第50天時100 尾·m-3的腦組織；5)d2B指第50天時500尾·m-3的腦組織。

1.3 RNA的提取與檢測、文庫構(gòu)建及測序

取出冷凍的綠鰭馬面鲀腦組織，采用Trizol法對腦組織進行RNA提取，檢測合格后得到該魚腦組織的總RNA。通過帶有Oligo(dT)的磁珠富集具有polyA尾巴的真核mRNA后，用緩沖液把mRNA打斷。以片段化的mRNA為模版、隨機寡核苷酸為引物，在M-MuLV逆轉(zhuǎn)錄酶體系中合成cDNA第一條鏈，隨后用RNaseH降解RNA鏈，并在DNA polymerase I 體系下,以dNTPs為原料合成cDNA第二條鏈。純化后的雙鏈cDNA經(jīng)過末端修復(fù)、加A尾并連接測序接頭，用AMPure XP beads篩選200 bp左右的cDNA，進行PCR擴增并再次使用AMPure XP beads純化PCR產(chǎn)物，最終獲得文庫。構(gòu)建好的文庫利用Illumina HiSeq平臺進行測序。

1.4 測序數(shù)據(jù)組裝及基因功能注釋

將測序獲得的原始數(shù)據(jù)(raw reads)利用fastp[14]進行質(zhì)控、去除含接頭reads、去除含N比例大于10%的reads、去除全部都是A堿基的reads、去除低質(zhì)量reads(質(zhì)量值Q≤20的堿基數(shù)占整條read的50%以上)。得到高質(zhì)量干凈序列(clean reads)進行后續(xù)分析。將質(zhì)控后的過濾數(shù)據(jù)利用HISAT2軟件與綠鰭馬面鲀的參考基因組進行比對[15]，以便進行后續(xù)的基因表達和差異分析。

1.5 基因表達和差異分析

本研究采用 FPKM(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced)方法計算基因的表達量[16]。再用DESeq2[17]軟件包篩選差異表達基因, 篩選閥值為Fold change>2和Q值<0.05，且Q值越小基因表達量差異越顯著[18]。最后，進行差異表達基因的GO和 KEGG富集分析(false discovery rate, FDR< 0.05)。

2 結(jié)果與分析

2.1 腦轉(zhuǎn)錄組測序數(shù)據(jù)

運用 Illumina HiSeq 高通量測序平臺分別對不同密度組綠鰭馬面鲀腦組織進行了轉(zhuǎn)錄組測序，共得到100.3 Gb RawData，過濾后得到98.9 Gb CleanData。堿基質(zhì)量及組成分析顯示，每組GC含量區(qū)間為48.46%～50.38%，各樣品Q20的堿基質(zhì)量值比例均不小于96.67%，各樣品Q30的堿基質(zhì)量值比例均不小于91.36%(表1)。

表1 綠鰭馬面鲀腦轉(zhuǎn)錄組測序數(shù)據(jù)統(tǒng)計

2.2 基因表達量總體分布

對照組和密度脅迫組綠鰭馬面鲀腦組織間比較能反映它們在基因表達水平上的基本差異(圖1-a～d)。圖1中所示每個點代表一個基因，橫坐標(biāo)表示兩個分組間的差異倍數(shù)對數(shù)值，縱坐標(biāo)表示兩個分組差異的FDR的-Lg值，紅色(group_2相對于group_1表達量上調(diào))和藍(lán)色(表達量下調(diào))的點表示基因的表達量有差異(判斷標(biāo)準(zhǔn)為 FDR <0.05, 且差異倍數(shù)在兩倍以上)，黑色的點為沒有差異。其中，越靠近兩端的基因，差異程度越大。在CB-vs-D1B中，有197個差異基因上調(diào)，418個差異基因下調(diào)；在CB-vs-D2B中，有114個差異基因上調(diào)，301個差異基因下調(diào)；在CB-vs-d1B中，有317個差異基因上調(diào)，491個差異基因下調(diào)；在CB-vs-d2B中，有311個差異基因上調(diào)，608個差異基因下調(diào)(圖1-e)。

圖1 密度脅迫組和對照組差異基因表達狀況

2.3 GO功能富集分析

將轉(zhuǎn)錄組信息與 GO(Gene Ontology)數(shù)據(jù)庫作比對，分別注釋到了生物過程(biological process, BP)、分子功能(molecular function, MF)、細(xì)胞組分(cellular component, CC)3個亞類(圖2)，其中紅色代表上調(diào)表達，藍(lán)色代表下調(diào)表達。在生物過程亞類中，得到最顯著注釋的條目依次為：細(xì)胞過程、單組織過程、代謝過程、生物調(diào)控、生物過程的調(diào)節(jié)；在分子功能亞類中, 得到最顯著注釋的條目依次為：結(jié)合、催化活性；在細(xì)胞組分亞類中, 得到最顯著注釋的條目依次為：細(xì)胞組分、細(xì)胞、細(xì)胞膜。但在CB-vs-D2B中，顯著差異注釋的條目相較其他3組更少。

圖2 差異表達基因的 GO 功能富集分析

2.4 KEGG 功能富集分析

用KEGG數(shù)據(jù)庫富集出在密度脅迫后表現(xiàn)出顯著差異的代謝途徑繪制KEGG富集氣泡圖(圖3)，其中Q值為多重假設(shè)檢驗校正之后的P值。利用Q值最小的前20個pathway來作圖，縱坐標(biāo)為pathway，橫坐標(biāo)為富集因子(該pathway中差異基因數(shù)量除以該pathway中所有數(shù)量)，大小表示數(shù)量多少，顏色越紅Q值越小。富集分析的結(jié)果表明：在CB-vs-D1B中，差異表達的基因富集到267條通路上，其中富集程度最高的通路依次是：氧化磷酸化、晝夜節(jié)律、產(chǎn)熱、病毒蛋白與細(xì)胞因子和細(xì)胞因子受體的相互作用(圖3-a)。在CB-vs-D2B中，差異表達的基因富集到230條通路上，其中富集程度最高的通路是：晝夜節(jié)律(圖3-b)。在CB-vs-d1B中，差異表達的基因富集到278條通路上，其中富集程度最高的通路是：晝夜節(jié)律、蛋白質(zhì)消化吸收(圖3-c)。在CB-vs-d2B中，差異表達的基因富集到269條通路上，其中富集程度最高的通路是：晝夜節(jié)律、賴氨酸降解、MAPK信號通路、神經(jīng)活性配體-受體相互作用(圖3-d)。

2.5 趨勢分析

趨勢分析(trend analysis)是針對多個連續(xù)型樣本(至少3個)的特點(樣本間包含特定的時間、空間或處理劑量大小順序)而對基因的表達模式(在多階段中表達曲線的形狀)進行聚類的方法。然后從聚類結(jié)果中挑選符合一定生物學(xué)特性(如表達量持續(xù)上升)的基因集。

2.5.1 100尾·m-3組的趨勢分析

來自CB_D1B_d1B的1 520個基因被分配到8個趨勢中(圖4)，其中 Profile 2、Profile 1、Profile 0和Profile 6為顯著趨勢。在Profile 2中，有361個基因表達先下調(diào)后恢復(fù)到脅迫前的狀態(tài)，KEGG通路分析顯示其中最主要富集的通路是白細(xì)胞跨內(nèi)皮遷移和細(xì)胞凋亡。在Profile 1中，有273個基因表達先下調(diào)后趨平，KEGG通路分析顯示其中最主要富集的通路是晝夜節(jié)律。在Profile 0中，有271個基因表達持續(xù)性下調(diào)，KEGG通路分析顯示其中最主要富集的通路是MAPK信號通路、晝夜節(jié)律、皮質(zhì)醇合成與分泌、鈣信號通路、ECM-受體相互作用。在Profile 6中，有234個基因表達情況先上調(diào)后趨平，KEGG通路分析顯示其中最主要富集的通路是氧化磷酸化、產(chǎn)熱、代謝途徑、糖酵解/糖異生、碳代謝、氨基酸的生物合成、核糖體和RNA降解。

圖4 CB_D1B_d1B所有趨勢總圖

2.5.2 500尾·m-3組的趨勢分析

來自CB_D2B_d2B的1 547個基因被分配到8個趨勢中(圖5)，其中Profile 2、Profile 0、Profile 4和Profile1為顯著趨勢。在Profile 2中，有363個基因表達先下調(diào)后恢復(fù)到脅迫前的狀態(tài)，KEGG通路分析顯示其中最主要富集的通路是Toll樣受體信號通路。在Profile 0中，有349個基因表達持續(xù)性下調(diào)，KEGG通路分析顯示其中最主要富集的通路是MAPK信號通路、cAMP 信號通路、多巴胺能突觸、鈣信號通路和谷氨酸能突觸。在Profile 4中，有178個基因表達先不變再上調(diào)，KEGG通路分析顯示其中最主要富集的通路是PPAR信號通路、過氧化物酶體、氮代謝和類固醇生物合成。在Profile 1中，有168個基因表達先下調(diào)再趨平，KEGG通路分析顯示其中最主要富集的通路是晝夜節(jié)律。

圖5 CB_D2B_d2B所有趨勢總圖

3 討論

為了滿足人類社會需求，水產(chǎn)動物高密度養(yǎng)殖逐漸成為降低養(yǎng)殖成本、提高養(yǎng)殖產(chǎn)量的重要養(yǎng)殖模式，然而過度追求高密度養(yǎng)殖會造成種群內(nèi)對資源的掠奪加劇，最終可能造成水產(chǎn)動物生存能力的降低[19]。當(dāng)密度發(fā)生改變時，水產(chǎn)動物體內(nèi)并非單個的基因或者信號通路發(fā)生變化，而是整個調(diào)控網(wǎng)絡(luò)均會發(fā)生改變[20]。目前，關(guān)于綠鰭馬面鲀的的研究多集中在生長發(fā)育、攝食習(xí)性、養(yǎng)殖技術(shù)、營養(yǎng)成分和遺傳結(jié)構(gòu)分析[21-25]。近些年，轉(zhuǎn)錄組技術(shù)在多種水產(chǎn)動物研究中均得到廣泛應(yīng)用，其在生物進化研究、免疫應(yīng)答反應(yīng)和毒理學(xué)等方面取得了大量的研究進展[26]。為了闡明綠鰭馬面鲀腦組織的基因表達模式，探究明顯富集通路的關(guān)鍵基因，本研究首次構(gòu)建了綠鰭馬面鲀腦組織的cDNA文庫，并通過Illumina平臺進行了高通量測序分析。

3.1 腦組織差異表達基因數(shù)目和GO富集

對綠鰭馬面鲀腦組織的不同組進行轉(zhuǎn)錄組測序、組裝、注釋，發(fā)現(xiàn)在CB-vs-D1B中，有197個差異基因上調(diào)，418個差異基因下調(diào)；在CB-vs-D2B中，有114個差異基因上調(diào)，301個差異基因下調(diào)；在CB-vs-d1B中，有317個差異基因上調(diào)，491個差異基因下調(diào)；在CB-vs-d2B中，有311個差異基因上調(diào)，608個差異基因下調(diào)。表明綠鰭馬面鲀在受到密度脅迫后可能激活了細(xì)胞的代謝，以應(yīng)答應(yīng)激環(huán)境下對綠鰭馬面鲀體內(nèi)造成的傷害。

在GO功能富集分析中，注釋得到的最顯著的條目有細(xì)胞過程、單組織過程、代謝過程、結(jié)合和催化活性等部分，表明綠鰭馬面鲀腦組織的此類生理過程受到了顯著影響。但與其他組相比，CB-vs-D2B中的顯著差異注釋的條目相較于其他3組相對較少，可能是相對其他組來說，D2B在密度脅迫下受到的影響更小。

3.2 各組中均得到富集的通路

在對綠鰭馬面鲀腦組織進行KEGG分析后發(fā)現(xiàn)，各個組中的晝夜節(jié)律通路都得到了顯著富集。晝夜節(jié)律又稱日夜節(jié)律，是生物體內(nèi)的生物鐘面對外界信號影響時所作出的一種內(nèi)源性調(diào)控應(yīng)答反應(yīng)機制，從而使生物的行為節(jié)律以及內(nèi)分泌系統(tǒng)與外界信號同步, 因此是機體對外界環(huán)境信息的一種適應(yīng)性應(yīng)答[27]。clock基因編碼的蛋白質(zhì)在調(diào)節(jié)晝夜節(jié)律中起到核心作用。該蛋白具有組蛋白乙酰轉(zhuǎn)移酶活性，是堿性-螺旋-環(huán)-螺旋家族的轉(zhuǎn)錄因子，在DNA損傷、細(xì)胞凋亡等過程具有重要作用[28]。cry是一種蛋白質(zhì)編碼基因，與cry相關(guān)的疾病包括延遲睡眠階段障礙和睡眠障礙，其相關(guān)通路包括BMAL1-CLOCK、NPAS2，激活晝夜節(jié)律基因的表達與該基因相關(guān)的GO注釋包括核苷酸結(jié)合和轉(zhuǎn)錄因子結(jié)合[29]。這些與生物體晝夜節(jié)律相關(guān)的基因異常表達可能代表著密度脅迫影響了綠鰭馬面鲀腦組織對于生物鐘的調(diào)控，且隨著時間改變未恢復(fù)到正常水平。

3.3 趨勢分析中的關(guān)鍵通路

在綠鰭馬面鲀腦組織的趨勢分析中，發(fā)現(xiàn)MAPK信號通路在CB_D1B_d1B和CB_D2B_d2B均持續(xù)性下調(diào)表達，是真核細(xì)胞介導(dǎo)細(xì)胞外信號到細(xì)胞內(nèi)反應(yīng)的重要信號轉(zhuǎn)導(dǎo)系統(tǒng)，調(diào)節(jié)細(xì)胞的生長、分化、凋亡和死亡等多種生理過程。其中bdnf基因是編碼神經(jīng)生長因子蛋白質(zhì)家族的成員，選擇性剪接產(chǎn)生多個轉(zhuǎn)錄變體，其中至少一個編碼前體蛋白，該前體蛋白經(jīng)過蛋白水解加工產(chǎn)生成熟蛋白[30]。rps6ka3是一種蛋白質(zhì)編碼基因，該基因是編碼絲氨酸/蘇氨酸激酶的 RSK家族的成員[31]。il1r1基因編碼屬于白細(xì)胞介素-1 受體家族的細(xì)胞因子受體，它是參與許多細(xì)胞因子誘導(dǎo)的免疫和炎癥反應(yīng)的重要介質(zhì)[32]。由此看來，MAPK信號通路的持續(xù)性下調(diào)表達對綠鰭馬面鲀的蛋白質(zhì)編碼造成了不利影響，并可能造成該魚的抗病免疫功能降低。

另外，鈣離子信號通路也在CB_D1B_d1B和CB_D2B_d2B中持續(xù)性下調(diào)表達，鈣離子信號通路可以調(diào)節(jié)有氧代謝[33]，其中htr5a是一種蛋白質(zhì)編碼基因，起神經(jīng)遞質(zhì)、有絲分裂原和激素的作用，G蛋白介導(dǎo)這種受體的活性，它可能通過調(diào)節(jié)細(xì)胞內(nèi)Ca2+水平而發(fā)揮作用，可能與精神疾病有關(guān)[34]。ryr2基因編碼主要位于內(nèi)質(zhì)網(wǎng)膜上的RYR-2蛋白，RYR-2作為同四聚體的一個亞基，構(gòu)成了調(diào)節(jié)細(xì)胞內(nèi)鈣離子濃度的鈣通道，其在心臟和大腦中廣泛分布和高表達[35]。pln基因可逆地抑制心臟肌質(zhì)網(wǎng)中ATP2A2的活性，通過對ATP2A2的影響，調(diào)節(jié)心肌對生理刺激的收縮性，在肌肉放松期間調(diào)節(jié)鈣的再攝取，并在心肌中的鈣穩(wěn)態(tài)中發(fā)揮重要作用[36]。

除去以上幾個通路出現(xiàn)顯著富集外，筆者還關(guān)注到PPAR信號通路、Toll樣受體信號通路。PPAR信號通路和動物免疫防御和抗炎方面有關(guān)，PPARα調(diào)節(jié)因子影響脂蛋白和長鏈脂肪酸代謝[37]。推測本研究中綠鰭馬面鲀在密度脅迫下能量消耗變大，隨后糖酵解/糖異生增強，長時間的密度脅迫導(dǎo)致該魚免疫抗病能力的降低。Toll樣受體可識別不同的病原體相關(guān)分子模式，并在先天性免疫應(yīng)答中扮演著不可或缺的角色。它們是抵御病原體入侵的第一道防線，在炎癥、免疫細(xì)胞調(diào)控、存活和增殖方面發(fā)揮著關(guān)鍵作用[38]。

致謝：綠鰭馬面鲀養(yǎng)殖實驗過程中，得到本所吳艷慶、劉威等同志的幫助;綠鰭馬面鲀?nèi)蚪M序列由水科院黃海所陳四清研究員、邊力博士提供；謹(jǐn)致謝忱。