不同生長速率花鱸肌肉轉(zhuǎn)錄組分析及生長相關(guān)基因篩選

范嗣剛，楊文燕，黃 皓,3，王鵬飛，趙 超，閆路路，張 博，邱麗華,3

（1.中國水產(chǎn)科學(xué)研究院南海水產(chǎn)研究所/廣東省漁業(yè)生態(tài)環(huán)境重點實驗室，廣東 廣州 510300；2.廣東省農(nóng)業(yè)技術(shù)推廣中心，廣東 廣州 510520；3.上海海洋大學(xué)水產(chǎn)與生命學(xué)院，上海 201306）

除環(huán)境外，魚類生長還受許多基因的影響和控制。目前已挖掘出一些生長相關(guān)基因。生長激素（Growth hormone，GH）-胰島素樣生長因子（Insulinlike growth factor，IGF）生長軸是魚類生長發(fā)育的主要作用通道[1]。增加魚類體內(nèi)的GH 含量可促進(jìn)生長[2]。IGFI和IGFII可促進(jìn)葡萄糖進(jìn)入肌肉，加速肌細(xì)胞的分裂和分化[3-4]。肌肉生長抑制素（Myo‐statin，MSTN）是轉(zhuǎn)化生長因子-β (Transforming growth factor-beta，TGF-β)超家族的成員，在動物肌肉生長和發(fā)育過程中起負(fù)調(diào)控作用[5]。骨骼肌中的α-肌動蛋白基因在大、小規(guī)格虹鱒（Oncorhynchus mykiss）中差異性表達(dá)[6]，在快速生長鲇（Ictalurus punctatus）中表達(dá)量顯著升高[7]。成纖維細(xì)胞生長因子（Fibroblast growth factor，F(xiàn)GF）可調(diào)控細(xì)胞增殖、遷移和分化[8]。轉(zhuǎn)錄組測序技術(shù)（RNA-sequenc‐ing，RNA-seq）可準(zhǔn)確、有效、快速地獲得組織轉(zhuǎn)錄組序列及表達(dá)情況，已廣泛用于水產(chǎn)動物生長基因的篩選和鑒定等研究。孫雪等[9]對草魚（Ctenopharyngodon idella）快長組和慢長組個體的肌肉進(jìn)行RNA-seq，獲得552 個差異表達(dá)基因（Differentially expressed genes，DEGs）。Zhang 等對生長快和生長慢的青魚（Mylopharyngodon piceus）進(jìn)行轉(zhuǎn)錄組分析，獲得1 913 個DEGs，在生長和發(fā)育相關(guān)的代謝途徑中富集[10]。Luo 等[11]用轉(zhuǎn)錄組測序技術(shù)研究生長快和生長慢的鳙（Hypophthalmichthys nobilis），鑒定出15 個關(guān)鍵信號途徑和20 個DEGs。Cleveland等[12]通過比較虹鱒（O.mykiss）快、慢生長組的轉(zhuǎn)錄組數(shù)據(jù)，獲得36 個DEGs。王玉梅等[13]研究2 種不同生長速率的福瑞鯉（Cyprinus carpio）轉(zhuǎn)錄組，獲得749 個DEGs，鑒定出肌紅蛋白（Myoglobin,mb）、肌球蛋白輕鏈2（Myosin light chain 2,myl2）等與肌肉生長相關(guān)的關(guān)鍵基因。

花鱸（Lateolabrax maculatus）隸屬于鱸形目（Perciformes）鮨 科（Serranidae）花鱸屬（Lateolabrax），俗稱鱸魚、七星鱸、花寨等，廣泛分布于我國渤海至南海的近岸水域以及通海江河水域，是一種廣溫廣鹽性經(jīng)濟(jì)魚類，有生長較快、適應(yīng)性廣、味道鮮美等特點[14-15]。目前，已有花鱸生長的研究。Liu 等[16]構(gòu)建了花鱸首個高密度連鎖圖譜，篩選出24個與生長相關(guān)的數(shù)量性狀位點（QTL），獲得30個生長有關(guān)候選基因。胰島素樣生長因子結(jié)合蛋白基因（IGFBP）、IGF、生長激素受體基因（GHR）、甲狀腺激素受體（Thyroid hormone receptor，TR）和瘦素A（LeptinA，LepA）等與花鱸生長有關(guān)的基因序列及其組織表達(dá)也有研究[17-18]。但尚未見對花鱸生長方面的轉(zhuǎn)錄組研究報道。魚類肌肉的生長發(fā)育直接決定魚類的總質(zhì)量。本研究擬對同一養(yǎng)殖環(huán)境下生長快與慢的花鱸肌肉組織進(jìn)行轉(zhuǎn)錄組測序，研究二者在轉(zhuǎn)錄組水平上的差異，篩選出與生長相關(guān)的差異表達(dá)基因和信號通路，為解析花鱸生長的分子機(jī)制奠定基礎(chǔ)。

1 材料與方法

1.1 材料

花鱸同批次繁殖苗種（體長2～3 cm）1 萬尾，2020 年3 月取自寧德市富發(fā)水產(chǎn)有限公司，在南海水產(chǎn)研究所珠海試驗基地池塘中養(yǎng)殖9個月。取花鱸500 尾，其中461 尾生長快速，達(dá)到商品規(guī)格（500～1 000 g），39 尾生長緩慢，個體短?。?0～150 g）。隨機(jī)選取4 尾生長快[體長(34.22±2.53)cm]和4尾生長慢[體長(17.46±1.78)cm]的花鱸個體，取肌肉組織液氮保存?zhèn)溆谩?/p>

1.2 RNA提取、cDNA文庫構(gòu)建及測序

用Tirzol 法（Thermofisher scientific，美國）提取花鱸肌肉組織的RNA，用DNase（Takara，日本）消化DNA。用帶有Oligo (dT) 的磁珠（杭州尚炫生物科技有限公司）富集mRNA，加入打斷試劑將mRNA打斷成短片段。以打斷后的mRNA 為模板，用六堿基隨機(jī)引物合成一鏈cDNA，配制二鏈合成反應(yīng)體系合成二鏈cDNA，并用試劑盒純化雙鏈cDNA。對純化的雙鏈cDNA 進(jìn)行末端修復(fù)、加A 尾并連接測序接頭、片段大小選擇、PCR 擴(kuò)增。構(gòu)建的文庫用2100 生物分析儀（Agilent，美國）質(zhì)檢合格后，用Illumina HiSeq X 測序儀（Illumina，美國）測序，產(chǎn)生150 bp的雙端數(shù)據(jù)。文庫構(gòu)建和轉(zhuǎn)錄組測序均由上海生工生物工程股份有限公司完成。

1.3 測序數(shù)據(jù)組裝

用Trimmomatic 軟件[19]分析原始數(shù)據(jù)（Raw reads），去除reads的接頭序列、低質(zhì)量reads（質(zhì)量值Qphred＜20）、N 堿基比例大于5%的reads 等，得到純凈讀數(shù)（Clean reads）。

用Trinity[20]將干凈讀長從頭組裝成轉(zhuǎn)錄本（Transcript），參數(shù)為min_kmer_cov 2，其余參數(shù)為默認(rèn)值。轉(zhuǎn)錄本去冗余后，取每個轉(zhuǎn)錄本聚類中最長的轉(zhuǎn)錄本作為unigenes，以此作為后續(xù)分析的參考序列。

1.4 基因功能注釋

將獲得的unigenes 在NT (NCBI nucleotide sequences)、Nr(NCBI non-redundant protein sequences)、KOG (euKaryotic Ortholog Groups)和Swiss-Prot 等數(shù)據(jù)庫進(jìn)行比對和注釋，e值均小于1e-5。

用Blast 2 GO 軟件（version 2.5）對unigenes 進(jìn)行GO 功能注釋（e-value ＜1e-6），用KAAS[21]軟件對unigenes進(jìn)行KEGG注釋（e-value ＜1e-10）。

1.5 差異表達(dá)基因的鑒定與分析

用Salmon[22]計算基因的表達(dá)量。用RSEM 軟件將每個樣品的clean reads 比對到拼接轉(zhuǎn)錄本中，獲得讀長數(shù)量（read count）。用TMM 軟件對read count 數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，用DESeq2 1.12.4 軟件包分析生長快和生長慢花鱸的差異表達(dá)基因。篩選條件設(shè)為：q ＜0.05 且|log2FC|＞2，F(xiàn)C 為差異倍數(shù)（Fold change），q為多重假設(shè)檢驗校正后的P值。

用topGO（version 2.24）軟件包對DEGs 進(jìn)行GO 富集分析，用clusterProfiler（version 3.0.5）軟件包進(jìn)行KEGG通路富集分析。

1.6 差異表達(dá)基因的驗證

隨機(jī)選取轉(zhuǎn)錄組中的13 個差異表達(dá)基因進(jìn)行qRT-PCR驗證（表1）。用與轉(zhuǎn)錄組測序同批次的RNA進(jìn)行反轉(zhuǎn)錄，獲得相應(yīng)的cDNA。β-actin為內(nèi)參基因。用SYBR Green Pro Taq HS預(yù)混型qPCR試劑盒（艾科瑞，中國）在LightCycler 480熒光定量PCR儀（Roche，美國）上進(jìn)行qRT-PCR。反應(yīng)體系為12.5 μL，包含cDNA（40 ng/μL）5.25 μL，上下游引物（10 μmol/L）各0.5 μL，2X SYBR Green Pro Taq HS Premix 6.25 μL。反應(yīng)條件：95 ℃30 s；94 ℃5 s，60 ℃30 s，共40 個循環(huán)。反應(yīng)結(jié)束后，對PCR 產(chǎn)物進(jìn)行溶解曲線檢測，以驗證擴(kuò)增產(chǎn)物的特異性。每次試驗做3 個樣品，每個樣品重復(fù)檢測3 次。用2-ΔΔCt法計算基因的相對表達(dá)量。

表1 用于實時定量PCR驗證的引物Table 1 Primers used in Real-time PCR confirmation

2 結(jié)果與分析

2.1 花鱸肌肉轉(zhuǎn)錄組測序與組裝結(jié)果

對生長快和生長慢花鱸的肌肉組織分別獲得原始序列44 434 288、48 256 446 個，Clean reads 42 928 714、46 975 694 個，Q30值分別為91.65%、91.84%（表2）。測序結(jié)果已提交到NCBI SRA 數(shù)據(jù)庫（PRJNA811720）。

表2 花鱸肌肉轉(zhuǎn)錄組數(shù)據(jù)的基本信息Table 2 Summary statistics of Lateolabrax maculatus muscle transcriptome sequencing

兩個Clean reads拼接后獲得137 960個unigenes；N50為637 bp，N90234 bp；總長70 483 974 bp，最大長度19 775 bp，最小長度201 bp，平均長度510.9 bp；長度≥500 bp unigenes有33 908個，長度≥1 000 bp的有13 373個（圖1）。

圖1 花鱸unigene長度分布Fig.1 Unigene length distribution in transcriptome of Lateolabrax maculatus

2.2 基因功能注釋

unigenes 注釋結(jié)果如表3 所示。有68 334 個unigenes（49.53%）在Nr 中注釋，68 305 個unigenes（49.51%）在Swissprot 中注釋，在GO 數(shù)據(jù)庫注釋到70 946 個unigenes（51.43%）。有105 310 個unigenes在至少一個數(shù)據(jù)庫中獲得注釋信息，有5 215 個unigenes在所有數(shù)據(jù)庫中均有注釋。

表3 花鱸肌肉轉(zhuǎn)錄組測序結(jié)果統(tǒng)計Table 3 Statistics of Lateolabrax maculatus muscle transcriptome sequencing

2.3 差異表達(dá)基因鑒定與分析

根據(jù)q ＜0.05 且|log2FC|＞2 的標(biāo)準(zhǔn)，快速生長組與慢速生長組比較，有10 552 個差異表達(dá)基因（圖2）。其中，2 064 個unigenes 表達(dá)量上調(diào)（1 661個unigenes 在Nr 中成功注釋），8 488 個unigenes 表達(dá)量下調(diào)（3 793 個unigenes 在Nr 中成功注釋）。在這些DEG中，有一些與肌肉生長和發(fā)育相關(guān)的基因被篩選出來，如胰島素樣生長因子結(jié)合蛋白1（Insu‐lin-like growth factor-binding protein 1，IGFBP1）、肌球蛋白重鏈（Myosin heavy chain，MHC）、肌球蛋白輕鏈（Myosin light chain，MLC）、成纖維細(xì)胞生長因子結(jié)合蛋白3（Fibroblast growth factor-binding pro‐tein 3，F(xiàn)GFBP3）、多表皮生長因子樣結(jié)構(gòu)域蛋白9（Multiple epidermal growth factor-like domain pro‐tein 9，MEGFL 9）、成纖維細(xì)胞生長因子（Fibroblast growth factor，F(xiàn)GF）、生長激素受體1（Growth hor‐mone receptor 1，GHR1）、肌生成抑制蛋白（Myo‐statin，MSTN）、生肌因子（Myogenic factor，MYF）、肌凝蛋白結(jié)合蛋白（Myosin binding protein，MYBP）、生長抑制蛋白（Inhibitor of growth protein，ING）等。

圖2 慢速、快速生長組花鱸肌肉的差異表達(dá)基因M-A圖Fig.2 M-A plot for DEGs comparison between the slowand fast-growing group of Lateolabrax maculatus

2.4 差異表達(dá)基因功能注釋

對10 552個DEGs進(jìn)行GO富集分析，在生物過程富集的功能分支（Term）最多，有9 493 個，其次在分子功能中，富集2 345個，在細(xì)胞成分中富集1 358個。在每個類別（Ontology）中，顯著富集的前10 個功能分支如表4所示。在生物過程中顯著富集的前10 個功能分支中，多為前體、多肽等物質(zhì)的代謝。在細(xì)胞成分中，有6 個與核糖體有關(guān)，2 個與肌肉有關(guān)。在分子功能排名前5 的是核糖體的結(jié)構(gòu)成分（Structural constituent of ribosome）、結(jié)構(gòu)分子活性（Structural molecule activity）、rRNA 結(jié) 合（rRNA binding）、RNA結(jié)合（RNA binding）、細(xì)胞色素c氧化酶活性（Cytochrome-c oxidase activity）。

表4 不同規(guī)格花鱸差異表達(dá)基因在GO功能三大類別中富集前10的條目分析Table 4 Analysis of the first 10 terms about DEGs for large and small size of Lateolabrax maculatus in the three majorcategories of GO function

在KEGG 富集分析中，共富集到255 個通路。上調(diào)的DEGs 富集最多的前5 個KEGG 通路為催產(chǎn)素信號通路、胰高糖素信號通路、脂肪酸降解、腺苷酸活化蛋白激酶（AMPK）信號通路、碳代謝，下調(diào)的DEGs富集最多的前5個KEGG通路為核糖體、氧化磷酸化、心臟肌肉收縮、抗原加工和提呈、氨基酸生物合成（圖3）。

圖3 差異表達(dá)基因顯著富集的前10個KEGG通路Fig.3 Top 10 enriched KEGG pathways in DEGs

2.5 差異表達(dá)基因的驗證

隨機(jī)挑選13個DEGs中，有8個上調(diào)基因和5個下調(diào)基因。圖4 表明，所有的差異表達(dá)基因的表達(dá)模式與測序結(jié)果一致，表明測序結(jié)果的可信度高。

圖4 差異表達(dá)基因的實時定量PCR驗證Fig.4 qRT-PCR validation analysis of DEGs

3 討論

本研究中，同一批花鱸苗種（體長2～3 cm）放入池塘中養(yǎng)殖，在養(yǎng)殖水體、飼料等條件一致的情況下，9 個月后獲得兩種規(guī)格不同的花鱸個體，可排除環(huán)境因素對生長的影響，確定某些基因和信號通路對生長的調(diào)控作用。

魚類的肌肉生長包括肌纖維增生和肥大，受多種基因和信號通路的調(diào)控[23]。GH/IGF 軸在調(diào)節(jié)魚類生長方面起重要作用[1]。本研究篩選出該軸上的GHR1、IGFBP1、IGFBP7和IGFBP4等 DEGs。GHR1是GH的受體，能促進(jìn)細(xì)胞生長、增殖，減少蛋白降解[24]。魚類的GHRs 是單一跨膜蛋白，包括胞外結(jié)構(gòu)（包含6 或7 半胱氨酸殘基）、跨膜結(jié)構(gòu)和胞內(nèi)結(jié)構(gòu)[25]。虹鱒（O.mykiss）GHR1 主要通過信號轉(zhuǎn)導(dǎo)及轉(zhuǎn)錄激活因子5（STAT5）起作用[26]?；|GH序列已被克隆出來[17]，但未見花鱸GHR的研究報道。本研究發(fā)現(xiàn)，生長快的花鱸GHR1表達(dá)量顯著上調(diào)，說明GHR1能促進(jìn)花鱸生長。IGFBP1受胰島素的負(fù)調(diào)控，與魚類的生長負(fù)相關(guān)[27]。IGFBP4在體內(nèi)是一種抑制劑，通過抑制IGF的活性，阻止受體與IGF 結(jié)合[28-29]。IGFBP7，又 名IGFBP 相關(guān)蛋白1（IGFBP-related protein 1，IGFBP-rP1），與IGF 結(jié)合能力弱，與胰島素結(jié)合能力強(qiáng)[30-31]。在生長快和生長慢的大口黑鱸（Micropterus salmoides）轉(zhuǎn)錄組比較研究中，IGFBP1也被篩選出來[32]。肌球蛋白重鏈和肌球蛋白輕鏈均為肌球蛋白的重要組成部分，參與了肌肉的生長發(fā)育過程[33]。生長抑制蛋白家族包括ING1、ING2、ING3、ING4和ING5等5 個成員[34]。本研究檢測出ING3、ING4和ING5等DEGs。生長抑制蛋白能抑制細(xì)胞增殖、遷移和擴(kuò)散[35]。Liu等[16]基于花鱸高密度遺傳連鎖圖譜，發(fā)現(xiàn)FGFR4、FGF12a、FGF18 等3 個FGF 家族成員與花鱸生長有關(guān)。本研究DEGs 中有4 個FGF 家族成員。FGFBP3 是一種分泌分子伴侶，能影響碳水化合物和脂質(zhì)代謝，在老鼠中敲除FGFBP3后，能改變脂質(zhì)代謝[36]。FGF 6 的功能主要有肌源干細(xì)胞遷移，肌肉分化、增生[37]。草魚的FGF6 表達(dá)量與肌纖維直徑正相關(guān)[38]。FGFRL1 通過結(jié)合FGF，對細(xì)胞增殖起負(fù)調(diào)控作用[39]，對骨形成和細(xì)胞分化起正調(diào)控作用[40]。在草魚饑餓再投喂實驗中，F(xiàn)GFRL1表達(dá)量顯著上調(diào)[41]。在哺乳動物中，F(xiàn)GF2是一種原型血管生成生長因子，能誘導(dǎo)血管生成和新血管形成[42]。在斑馬魚（Danio rerio）中也發(fā)現(xiàn)FGF2 能促進(jìn)血管生成[43]。

在GO的細(xì)胞組分中，有肌原纖維和肌節(jié)等2個DEGs 顯著富集，KEGG 途徑中有心肌收縮的DEGs顯著富集。在福瑞鯉、青魚和鳙等魚類中也發(fā)現(xiàn)有DEGs 匹配到肌原纖維和心肌收縮等GO 和KEGG富集中[10-11,13]。由此可見，在不同魚類中，肌原纖維和心肌收縮與肌肉生長都緊密相關(guān)。

研究表明，物質(zhì)和能量代謝能影響生物生長[44]。本研究也證實這一點：DEGs的GO 和KEGG 富集結(jié)果多集中在物質(zhì)和能量代謝方面。在GO 富集的生物過程中，排在前10 的是蛋白、核酸等物質(zhì)代謝。在生長快速的草魚和鳙中，富集在代謝過程的DEGs表達(dá)量顯著上調(diào)[11,45]。

本研究部分KEGG 分析結(jié)果與其他魚類的研究結(jié)果有相同之處，如糖酵解/糖異生途徑[32]、脂肪酸降解途徑[11,13,46]、AMPK 信號通路[10]、碳代謝和脂肪酸代謝[10-11]、精氨酸和脯氨酸代謝[10-11,44]、D-谷氨酰胺和D-谷氨酸鹽代謝[47]。這些KEGG 途徑均與肌肉生長發(fā)育有關(guān)，如糖酵解是糖原產(chǎn)生ATP 的過程，糖異生是非糖有機(jī)物轉(zhuǎn)變成葡萄糖的過程。虹鱒再投喂試驗發(fā)現(xiàn)有一些基因與糖酵解/糖異生有關(guān)[48]。AMPK在調(diào)控肌肉生長和再生方面起重要作用，能促進(jìn)肌肉匯總葡萄糖的吸收和脂肪酸的氧化，為肌肉收縮提供能量[49]。攝入精氨酸和脯氨酸等氨基酸能增強(qiáng)肌肉生長，減少過多的脂肪[50]。今后可進(jìn)一步研究這些KEGG 途徑與魚類肌肉生長發(fā)育的關(guān)系。

本研究中，分別有6 個細(xì)胞組分和2 個分子功能與核糖體有關(guān)，KEGG 富集中下調(diào)最顯著的是核糖體，可見核糖體對魚類生長的重要性。事實上，核糖體是細(xì)胞合成蛋白質(zhì)的主要場所，RNA 翻譯到蛋白質(zhì)這一過程即在核糖體中完成。核糖體蛋白和rRNA 組成兩個大小不同的核糖體亞基。研究發(fā)現(xiàn)，核糖體蛋白能控基因轉(zhuǎn)錄、mRNA 翻譯，影響細(xì)胞分化、增殖和凋亡，對肌肉生長有重要作用[51]。