基于皮膚轉錄組的半滑舌鰨無眼側黑化相關基因的篩選鑒定

趙娜 ，何曉旭 ，賈磊，朱春華，張博*

(1.廣東海洋大學 水產學院，廣東 湛江 524025；2.天津市水產研究所，天津 300457；3.上海海洋大學 水產種質資源發掘與利用教育部重點實驗室，上海 201306；4.天津市農學院，天津 300384)

1 引言

動物體色具有迷惑或躲避天敵等自我保護功能，不同體色對魚類的生存也有著特殊的意義。魚類在生長發育過程中體色的變化主要受遺傳因素的影響，此外，外部環境[1-3]、飼料營養[4]和神經內分泌物[5-7]也會影響魚類的體色變化[8-9]。多數色素細胞的遺傳是保守的[10-11]，如一些與色素相關的基因先后在小鼠（Mus musculus）的皮膚和斑馬魚（Danio rerio）的皮膚中均有發現，而mitf、tyr、dct、sox10和ednr等[12-15]基因在皮膚中的生長發育表達特征也非常相似。此外，在硬骨魚中發現的slc24a5基因也參與哺乳動物體色或毛發中色素細胞的生長和發育[16]。在魚類體色生成過程中，黑色素的產生和擴散起主導作用，其中酪氨酸酶基因家族中的酪氨酸酶（TYR）是主效基因[17]。TYR 受多個基因控制，酪氨酸酶受體基因1（tyr1）、性別決定區Y-box 10 基因（sox10）[18]、多巴素異構酶基因（dct）[19]、黑皮質素受體-1（mcir）對酪氨酸酶的產生均具有調節作用[20-22]。

研究表明，飼養密度[23]、水溫[24]、光照強度[24]能顯著地影響比目魚（Pleuronectiformes）無眼側黑化的發生。科學的養殖管理方法如營養成分的攝入控制和親本篩選[25]也能有效地減少比目魚體色異常的發生。Nakamura 等[26]認為，缺乏光敏物質會導致黑色素合成不足，如核黃素、類胡蘿卜素、維生素A 和維生素D，而核黃素被認為是其中最重要的因素。研究表明[27]，牙鲆（Paralichthysolivaceus）的白化病是視紫紅質缺乏的結果，其產生依賴于DHA。而Devresse 等[28]認為，色素沉著似乎與DHA 無關，而是與DHA 和二十碳五烯酸的比例有關。胰蛋白酶、維甲酸相關酶、維甲酸受體及其結合受體（TR、VDR、PPAR 等）的突變及其對仔魚色素細胞分化的調控，可能影響魚體兩側色素細胞的命運，進而影響體細胞色素異常[29-30]。不同環境背景色養殖實驗中發現，條斑星鰈（Verasper moseri）在淺色養殖環境下mch-R2的mRNA 表達量低于黑色養殖環境，但是兩種環境中mch-R1的mRNA 沒有變化，表明了mch-R2參與了不同環境背景下體色變化過程，mch-R可能對黑色素聚集激素（MCH）有抑制調節作用[31]。無眼側黑化消褪的牙鲆垂體MCH mRNA表達水平顯著高于無眼側黑化個體，同時與無眼側正常個體無顯著差異，表明了MCH 及其基因均參與了無眼側黑化消褪的調控過程[32]。在魚類中，tyr基因突變導致青鳉（Oryzias latipes）的黑色素細胞異常，皮膚呈現橙色[33]。虹鱒魚（Oncorhynchusmykiss）中的tyr基因突變顯示嵌合白化（偽白化）和完全白化現象[34]。黃顙（Pelteobagrus fulvidraco）的tyr基因突變則顯示黑色素缺乏，呈現體色部分白化[35]。綜上所述，遺傳和環境因素的相互作用影響了魚類的體色異常。

半滑舌鰨（Cynoglossus semilaevis）的工廠化養殖過程中經常出現高比例體色異常的苗種，即無眼側出現黑色素的黑化現象[36-37]，主要表現為無眼側部分覆蓋斑塊色素，有時甚至完全呈黑底，不同的養殖環境會存在不同程度的黑化現象。體色異常極大地影響半滑舌鰨的銷售價格，長期困擾養殖企業，卻一直未能得到有效解決[38]。針對半滑舌鰨無眼側黑化，影響體色變化的新基因有待挖掘和鑒定。轉錄組是目前研究基因表達的重要手段，是篩選新功能基因的有效方法。本文選取半滑舌鰨體色異常皮膚樣本和正常皮膚樣本，進行轉錄組測序的比較分析。以期找出與半滑舌鰨無眼側黑化相關的表達顯著變化的功能基因，為半滑舌鰨體色異常的機制提供分子水平的理論支持。

2 材料和方法

2.1 魚類取樣準備

從天津乾海源水產養殖有限公司獲得3 尾無眼側黑化的工廠化養殖半滑舌鰨，平均體重為（0.61±0.02）kg，平均體長為（28.2±0.1）cm。用濃度為120 mg/L的間氨基苯甲酸乙酯甲磺酸鹽（MS-222）（北京格林恒興生物科技有限公司）將魚麻醉后，使用滅菌的手術剪和鑷子剪開實驗魚的皮膚并采集樣本，獲取無眼側色素沉著的皮膚組織和非色素沉著的正常皮膚組織（圖1），并在-80℃下保存，直至提取RNA。

圖1 半滑舌鰨無眼側黑化魚體色特征Fig.1 The body color feature of the Cynoglossus semilaevis with no eye side blackened fish

2.2 RNA 測序、數據過濾及生物信息學分析

使用Trizol 試劑（日本Takara）提取總RNA。安捷倫2100 生物分析儀（安捷倫科技公司，美國圣克拉拉）用于評估RNA 質量和濃度。無眼側黑化皮膚組織中的樣品標記為csc1、csc2 和csc3；正常皮膚樣品標記為csh1、csh2 和csh3。測序平臺為Illumina HiSeq?2500（Illumina 公司，圣地亞哥，加利福尼亞州，美國）。用帶Oligo（dT）的磁珠富集mRNA，再以片段化后的mRNA 為模板，用六堿基隨機引物合成cDNA 第一鏈，并加入緩沖液、dNTPs（A、U、G、C）、核糖核酸酶H 和DNA 聚合酶I 合成cDNA 第二鏈；再經過磁珠純化并加EB 緩沖液洗脫，對洗脫之后的雙鏈cDNA進行末端修復、加堿基A、加測序接頭處理，然后用磁珠進行片段大小選擇、降解含U 鏈，并進行PCR 擴增，從而完成整個文庫制備工作。根據物種的參考基因組和基因信息（GCA_000523025.1），將轉錄組測序所產生的數據比對到參考基因組上。獲得過濾后數據后，將其與參考基因組進行序列比對，獲取測得條目的定位信息。

2.3 差異表達分析和樣品相關性分析

用bowtie2 軟件將過濾后數據比對到參考轉錄本序列，而后利用RSEM 軟件，調用bowtie2 的比對結果進行統計，得到每個樣品比對到每個轉錄本上的測得條目數目，并對其進行FPKM 轉換。FPKM 是指每百萬測得片段中來自某一基因平均每1 000 堿基長度的測得片段數目，其同時考慮了測序深度和基因長度對測得片段計數的影響，是目前最為常用的基因表達水平估算方法。用R 語言包edgeR 進行差異分析，篩選閾值為錯誤發現率（False Discovery Rate，FDR）小于0.05，log2FC（FC 為實驗組與對照組可測量比值）大于1 或log2FC 小于-1。利用皮爾森相關系數進行相關性分析，相關系數越接近1，表明樣品之間表達模式的相似度越高。相關性系數是介于-1～1 之間的實數，當相關性系數介于-1～0 時，表明變量之間存在負相關關系；當相關性系數介于0～1 時，表明變量之間存在正相關關系；當相關性系數為0 時，二者之間不存在相關性。

2.4 差異基因GO 和KEGG 富集及功能篩選

根據分析目的篩選出差異基因后，研究差異基因在Gene Ontology（簡稱GO，http://www.geneontology.org/）中的分布狀況。KEGG（Kyoto Encyclopedia of Genes and Genomes）是有關通路的主要公共數據庫。通路顯著性富集分析以KEGG 通路為單位，應用超幾何檢驗，找出差異基因相對于所有有注釋的基因顯著富集的通路。GO 富集分析方法為Hyper-geometric 分布，與KEGG 富集分析方法一樣，我們選取FDR≤0.05的GO 條目作為顯著富集的GO 條目。FDR≤0.05 的通路定義為在差異表達基因中顯著富集的通路。我們挑選了富集最顯著的20 條通路，繪制差異表達基因KEGG 富集散點圖。KEGG 富集程度通過富集因子、Qvalue 和富集到此通路上的基因個數來衡量。其中富集因子指該通路中富集到的差異基因個數與注釋基因個數的比值。富集因子越大，表示富集的程度越大。Qvalue 是做過多重假設檢驗校正之后的Pvalue，Qvalue 的取值范圍為0～1，越接近0，表示富集越顯著。差異基因功能篩選集中于比目魚體色異常的相關6 條代謝通路，這6 條KEGG 代謝通路如下：花生四烯酸代謝通路、甲狀腺激素合成通路、光轉導通路、黑色素生成通路、視黃醇代謝通路和PPAR 信號通路。

2.5 候選基因的qPCR 定量驗證

用qPCR 方法檢測黑化魚正常皮膚樣本（csh1、csh2 和csh3）和黑化皮膚樣本（csc1、csc2 和csc3）中的候選mRNA 的表達。采用北京百奧創新科技有限公司核酸提取試劑盒和PowerUp-SYBR-Green-Master-Mix 進行mRNA 提取和qPCR。在42℃加熱2 min，在冰浴中孵育2 min，加入提取RNA 的試劑，在42℃加熱15 min，提取總RNA，然后在85℃下保持5 min。將RNA 儲存在-80℃下。cDNA 合成按照PrimeScript?RT reagent Kit 試劑盒操作方法進行，cDNA 儲存在-20℃下備用。使用QuantStudio6 Flex 實時PCR 系統（美國加利福尼亞州賽默飛世爾科學公司）在96 孔板中進行qPCR，反應體積為20 μL：AceQ?qPCR SYBR?Green Master Mix (Without ROX) 10.0 μL，Primer 1(10 μmol/L) 0.4 μL，Primer 2 (10 μmol/L) 0.4 μL，cDNA 1.0 μL，ddH2O 8.2 μL，配制過程在冰上完成。反應程序包括：50℃預熱2 min，95℃預熱2 min，隨后在95°C 下，持續30 s，58°C 持續30 s，95°C 持續15 s，進行40 個循環。選擇β-肌動蛋白的表達水平作為正常化的內源性對照。每個樣品設置3 個技術重復，并使用2-ΔΔCt方法計算mRNA 的相對表達量。

3 結果

3.1 RNA-seq 數據披露

通過RNA-seq，6 個樣本共產生416 464 144 個過濾后數據，平均數據量在10 G 左右（9 208～11 897 M），Q20 和Q30 分別在92%和97%左右。數據結果表明，6 個樣品的測序質量良好，數據量充足，詳細數據如表1 所示，轉錄組數據已經上傳至美國國立生物技術信息中心，數據獲取編號為（PRJNA665385）。

表1 6 個半滑舌鰨皮膚轉錄組測序樣品過濾后數據產出及質控分析Table 1 Clean reads output and quality control analysis of six Cynoglossus semilaevis skin transcriptome sequencing samples

對兩組6 個樣品的RNA-seq 的過濾后數據進行基因組比對統計，6 個樣品的過濾后數據數量分布在30 772 761～39 743 985 之間，唯一比對到參考基因組上的測得條目百分比分布在76.47%～81.36%之間，不能比對到參考基因組的測得條目數目（雙端統計）普遍低于20%（表2）。以上結果進一步說明6 個樣品的轉錄組測序數據質量過關，可作進一步分析。

表2 6 個半滑舌鰨皮膚轉錄組測序樣品基因組比對統計表Table 2 Genomic comparison statistics of six Cynoglossus semilaevis skin transcriptome sequencing samples

3.2 差異表達分析及樣品相關性分析

csc1、csc2 和csc3 樣本之間的相關系數處于0.97～1之間，接近1，表明變量之間存在正相關關系，樣品之間表達模式為強相關。csh1、csh2 和csch3 樣本之間的相關系數處于0.76～1 之間，表明變量之間存在正相關關系，樣品之間表達模式為強相關，但是相關性弱于csc 樣本。csc1、csc2 和csc3 樣本和csh1、csh2和csch3 樣本的相關系數處于0.18～0.45 之間，表明變量之間存在正相關關系，樣品之間表達模式處于弱相關和中度相關（圖2）。

圖2 測序樣品相關性分析熱圖Fig.2 Correlation analysis heat map of sequencing samples

圖3a 顯示，樣本無眼側黑化皮膚csc1、csc2、csc3和無眼側正常皮膚csh1、csh2、csch3 之間差異基因總數1 665 個，其中上調基因497 個，下調基因1 168 個。通過轉錄組測序得到黑化半滑舌鰨的無眼側正常和無眼側黑化皮膚的轉錄組序列信息，對轉錄組數據分析共得到837 540 個位點，其中A/G、C/G、C/T、G/A之間轉換數量最多。通過圖3b 也可以觀察到，csc與csh 兩組之間部分基因的表達量存在明顯差異（圖3）。

圖3 半滑舌鰨正常皮膚與黑化皮膚樣品間差異表達基因火山圖和差異表達基因豐度熱圖Fig.3 Volcano map of the number of differentially expressed genes and heat map of abundance of differentially expressed gene between normal and hypermelanotic skin samples of Cynoglossus semilaevis

3.3 差異表達基因的GO 及 KEGG 富集分析及體色功能差異基因篩選

對無眼側黑化組和無眼側正常組皮膚進行GO富集分析時發現，在生物過程水平上差異基因主要富集在細胞過程、信號組織過程、代謝過程、生物調控、生物學過程這5 個GO 通路；在細胞組分上主要富集在膜、細胞、細胞組分、膜組分、細胞器這5 個GO 通路；在細胞功能上主要富集在連接、催化活性這2 個GO 通路（圖4a）。由KEGG 富集散點圖（圖4b）可見，無眼側黑化和無眼側正常皮膚表達差異顯著的基因參與的通路主要有：丙酮酸代謝、C5-支鏈二元酸代謝、甲烷代謝及苯丙氨酸、酪氨酸和色氨酸的生物合成等。

圖4 半滑舌鰨正常與黑化皮膚轉錄組測序差異基因富集注釋的GO 通路和KEGG 通路Fig.4 GO pathways and KEGG pathways analysis of differential expressed genes enrichment annotations of transcriptome sequencing of normal and hypermelanotic skin of Cynoglossus semilaevis

3.4 定量PCR 驗證結果

從轉錄組富集結果中選取的6 條KEGG 代謝通路如下：花生四烯酸代謝通路、甲狀腺激素合成通路、光轉導通路、黑色素生成通路、視黃醇代謝通路和過氧化物酶體增殖物激活受體信號通路。在以上6 條通路中選取差異最顯著的前10 個基因進行定量表達驗證。對黑化半滑舌鰨皮膚中差異表達基因進行相對表達水平定量分析發現，黑化半滑舌鰨皮膚差異表達基因的表達趨勢與RNA-seq 定量結果一致，表明轉錄組測序結果可靠。針對6 個測序樣本，對6 個候選體色相關mRNA 定量qPCR 表達檢測結果繪制箱線圖（圖5）。其中5 個基因txndc、alox15b、ptgs2、ptgis、atp1a2a在兩組樣本中均有顯著差異（p＜0.05），且有4 個基因在黑化組中表達均高于對照組，只有atp1a2a表達在黑化組表達較低，而acbd7基因則在兩組中無顯著差異（p＞0.05）。

圖5 6 個候選體色相關mRNA 在半滑舌鰨6 個測序樣品中的qPCR 表達檢測結果箱線圖Fig.5 Box plot of qPCR expression results of 6 candidate body color-related mRNAs in 6 sequencing samples of Cynoglossus semilaevis

4 討論

色素的異常沉著是比目魚養殖過程中的一個主要問題，體色異常顯著降低了比目魚的市場價值[39]。近年來，關于半滑舌鰨體色異常的研究多有文獻報道[40-41]，其主要集中于已知體色調控基因的表達模式、著色劑飼料增色效果、理化因子對體色的影響等[42-43]，而從組學角度挖掘半滑舌鰨體色調控新基因的研究則比較缺乏。本研究運用二代測序技術對黑化半滑舌鰨的皮膚轉錄組進行測序分析，得到了半滑舌鰨無眼側黑化皮膚組織和非色素沉著的正常皮膚組織的轉錄組數據，從6 條代謝通路（花生四烯酸代謝通路、甲狀腺激素合成通路、光轉導通路、黑色素生成通路、視黃醇代謝通路和過氧化物酶體增殖物激活受體信號通路）中篩選出半滑舌鰨體色異常相關的差異基因，為今后半滑舌鰨進行體色異常排查提供了參考。選取其中差異最顯著的10 個基因進行驗證后發現有5 個差異基因與測序結果趨勢相一致，分別為txndc、alox15b、ptgs2、ptgis、atp1a2a，且前4 個基因均在黑化組中上調，只有atp1a2a在黑化組中下調表達，暗示了這5 個差異基因可能參與了無眼側皮膚黑化的調控。

TXNDC 是含硫氧還原蛋白，此基因編碼內質網蛋白質的二硫鍵異構酶家族成員，該酶催化蛋白質折疊和巰基-二硫鍵互換反應。硫氧還原蛋白5（TXNDCS），是近年發現的PDI 家族成員之一，具有抗凋亡、抗氧化損傷和促進細胞增殖、促進血管形成、參與細胞炎癥、能量代謝等多種生物功能[44-45]。ALOX15B 是花生四烯酸15 脂氧合酶，編碼該酶的基因編碼參與脂肪酸氫過氧化物生產的結構相關的非血紅素鐵雙加氧酶的脂氧合酶家族成員，其所編碼的蛋白質僅將花生四烯酸轉化為15S-氫過氧二十碳四烯酸[46-47]。PTGS2 是前列腺素內過氧化物合酶2，又稱環加氧酶2，是花生四烯酸（AA）轉化為前列腺素（PGs）的關鍵酶。在前列腺素合成初期起催化作用，可將AA 代謝成各種PGs 類產物[48]，這些PGs 參與維持機體內多種生理過程，是胚胎移植初期時胚胎和子宮信息交流的重要介質之一，在調節發情周期、妊娠維持和分娩中起著關鍵作用[49-50]。PTGIS 是前列環素合酶，是細胞色素P450（CYP450）超家族的8 族（CYP8）成員，也是AA 的主要產物[51]。研究發現，PTGIS 在許多生理和病理過程中起著重要作用，該基因編碼前列腺素12 合酶，并催化前列腺素I2（PGI2）的合成[52]。atp1a2是ATP 酶A2 亞單位，該基因編碼的蛋白質屬于P 型陽離子，轉運ATP 酶家族和Na+/K+-ATP 酶亞家族[53]。

這5 個差異基因中有3 個基因均與AA 有一定的相關性，即ALOX15B 是AA15 脂氧合酶，PTGS2 是AA 轉化為PGs 的關鍵酶，PTGIS 是AA 的主要產物。AA 是一種不飽和脂肪酸，由磷脂酶 A2（Phospholipase A2，PLA2）催化而來，可通過環氧合酶（Cyclooxygenase，COX）、脂氧合酶（Lipoxygenase，LOX）與CYP450 途徑生成多種代謝產物。在海水魚類中，AA 主要在生長、存活、調節免疫、抗應激及繁育等方面發揮重要作用[54-56]。目前，關于AA 在魚類方面的研究主要集中在飼料中改變AA 含量，分析其對魚類的生長性能、脂質代謝、脂肪酸營養、腸道組織、色素沉積等方面的影響。Sargent 等[55]指出，在比目魚的開口餌料中添加一定量的AA，將對背部皮膚色素沉積產生較大影響。Estévez 等[57]認為，飲食中增加的AA 含量對色素沉積有負面影響。McEvoy 等[58]、Copeman 等[59]、Bell 等[60]的研究均表明，高AA 水平與各種比目魚的色素沉著有關[58-60]。Willey 等[61]用富含AA 的輪蟲(Rotifer)和鹵蟲（Artemia）喂養大西洋牙鲆（Paralichthys dentatus），發現增加AA 飼料和色素沉著發生率存在劑量依賴效應。Villalta 等[62]改變了塞內加爾舌鰨（Cynoglossus senegalensis）日糧中AA 的含量，發現隨著AA 的增多伴隨著PG 含量的升高，并且PG 含量與色素沉積存在著緊密的相關性。Lund 等[63]在歐洲鰨（Solea solea）中發現，體色異常只與AA 的絕對量相關，而與AA 和其他多不飽和脂肪酸（Polyunsaturated Fatty Acid，PUFA）的比例無關。Boglino 等[64]研究發現，PGE2 在塞內加爾舌鰨中也引起了體色異常的結果。以上研究都表明AA 的變化會影響比目魚類體表色素沉著。導致AA 誘導比目魚色素沉著的確切機制尚不清楚，但已存在3 種假說：第一，改變飼料中AA 含量可能導致神經組織膜的次優生化組成，而神經組織參與控制變態、黑色素合成和色素團分化[57]；第二，AA 衍生的二十碳五烯酸的過量產生會導致魚類經歷非生物化學誘導的應激，從而導致色素失調[55]；第三，酪氨酸酶的合成和降解是調節色素沉著所必需的，飲食中n-6多不飽和脂肪酸可以改變酪氨酸酶的合成和降解[65]。有研究將正常體色個體與由AA 修飾引起的變色日本牙鲆進行對比研究，發現黑素皮質激素（a-黑素細胞刺激激素a-MSH 和腎上腺皮質激素ACTH）在色素形成過程中起重要作用[66]。而MSH 的作用主要為激活酪氨酸酶，并促進酪氨酸酶合成，從而促進黑色素合成，使皮膚顏色加深。

本研究篩選出半滑舌鰨體色異常相關的5 個差異基因，其中alox15b、ptgs2可以將AA 轉變或者使其代謝為其他物質，從而改變半滑舌鰨中AA 的含量。而ptgis不僅是AA 的主要產物，還是CYP450 中的成員，飼料中營養素對魚類體色的影響主要就是通過魚體色素細胞或色素體起作用的。因此，本文推測alox15b、ptgs2、ptgis基因是通過影響AA 進而影響了半滑舌鰨無眼側體色的色素沉積。本研究為半滑舌鰨無眼側體色異常的進一步研究提供了分子線索，關于AA影響魚類體色異常的分子機制仍有待進一步研究。