生長異質大西洋鮭轉錄組學分析

章 昭，聶苗苗，高 強，劉 丹，王振吉，祁得林*

(1．省部共建三江源生態(tài)與高原農牧業(yè)國家重點實驗室，青海大學，青海 西寧 810016；2．青海大學農牧學院，青海 西寧 810016；3．青海省漁業(yè)環(huán)境監(jiān)測站(青海省高原水生生物及生態(tài)環(huán)境重點實驗室),青海 西寧 810012)

大西洋鮭(Salmosalar), 俗稱三文魚，主要分布于大西洋北部及其附近河流流域，隸屬于硬骨魚綱(Osteichthyes)、輻鰭魚亞綱(Actinopterygii)、鮭形目(Salmoniformes)、鮭科(Salmonidae)，因其營養(yǎng)豐富，經濟價值高，成為人們最喜愛的高檔水產品。根據FAO統(tǒng)計數據顯示，2017年至2018年全球養(yǎng)殖大西洋鮭總產量達244.5萬t，預計2020年達到260萬t[1]。但由于大西洋鮭在中國沒有自然分布，且對生長環(huán)境要求較高，人工育苗和養(yǎng)殖難度較大，所以我國利用工廠化養(yǎng)殖大西洋鮭的技術尚處于起步階段。2016年4月，青海省首次從冰島引進淡水大西洋鮭發(fā)眼卵24 000粒進行孵化，經過2年的工廠化循環(huán)水養(yǎng)殖實驗發(fā)現，同一批次大西洋鮭最大個體體重達到4.6 kg，而最小個體體重不足0.5 kg，養(yǎng)殖群體表現出顯著地生長異質性。目前，對魚類個體生長異質相關的研究多集中于遺傳因素[2]、養(yǎng)殖密度[3]、社會等級因素[4]等方面。但對相同飼養(yǎng)管理條件下，產生生長異質的分子機制尚不明確。

基于此，文中以大西洋鮭為研究對象，利用轉錄組測序技術獲得生長異質大西洋鮭的基因表達差異情況，并對差異基因進行GO功能富集、KEGG通路富集分析，期望能從轉錄組學研究中解析大西洋鮭生長異質的原因，從而為大西洋鮭的遺傳繁育提供參考依據。

1 材料和方法

1.1 試驗材料

本研究以同一養(yǎng)殖環(huán)境下的24月齡雄性大西洋鮭作為研究對象，隨機選擇生長速度較慢，平均體重為0.51 kg的3尾魚，作為慢速生長組；隨機選擇生長速度較快，平均體重為4.20 kg的3尾魚，作為快速生長組(表1)。收集樣品時，用MS-222麻醉大西洋鮭，解剖并分離其肌肉組織，保存于液氮中備用。

表1大西洋鮭測序樣品采集表Tab.1 Collection of samples of Salmo salar sequencing

1.2 試驗方法

1.2.1 總RNA提取及質量檢測

取100 mg肌肉組織樣品，采用 Ambion MagmaxTM-96 total RNA isolation kit (Ambion，美國)試劑盒，按照說明提取總RNA，用1%瓊脂糖凝膠檢測污染情況。用 Nano Photometer?(IMPLEN，USA)分光光度計檢測RNA純度，用 Qubit?RNA Assay Kit(Thermo Fisher，中國)試劑盒檢測RNA濃度，最后使用安捷倫生物分析儀2100系統(tǒng)的 RNA Nano 6000檢測試劑盒檢測RNA完整性。

1.2.2 文庫構建及測序

從6個組織樣品中分別取3 μg RNA，使用Illumina公司的NEBNext UltraTM RNA試劑盒(NEB，USA)構建大西洋鮭測序文庫。大西洋鮭轉錄組文庫的測序由北京百邁客生物科技有限公司通過Illumina HiSeq 2500高通量測序平臺完成。

1.2.3 原始數據質量控制

基于邊合成邊測序(Sequencing by Synthesis，SBS)技術，完成文庫測序。對測序所得的原始數據進行初步處理，去除包含接頭、ploy-N含量大于5%和低質量的數據，從而得到高質量數據(clean reads)并保存為 Fastq 格式[5]。利用Trinity軟件[6]進行轉錄本組裝，最終獲得 Unigenes。

1.2.4 轉錄本注釋

將篩選出的 Unigenes 注釋到以下數據庫中，包括Cluster of Orthologous Groups of proteins (COG)數據庫、Gene Ontology(GO)數據庫、Kyoto Encyclopedia of Genes and Genomes (KEGG)數據庫、clusters of orthologous groups for eukaryotic complete genomes(KOG)數據庫、Pfam 數據庫、Swiss-prot 數據庫、evolutionary genealogy of genes：non-supervised orthologous groups(eggNOG)數據庫和non-redundant protein sequences(NR)數據庫。

1.2.5 Unigenes差異表達分析

以 Bowtie[7]軟件為基礎，將轉錄組測序得到的 clean reads 比對到參考基因組上(https:∥www.ncbi.nlm.nih.gov/genome/term=Salmo)，隨后用 DESeq R包[8](依據負二項分布原理)計算基因和轉錄本的表達水平。比較慢速生長組和快速生長組各三個重復樣品并計算其基因表達平均值。采用 Benjamini-Hochberg方法對P值(P-value)進行校正，并采用錯誤發(fā)現率(False Discovery Rate，FDR)作為差異表達基因篩選的關鍵指標，而FDR是通過對差異顯著性P值校正得到的，最終將差異倍數(Fold Change)≥2且FDR<0.05作為篩選標準，并對篩選得到差異表達基因進行GO功能富集和KEGG通路富集等分析。

2 結果與分析

2.1 大西洋鮭轉錄組測序質量評估

由表2可知，通過 Illumina HiSeq 2500 高通量測序平臺測序，6個樣品共產生81.27 Gb Clean Data，每個樣品的堿基量均達到12.61 Gb，去除帶接頭、低質量和短序列后，共獲得271 841 994大西洋鮭高質量序列。樣品測序的堿基錯誤率均低于0.05%，GC含量均在50.7%左右，Q30堿基百分比均大于92.55%，說明測序質量較高。

表2大西洋鮭樣品轉錄組測序數據統(tǒng)計Tab.2 Data statistics of transcriptome sequencing of Salmo salar samples

2.2 基因功能注釋

為獲得全面的基因功能信息，將Unigenes注釋到8個數據庫中，其中注釋基因最多的數據庫是NR數據庫，其次是eggNOG數據庫，最少的是COG數據庫(表3)。

表3大西洋鮭基因功能注釋結果統(tǒng)計Tab.3 Statistical results of gene functional annotation of Salmo salar

2.3 基因表達水平分析

使用Bowtie軟件將每一個樣品的clean reads比對到參考基因組，比對統(tǒng)計結果見表4。從比對結果來看，各樣品的clean reads與參考基因組的比對效率為82.15%～87.04%。

表4 clean reads與參考序列比對統(tǒng)計Tab.4 Clean reads and comparison of reference sequence

2.4 差異表達基因分析

通過對慢速生長組和快速生長組6個樣品進行樣品間相關性分析發(fā)現，Xiaoxiong3-X3和Daxiong1-S1兩個樣品組內相關性較低，因此后續(xù)對兩組樣品進行差異表達基因分析時剔除這兩個樣品。最終篩選出948個差異表達基因，其中包括363個上調基因，585個下調基因(快速生長組中的表達水平高于慢速生長組中的表達水平稱為上調基因，反之稱為下調基因),如圖1所示。

圖1 差異表達基因火山圖Fig.1 Volcano plot of differentially expressed genes

2.5 差異表達基因的GO富集和KEGG通路富集分析

GO功能富集分析顯示，共有543個差異表達基因得到歸類注釋。在生物過程(biological process)、細胞組分(cellular component)、分子功能(molecular function)3大類中大部分Term都富集到了基因(圖2)。KEGG通路富集分類中共有229個差異表達基因參與到了6大類50小類114個代謝通路中， KEGG通路富集分析發(fā)現，6個分支所富集的基因占比例不同：環(huán)境信息處理(environmental information processing，56.33%)、代謝(metabolism，51.53%)、細胞過程(cellular processes，33.18%)、生物體系統(tǒng)(organismal system，27.51%)、人類疾病(human diseases，6.11%)和遺傳信息處理(genetic information processing，3.93%)。注釋基因最多的3個通路依次是神經活性受體-配體交互作用信號通路(neuroactive ligand-receptor interaction pathway)、鈣信號通路(calcium signaling pathway)和緊密連接信號通路(tight junction pathway)(圖3)。

圖2 差異基因的GO功能分類Fig.2 Classification of GO function of differentially expressed genes

圖3 KEGG注釋和分類Fig.3 Annotation and classification of KEGG

3 討論與結論

(1)本研究采用 Illumina HiSeq 2500高通量測序平臺對生長異質大西洋鮭轉錄組文庫進行高通量測序分析，經過測序質量控制，共得到81.27 Gb Clean Data，各樣品Q30堿基百分比均大于92.55%，高于多個已發(fā)表的魚類轉錄組組裝結果，包括斑馬魚(Daniorerio)(91.93%)[9]、金錢魚(Scatophagusargus)(91.87%)[10]等，通過與參考基因組比對，各樣品的clean reads與參考基因組的比對效率為82.15%～87.04%，說明本試驗測序質量較好，可進行后續(xù)分析。

(2)肌球蛋白(myosin)是一種普遍存在的真核運動蛋白，也是肌肉的主要組成蛋白質，它與肌動蛋白相互作用，參與從細胞分裂到肌肉收縮等各種細胞運動的過程[11]。本研究發(fā)現，參與緊密連接代謝通路調控的肌球蛋白重鏈(myosin heavy chain，MyHC)基因出現顯著上調。相關的研究發(fā)現,MyHC基因的表達存在于整個魚類的肌肉發(fā)育過程中[12]。魚類的生長是蛋白質和脂質沉積的綜合結果[13]，Hevr?y發(fā)現大西洋鮭的生長速度與MyHC基因的表達量呈正相關[14]，在對翹嘴鱖(Sinipercachuatsi)早期生長階段的研究也發(fā)現MyHC在快速生長組中的表達量更高[15]，這一結果與本研究的結果相一致，說明魚類的生長速度可能與MyHC相關。

蘇氨酸作為一種必須氨基酸，在幼齡動物的生長發(fā)育中具有重要的調控作用，包括參與蛋白質的合成，促進生長等[16]。在美國紅魚(Sciaenopsocellatus)和歐洲黑鱸(Dicentrarchuslabrax)的研究中發(fā)現，飼料中蘇氨酸缺乏會導致生長下降和飼料系數升高[17-18]。竇秀麗等[19]對鱸魚的研究中發(fā)現當飼料中蘇氨酸含量為1.83%～1.87%時，鱸魚的生長性能和蛋白利用效率達到最佳水平。本研究中參與蘇氨酸代謝調控的基因在快速生長組中表達量顯著高于慢速生長組，說明在快速生長組中，較高的蘇氨酸代謝水平是大西洋鮭生長較快的原因之一。

過氧化物酶體增殖物激活受體(peroxisome proliferators-activated receptors，PPARs)通路在脂質代謝、葡萄糖代謝和調節(jié)能量平衡過程中發(fā)揮重要作用[20]。PPARs家族包括3種同工型：PPARα，PPARβ/δ和PPARγ[21]。其中PPARα在肝臟、心臟、骨骼肌等代謝活躍的組織中高表達，主要作用于脂肪酸代謝[22]。在斑馬魚(Daniorerio)和紅鯛魚(Pagrusmajor)的研究中發(fā)現PPARs的激活是導致魚類脂肪形成和肥胖的重要因素[23-24]。本研究中有9個基因富集到PPARs通路中，其中包括參與血脂代謝和利用的載脂蛋白AI(apolipoproteinAI，apoAI)和負責氨基酸、核苷酸、糖等吸收和運輸功能的溶質載體家族27(solute carrier family 27)，并且這兩個基因表達水平均顯著上調，這些基因主要參與PPARs轉錄的調控，當PPARs被激活后機體脂肪酸濃度增加[25]，蛋白脂肪酶的合成增加，進而調控脂肪細胞的分化[26]，這可能也是大西洋鮭個體生長速度較快的重要因素。

鈣/鈣調蛋白依賴性蛋白激酶Ⅱ(calcium/calmodulin-dependent protein kinase Ⅱ，CaMKⅡ )是一種多功能蛋白激酶，可磷酸化40多種蛋白質[27]。CaMKⅡ不僅參與Ca2+的代謝，而且在維持細胞能量代謝方面也發(fā)揮重要作用[28]。本研究中，CaMKⅡ表達量出現顯著下調，同時，與其調控相關的谷氨酸受體(glutamate receptor ionotropic，NMDA1和NMDA2D)基因也出現顯著的下調?？焖偕L組中CaMKⅡ和NMDA的表達量下調說明在慢速生長組中該基因的表達量相對較高，NMDA對Ca2+具有高度通透性，并可使胞內Ca2+濃度升高，從而激活CaMKⅡ[29]。也有研究表明細胞內游離Ca2+濃度升高可能引起肌肉收縮、激素分泌等生理過程的改變[30]。而維持機體細胞內外Ca2+濃度平衡是保證機體正常生長發(fā)育的首要條件，劉亞靜等[31]在不同Ca2+濃度對青海湖裸鯉(Gymnocyprisprzewalskii)幼魚生長的研究中發(fā)現高Ca2+和低Ca2+均不同程度的減緩了青海湖裸鯉的生長速度。Laining 等[32]在對紅鰭東方鲀(Takifugurubripes)幼魚喂食高鈣餌料后發(fā)現該魚生長速度明顯減慢。因此，慢速生長組中CaMKⅡ 和NMDA相對較高的表達水平引起胞內Ca2+濃度升高可能對一部分大西洋鮭的生長產生了不利的影響。

魚類生長過程受到環(huán)境、溫度、采食等多種復雜因素的影響，因此準確的闡明影響其生長異質的確切機制是十分困難的，本試驗通過對生長異質大西洋鮭的轉錄組測序分析獲得了大西洋鮭的轉錄組數據、基因注釋結果，為后續(xù)探究大西洋鮭生長異質的主要原因提供數據支撐。同時，通過差異表達基因富集分析，對與本研究相關的鈣信號通路、緊密連接通路、PPARs等通路調控機制進行分析，揭示了與生長異質相關的生物學過程和分子機制，為大西洋鮭人工養(yǎng)殖繁育提供科學依據。