慢性密度脅迫對綠鰭馬面鲀肌肉轉錄組的影響

程起群 段曉晨 宋丹

(中國水產科學研究院東海水產研究所,上海 200090)

綠鰭馬面鲀(Thamnaconusseptentrionalis)隸屬于鲀形目(Tetraodontiformes)、單角鲀科(Monacanthidae)、馬面鲀屬(Thamnaconus),又稱馬面魚、剝皮魚、老鼠魚等,是我國重要的海產經濟魚類,主要分布在我國沿海地區[1]。近年來,由于過度捕撈、環境污染等因素,綠鰭馬面鲀野生資源呈現衰退的趨勢[2]。開展綠鰭馬面鲀的人工養殖來代替捕撈,是緩解綠鰭馬面鲀野生資源壓力、保障市場供給的有效途徑。

2011年,綠鰭馬面鲀人工繁育技術取得突破性的進展,自此開始了工廠化養殖[3]。養殖密度是水產養殖過程中的關鍵性限制因素。慢性密度脅迫是指在一定密度養殖條件下,魚類長時間暴露在相對狹小的空間中,其生長受密度的慢性壓力影響,表現為生長速度減緩、生理狀況和健康發生變化以及應激反應激活,最終影響其福祉和生存能力[4-5]。有研究表明,養殖密度顯著影響了中國明對蝦(Fenneropenaeuschinensis)的生長,不僅增加了碳水化合物和脂類的消耗,還抑制了組織的有氧代謝[6];Sun等[7]的研究顯示,擁擠可以顯著誘導魚的免疫反應,而長期處于密度脅迫下會損害魚類的免疫能力;王楠楠等[8]研究顯示,養殖密度過高或過低都會影響點帶石斑魚(Epinephelusmalabaricus)幼魚的生長和代謝,當密度較高時,魚體攝食較多,排泄物對水環境的影響增加,從而導致成活率較低。目前,關于綠鰭馬面鲀肌肉組織對密度脅迫的分子響應機制尚不明確,有待開展深入研究。

轉錄組學是揭示生物分子機制的重要手段[9],是通過對水產動物進行高通量測序,開展基因功能注釋和基因表達情況的分析[10-11]。目前,轉錄組測序技術在水產動物的研究中已有較多應用。例如,曾爽等[12]對翹嘴鱖(Sinipercachuatsi)極小個體和極大個體的肌肉組織進行了轉錄組測序,發現mhc、mlc、igfr1、ghr2等多個基因在極小個體中的表達量上調,大量差異表達基因在蛋白生成和消化相關通路得到顯著富集,表明它們可能參與了生長調控;李玲等[13]發現,高密度養殖(300尾/m2)增加了斑節對蝦(Penaeusmonodon)額外的能量代謝,降低了蝦的生長性能;閆允君等[14]首次基于轉錄組技術在無脊椎動物中篩選出了由肌肉生長抑制素基因調控的與肌肉生長發育相關的基因;王蘭梅等[15]發現,生長速度較快的福瑞鯉(Cyprinuscarpio),其lamb3和pdk2的基因表達量高于生長速度較慢的,說明該基因能調控福瑞鯉的肌肉生長,且生長速度較快的魚體內脂肪酸和葡萄糖的代謝情況可能發生了改變。但目前尚未發現轉錄組技術應用于綠鰭馬面鲀的研究報道。

本研究利用Illumina平臺對不同養殖密度下的綠鰭馬面鲀的肌肉組織進行高通量測序,再對測序所得數據進行質控、比對、差異性分析,然后再進行GO(Gene Ontology)功能注釋和KEGG(Kyoto Encyclopedia of Genes and Genomes)功能富集分析,以期找到該魚肌肉組織響應密度脅迫的關鍵富集通路和差異表達基因(differentially expressed genes),為探索綠鰭馬面鲀的最適養殖密度,進而開展分子育種研究提供數據參考。

1 材料和方法

1.1 材料

試驗所用綠鰭馬面鲀為中國水產科學院東海水產研究所贛榆研究中心人工繁殖的魚種。選擇體格健壯、無病無傷、平均體質量(9.5±0.5)g的個體進行試驗。

在正式試驗開始前,將綠鰭馬面鲀放在水體容量為40 m3的室內水泥池中暫養1周,放養密度為50尾/m3。正式試驗時,設置2個試驗密度組,即中密度組(100尾/m3)和高密度組(500尾/m3),每個密度設3個重復,模擬在慢性密度脅迫下綠鰭馬面鲀的生存環境。試驗過程中采用微流水對圓柱形養殖桶(底面積0.9 m2、高1.2 m)進行飽和曝氣,日換水率為120%,溶解氧保持在(6.7±0.4)mg/L,水溫保持在(17±1)℃,pH為8.20±0.9。每天投喂3次餌料(8:00、12:00、15:00),投喂量為魚體質量的3%。試驗共持續50 d(2020年11月7日—2020年12月27日)。

1.2 組織樣本采集

試驗開始時,采集暫養池中綠鰭馬面鲀的背部肌肉組織作為對照樣本;試驗開始后,分別于第25天和第50天在每個密度組隨機取3尾魚作為生物學重復。采集樣本時,用丁香酚將試驗魚麻醉,然后將魚撈起置于干凈的解剖盤中,用酒精擦拭魚體表面,采集其背部肌肉組織。

將對照組的肌肉組織標記為CM,養殖第25天時中密度組和高密度組的肌肉組織分別標記為D1M和D2M,養殖第50天時中密度組和高密度組的肌肉組織分別標記為d1M和d2M。獲得肌肉組織后,快速將之處理成<0.5 cm厚的隨機切片并浸入5～10倍體積的RNA保存液(RNAwait)中,然后將樣品置于4 °C冰箱中過夜,再于-20 ℃冷凍保存。

1.3 RNA的提取與檢測、文庫構建及測序

取冷凍的綠鰭馬面鲀肌肉組織,用Trizol法從肌肉組織中提取RNA,檢測符合要求后得到肌肉組織總RNA。用Oligo(dT)磁珠富集具有polyA尾的真核mRNA后,用緩沖液破壞mRNA。將片段化的mRNA作模板,引物為隨機寡核苷酸,在M-MuLV逆轉錄酶系統中合成第一條cDNA鏈,然后用RNase H消化RNA鏈,以dNTP作為DNA聚合酶Ⅰ系統的起始材料,用于合成第二條cDNA鏈。將純化得到的雙鏈cDNA進行末端修復,加A尾,連接測序接頭,最終得到約200 bp的cDNA文庫。委托基迪奧公司使用Illumina HiSeq平臺對構建的文庫進行測序。

1.4 測序數據組裝及基因功能注釋

使用fastp[16]對測序得到的原始數據(raw reads)進行質量控制,刪除包含接頭的、N比例超過10%的、全部是A堿基的以及低質量(Q≤20的堿基數占總數據的50%以上)的數據,獲得定性凈讀數以供進一步分析。使用Hisat2程序[17]將質量控制后的過濾數據與綠鰭馬面鲀的參考基因組進行比較,從而進行基因表達分析和差異分析。

1.5 基因表達和差異分析

在本研究中,使用FPKM(expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced)方法計算基因表達[18]。然后使用DESeq2軟件包[19]篩選差異表達基因,篩選閾值:差異表達倍數(fold change,FC)>2,Q<0.05,Q值越小,基因表達差異越顯著[20]。最后,對GO和KEGG差異表達基因(FDR值<0.05)進行富集分析。

2 結果

2.1 肌肉轉錄組測序數據

使用Illumina HiSeq高通量測序平臺對各組綠鰭馬面鲀肌肉組織樣本進行轉錄組測序,共得到104.5 Gb原始數據,過濾后得到103.3 Gb高質量數據。堿基質量及組成分析表明,每組的GC(guanine-cytosine)含量在43.39%～54.06%,各樣品的Q20堿基質量值比例均不小于97.48%,各樣品Q30的堿基質量值比例均不小于92.91%。

2.2 基因表達量總體分布

綠鰭馬面鲀肌肉組織對照組與試驗組之間的比較,反映出它們在基因表達量上的基本差異(見圖1)。圖中的每個點均代表1個基因,橫坐標為兩個比較組差異倍數的對數值,縱坐標為兩個比較組差異的FDR的負log10值。以FDR值<0.05、差異倍數在2倍以上為判斷標準,紅色代表基因表達量上調,藍色代表表達量下調,灰色點為沒有差異,越靠近兩端的基因,其差異程度越大。在CM-vs-D1M中,有107個差異基因上調,829個差異基因下調;在CM-vs-D2M中,有87個差異基因上調,762個差異基因下調;在CM-vs-d1M中,有236個差異基因上調,970個差異基因下調;在CM-vs-d2M中,有342個差異基因上調,876個差異基因下調(見圖2)。

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

圖2 密度脅迫組與對照組差異基因表達匯總

2.3 GO功能富集分析

將轉錄組信息與GO數據庫作比對,分別注釋到了生物過程(biological process,BP)、分子功能(molecular function,MF)和細胞組分(cellular component,CC)3個亞類(見圖3)。圖3中,橙色代表上調表達,藍色代表下調表達。在“生物過程”亞類中,得到最顯著注釋的條目依次為:細胞過程、代謝過程、單組織過程、生物調控、生物過程的調節;在“分子功能”亞類中,得到最顯著注釋的條目依次為:結合、催化活性;在“細胞組分”亞類中,得到最顯著注釋的條目依次為:細胞、細胞組分、細胞器。在4個比較組中,下調基因的數目均遠大于上調基因的數目,其中D1M和D2M兩個組與對照組的差異尤為明顯。

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

2.4 KEGG功能富集分析

利用KEGG數據庫富集出在密度脅迫后表現出顯著差異的代謝途徑,繪制KEGG富集氣泡圖(見圖4),其中Q值為多重假設檢驗校正之后的P值。利用Q值最小的前20個通路來繪圖,縱坐標為通路,橫坐標為富集因子(即該通路中差異基因數量除以該通路中所有基因數量),大小表示數量多少,顏色越紅,則Q值越小。富集分析結果表明:在CM-vs-D1M中,差異表達的基因富集到288條通路上,其中富集程度最高的通路為:晝夜節律、AMPK信號通路、破骨細胞分化、內質網中的蛋白質加工和MAPK信號通路(見圖4-a)。在CM-vs-D2M中,差異表達的基因富集到290條通路上,其中富集程度最高的通路為:破骨細胞分化、MAPK信號通路、AMPK信號通路和Th17細胞分化(見圖4-b)。在CM-vs-d1M中,差異表達的基因富集到300條通路上,其中富集程度最高的通路為:蛋白質消化吸收、晝夜節律、黏著斑、真核生物中的核糖體生物發生、PI3K-Akt信號通路和HIF-1信號通路(見圖4-c)。在CM-vs-d2M中,差異表達的基因富集到293條通路上,其中富集程度最高的通路為:黏著斑、晝夜節律、松弛素信號通路、PI3K-Akt信號通路和MAPK信號通路(見圖4-d)。

圖4 差異表達基因的 KEGG 功能富集分析

2.5 趨勢分析

趨勢分析(trend analysis)是針對3個以上連續型樣本的特點,例如樣本間包含特定的空間、時間或處理劑量順序而對基因的表達模式進行聚類的方法。趨勢分析之后,從聚類結果中挑選符合一定生物學特性(如表達量持續上升)的基因集。

2.5.1 中密度組趨勢分析

來自CM_D1M_d1M的1 771個基因被分配到

8個趨勢中(Profile),其中Profile 1、Profile 0為顯著趨勢。在Profile 1中,有692個基因表達量先下調后趨平,KEGG通路分析顯示,其中最主要富集的通路是內質網中的蛋白質加工、破骨細胞分化、MAPK信號通路。在Profile 0中,有421個基因表達量持續性下調,KEGG通路分析顯示,其中最主要富集的通路是真核生物中的核糖體生物發生、黏著斑、PI3K-Akt信號通路、黏著連接(見圖5)。

注:圖中紅色背景的Profile 1和Profile 0代表顯著趨勢,分別包含692和421個基因,顯示表達量先下調后趨平(Profile 1)與持續性下調(Profile 0)。曲線反映基因表達隨時間的變化。

2.5.2 高密度組趨勢分析

來自CM_D2M_d2M的2 017個基因被分配到8個趨勢中(見圖6),其中Profile 1、Profile 2、Profile 0和Profile 3為顯著趨勢。在Profile 1中,有552個基因表達量先下調后趨平,KEGG通路分析顯示,其中最主要富集的通路是MAPK信號通路、Th17細胞分化、真核生物中的核糖體生物發生、松弛素細胞通路、RNA降解。在Profile 2中,有378個基因表達量先下調后恢復到脅迫前的狀態,KEGG通路分析顯示,其中最主要富集的通路是補體和凝血級聯、破骨細胞分化、病毒蛋白與細胞因子和細胞因子受體的相互作用。在Profile 0中,有364個基因表達量持續性下調,KEGG通路分析顯示,其中最主要富集的通路是黏著斑、蛋白質消化吸收、PI3K-Akt信號通路、cGMP-PKG信號通路。在Profile 3中,有271個基因表達量先不變再下調,KEGG通路分析顯示,其中最主要富集的通路是ECM-受體相互作用、黏著斑、PI3K-Akt信號通路。

注:圖中紅色背景的Profile 1與Profile 0代表顯著基因表達趨勢,Profile 1包含552個基因,表現為表達量下調后趨于穩定;Profile 0包含364個基因,其表達量持續下降。綠色背景的Profile 2和藍色背景的Profile 3分別包含378和271個基因。

3 討論

3.1 肌肉組織差異表達基因數目和GO富集

對不同組的綠鰭馬面鲀肌肉組織進行轉錄組測序、組裝、注釋,發現在CM-vs-D1M中,有107個差異基因上調,829個差異基因下調;在CM-vs-D2M中,有87個差異基因上調,762個差異基因下調;在CM-vs-d1M中,有236個差異基因上調,970個差異基因下調;在CM-vs-d2M中,有342個差異基因上調,876個差異基因下調。這一結果表明綠鰭馬面鲀在受到密度脅迫后可能激活或改變了細胞的生理活動,以應答脅迫環境對其機體內分子層面的影響。

在GO功能富集分析中,注釋得到的最顯著條目有:細胞過程、代謝過程、單組織過程、生物調控、生物過程的調節、結合、催化活性、細胞部分和細胞等,表明綠鰭馬面鲀肌肉組織的這些生理過程受到了顯著影響。在這4個比較組中,下調基因數目均遠大于上調基因數目,說明綠鰭馬面鲀在受到脅迫后通過一系列基因表達情況的改變來適應環境,且隨著時間推移,基因表達情況發生了改變。

3.2 富集通路分析

在對綠鰭馬面鲀肌肉組織進行KEGG分析后發現,在4個試驗組以及趨勢分析中,破骨細胞分化、MAPK信號通路、黏著斑、PI3K-Akt信號通路均得到了顯著富集。

破骨細胞分化通路的富集表明密度脅迫可能促進了骨細胞的分化和活躍化。本研究中發現,破骨細胞分化通路中,有多個關鍵基因得到了富集,包括nfatc1、tgfbr2、map3k14、socs3、jak1、pparg和fos等。nfatc1誘導了一系列與骨細胞活性相關的基因[21]。tgfbr2編碼的蛋白質在調控細胞周期、細胞增殖、傷口愈合和免疫抑制等相關基因的轉錄中發揮著關鍵作用[22]。綠鰭馬面鲀可能通過激活破骨細胞分化通路來適應環境變化和維持骨骼健康。map3k14介導刺激信號從細胞膜到細胞核或其他靶點,能夠調節細胞的增殖、分化和細胞周期等[23]。在密度脅迫下,這個通路可能會被激活,導致細胞生物學響應的變化。socs3通過抑制炎癥相關信號通路核轉錄因子和絲裂原活化蛋白激酶的過度活化,調節促炎細胞因子的過度釋放,進而抑制炎癥反應[24-26]。jak1在影響介導炎癥、上皮重塑的基因表達中起著至關重要的作用[[27]]。pparg基因則在調節脂肪組織中脂肪酸代謝、影響脂肪細胞內脂肪酸沉積和脂滴形成中扮演主要角色,與脂肪細胞的生長分化也密切相關[28-29]。fos基因家族編碼了FOS蛋白,這些蛋白被認為是細胞增殖、分化和轉化的調節因子,在某些情況下,fos基因的表達也與細胞凋亡相關,人體中與fos相關的疾病包括先天性全身性脂肪營養不良和增殖性筋膜炎[30]。這些基因的富集表明密度脅迫對綠鰭馬面鲀可能導致骨骼發育、細胞增殖、分化以及炎癥反應等多方面生物學過程的調節。因此猜測,在綠鰭馬面鲀中存在類似的機制,這些機制對其適應性反應至關重要。

黏著斑基于整合素與肌動蛋白的相互作用,將細胞連接到細胞外基質,并位于上皮細胞的基底面[31]。黏著斑通路的富集表明,密度脅迫可能通過影響細胞的黏附和形狀來調節細胞活動。在黏著斑通路中,prkca是一種蛋白質激酶編碼基因,這種激酶在許多不同的細胞過程中發揮作用,例如細胞黏附、細胞轉化、細胞周期檢查點和細胞體積控制等[32]。進一步的研究可能可以揭示,prkca是如何調節細胞黏附和形態維持從而影響整個組織的生物學功能的。此外,lna和flnc基因編碼的蛋白質對肌動蛋白的組裝和肌肉細胞骨架的重塑至關重要[33]。fn1基因編碼纖連蛋白,纖連蛋白與細胞表面和各種化合物(包括膠原蛋白、纖維蛋白、肌動蛋白等)結合,并參與細胞黏附、細胞運動、調理作用、傷口愈合和細胞形狀維持[34]。tln1基因編碼一種細胞骨架蛋白,該蛋白集中在細胞基質和細胞間接觸區域,編碼的蛋白質在肌動蛋白絲的組裝以及各種細胞類型(包括成纖維細胞和破骨細胞)的擴散和遷移中起重要作用[35]。黏著斑通路的富集暗示這兩個基因可能調節綠鰭馬面鲀肌肉組織中的細胞動力學,從而影響其生物學過程。由thbs1基因編碼的蛋白質是同源三聚體二硫鍵連接蛋白的亞基,這是一種介導細胞間和細胞間相互作用的黏附性糖蛋白[36]。這些富集的基因可能在細胞間的相互作用以及肌肉組織整體的結構和功能維護中發揮著關鍵的協同作用,可以進一步解釋綠鰭馬面鲀肌肉組織是如何應對密度脅迫等外部因素的。

PI3K-Akt信號通路是細胞中一個經典的信號通路,在胰島素的刺激以及脂肪細胞生存中發揮著重要作用[37]。fn1基因在PI3K-Akt信號通路中發揮著同樣的作用[34]。jak1基因編碼膜蛋白,該基因在影響介導炎癥、上皮重塑的基因表達中起著至關重要的作用[38]。thbs1基因編碼的蛋白質是二硫鍵連接的同源三聚體蛋白質的亞基。這種蛋白質是一種黏附性糖蛋白,可介導細胞與細胞和細胞與基質的相互作用,可以與纖維蛋白原、纖連蛋白、層黏連蛋白等結合[39]。

MAPK信號通路也在密度脅迫后的綠鰭馬面鲀的肌肉組織中被富集,它編碼一組進化上保守的絲氨酸/蘇氨酸蛋白激酶,這些激酶被一系列細胞外刺激信號激活并介導從細胞膜到細胞核的信號轉導,它們調節許多生理活動,如炎癥、細胞凋亡、癌變、腫瘤細胞的侵襲和轉移等[40]。在MAPK信號通路和破骨細胞分化通路中,il1r1基因編碼屬于家族性白細胞介素-1受體的細胞因子受體,該受體是參與細胞因子誘導的免疫和炎癥反應的關鍵介質。由此可見,慢性密度脅迫導致綠鰭馬面鲀肌肉組織多種生理生化活動發生改變。多種信號通路共同調節綠鰭馬面鲀的細胞形狀、細胞周期、炎癥反應、細胞黏附等生理活動。

4 小結

研究結果表明,在慢性養殖密度脅迫情形下,nfatc1、tgfbr2、map3k14、prkca等基因在綠鰭馬面鲀肌肉組織中呈現差異性表達,可能對該魚的細胞周期、細胞增殖、細胞分化和傷口愈合產生明顯影響。socs3、jak1、il1r1等基因的差異性表達可能對綠鰭馬面鲀肌肉組織的蛋白質編碼產生影響,進而影響該魚的炎癥反應、免疫抗病能力、存活等方面。pparg、fos、fosb、fosl2的差異性表達影響綠鰭馬面鲀肌肉組織的脂肪細胞生長、分化等方面。flna、flnc、fn1、tln1、thbs1的差異性表達則可能影響該魚肌肉組織的細胞形狀和遷移、細胞黏附。本研究為揭示綠鰭馬面鲀肌肉組織在密度脅迫下的基因差異表達及信號通路機制提供了基礎信息,可為今后更好地開展綠鰭馬面鲀的集約化養殖和分子育種提供理論積累。

致謝:在綠鰭馬面鲀養殖試驗過程中,得到本所吳艷慶、劉威等同志的幫助;綠鰭馬面鲀全基因組序列由水科院黃海所陳四清研究員、邊力博士提供。在此謹致謝忱。