章 昭，聶苗苗，高 強(qiáng)，劉 丹，王振吉，祁得林*

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

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

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

1 材料和方法

1.1 試驗(yàn)材料

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

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

1.2 試驗(yàn)方法

1.2.1 總RNA提取及質(zhì)量檢測

取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 文庫構(gòu)建及測序

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

1.2.3 原始數(shù)據(jù)質(zhì)量控制

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

1.2.4 轉(zhuǎn)錄本注釋

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

1.2.5 Unigenes差異表達(dá)分析

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

2 結(jié)果與分析

2.1 大西洋鮭轉(zhuǎn)錄組測序質(zhì)量評估

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

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

2.2 基因功能注釋

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

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

2.3 基因表達(dá)水平分析

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

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

2.4 差異表達(dá)基因分析

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

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

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

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

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

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

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

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

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

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