王子予，程 超，朱敏杰，肖 敏，查云飛，簡少卿，趙大顯*

（1.南昌大學(xué) 生命科學(xué)學(xué)院/江西省水產(chǎn)動物資源與利用重點實驗室，江西 南昌 330031；2.江西省吉安市井岡山農(nóng)業(yè)科技園，江西 吉安 343000；3.井岡山沐云特種動物養(yǎng)殖有限公司，江西 吉安 343000）

【研究意義】西伯利亞雜交鱘是以西伯利亞鱘為母本、施氏鱘為父本雜交培育出來的一個養(yǎng)殖品種[1]，它具有生長速度快、抗病能力強、耐運輸?shù)忍攸c，因此在生產(chǎn)實踐中深受廣大養(yǎng)殖戶的歡迎。但隨著集約化養(yǎng)殖的發(fā)展，為了追求更高的利潤，養(yǎng)殖戶往往通過增加養(yǎng)殖密度來提高養(yǎng)殖產(chǎn)量，養(yǎng)殖密度是一個復(fù)雜且關(guān)鍵的環(huán)境影響因子，適當?shù)奶岣唣B(yǎng)殖密度是可以達到增加養(yǎng)殖產(chǎn)量的目的，但盲目增加養(yǎng)殖密度容易導(dǎo)致水質(zhì)惡化以及種內(nèi)對餌料和空間的競爭，進而導(dǎo)致魚類生長水平、存活率以及抗病能力的降低。養(yǎng)殖密度也會造成個體生長的差異，且隨著養(yǎng)殖密度的增加，差異的趨勢越明顯[2]，因此，研究養(yǎng)殖密度對西伯利亞雜交鱘生長性能的影響及其作用的分子機制對該品種的健康養(yǎng)殖顯得尤為重要。【前人研究進展】高通量測序技術(shù)也被稱為下一代測序（NGS）技術(shù)，近些年有一種低成本的NGS 技術(shù)受到追捧，它可以用來分析動態(tài)轉(zhuǎn)錄組，該技術(shù)即為RNA 測序（RNA-Seq）。RNA-Seq在硬骨魚類的轉(zhuǎn)錄組學(xué)研究中能夠更好的幫助了解魚類的許多生物學(xué)過程，如生長發(fā)育、適應(yīng)性進化、免疫反應(yīng)以及應(yīng)激反應(yīng)[3]。李永娟[6]利用RNA-Seq 技術(shù)分析了熱應(yīng)激條件下虹鱒肝臟和頭腎組織轉(zhuǎn)錄組數(shù)據(jù)，獲得了2 種組織共有差異基因36 個，發(fā)現(xiàn)了虹鱒熱適應(yīng)相關(guān)的調(diào)控通路主要為內(nèi)質(zhì)網(wǎng)的蛋白加工過程和免疫相關(guān)的“NOD 樣受體信號通路”，為進一步揭示虹鱒熱應(yīng)激反應(yīng)的分子調(diào)控機制提供了科學(xué)數(shù)據(jù)；張明洋等[7]利用RNA-Seq 對類志賀鄰單胞菌感染雜交鱘后的腸道組織進行了轉(zhuǎn)錄組分析，共得到差異表達基因13 542 個，其中上調(diào)基因9 774 個，下調(diào)基因為3 768 個，并篩選可能影響腸道黏膜免疫的相關(guān)通路和基因，為類志賀鄰單胞菌感染雜交鱘誘導(dǎo)腸道黏膜免疫應(yīng)答反應(yīng)的分子機制奠定了基礎(chǔ)；李營等[8]利用RNA-Seq 技術(shù)對人工養(yǎng)殖2 齡施氏鱘（Acipenser schrenckiiBrandt）精巢與卵巢進行轉(zhuǎn)錄組測序分析，獲得了雌雄性腺中共有19 690 個差異表達基因轉(zhuǎn)錄本，并發(fā)現(xiàn)了4 條與卵巢發(fā)育相關(guān)的通路，為研究鱘性腺分化和發(fā)育機制提供了基礎(chǔ)。因此，RNA-Seq 技術(shù)在冷水性魚類差異基因篩選和候選基因發(fā)掘上發(fā)揮重要功能?！颈狙芯壳腥朦c】本研究結(jié)合生產(chǎn)實踐，應(yīng)用RNA-Seq技術(shù)分析不同養(yǎng)殖密度條件西伯利亞雜交鱘幼魚轉(zhuǎn)錄組信息，篩選獲得與生長應(yīng)激相關(guān)的基因和調(diào)控通路，為井岡山地區(qū)西伯利亞雜交鱘健康養(yǎng)殖制定適宜的養(yǎng)殖密度提供理論依據(jù)?！緮M解決的關(guān)鍵問題】從分子水平闡明養(yǎng)殖密度對西伯利亞雜交鱘幼魚生長性能影響機制，為井岡山地區(qū)鱘魚流水健康養(yǎng)殖提供科學(xué)參考數(shù)據(jù)。

1 材料和方法

1.1 試驗條件

研究對象為西伯利亞雜交鱘，試驗用魚選自同一批次、規(guī)格一致、健康無損。養(yǎng)殖實驗在井岡山沐云特種動物養(yǎng)殖有限公司鱘魚養(yǎng)殖基地進行，實驗池為9 個四邊形室外水泥池，單個水泥池占地面積為28 m2，水深維持在40 cm 左右。流水池主要由進水口、魚池、出水口組成，每個進水口均為同一側(cè)流向，且3 個部位之間都有30 cm 左右的高度差，魚池的進、出水口均設(shè)置攔網(wǎng)，防止雜質(zhì)進入、水質(zhì)遭到污染以及可以有效的防止幼魚的逃逸。根據(jù)連通器的原理設(shè)計水交換系統(tǒng)。水源為井岡山山泉水，水體清澈無污染，透明度高，符合《中華人民共和國漁業(yè)水質(zhì)標準》（GB 11607—1989）。

1.2 試驗設(shè)計

試驗魚初始體質(zhì)量為（7.89±1.45）g/尾，初始體長為（10.41±0.70）cm/尾。隨機分配到9個相同地點養(yǎng)殖池中。本實驗設(shè)計了3個密度，即低密度（LD）、中密度（MD）、高密度（HD），分別為59 尾/m2、88 尾/m2、118 尾/m2，比例為2∶3∶4。每個密度設(shè)置3 個平行組。按照試驗設(shè)計低密度池每池投放約為1 650 尾，中密度池投放約為2 475 尾，高密度投放約為3 300 尾。養(yǎng)殖試驗時間從2019年7月下旬至10月下旬，試驗持續(xù)90 d。

1.3 飼養(yǎng)管理

試驗期間采用飽食投喂法投喂飼料，每天06:00，12:00，18:00 投喂鱘魚配合飼料。日投喂量為魚體質(zhì)量的1%～3%，根據(jù)魚體質(zhì)量變化，適時調(diào)整飼料的規(guī)格。每天觀察魚類的攝食情況，對投喂量實施調(diào)整。養(yǎng)殖期間，進水口水流流量為100 cm3/s。每天巡塘3 次，觀察魚類狀況，定期檢測水質(zhì)變化，養(yǎng)殖期間水質(zhì)情況見表1。

表1 不同密度組養(yǎng)殖水質(zhì)變化Tab.1 Changes of aquaculture water quality in different density groups

1.4 樣品采集

于養(yǎng)殖實驗的第90天采集鱘魚肌肉樣本。從每個養(yǎng)殖池中隨機撈取10尾魚進行采樣。解剖前，將魚放入100 mg/L 的MS-222 中麻醉，取背部肌肉組織樣品快速放入液氮罐中帶回實驗室保存在-80 ℃冰箱中用于轉(zhuǎn)錄組的測序。

1.5 RNA-Seq測序

1.5.1 總RNA 的提取與質(zhì)量檢測 整個RNA 提取及操作過程使用的槍頭和離心管均經(jīng)過滅菌，實驗器具經(jīng)過0.1%DEPC 水浸泡、高溫高壓滅菌。從凍存管中取出100 mg組織樣品，放入勻漿儀（Servicebio，武漢，中國）研磨至無組織塊，按照Trizol法提取總RNA。

使用Nanodrop 2000（Thermo，USA）檢測RNA 濃度及純度：OD260/280值應(yīng)在1.8～2.2 范圍內(nèi)，OD260/230值范圍為1.5～1.8。RNA 完整性檢測：配制1%的瓊脂糖凝膠，電泳檢測總RNA 的質(zhì)量。將濃度過高的RNA進行適當比例的稀釋。

1.5.2 cDNA 文庫的構(gòu)建和Illumina 測序 使用NEBNext UltraTM RNA Library Prep Ki（tNEB，USA）試劑盒生成測序文庫并對文庫質(zhì)量進行評價，合格后送至北京百邁客生物科技有限公司（BMK，北京，中國）完成轉(zhuǎn)錄組測序工作。。

1.6 測序數(shù)據(jù)處理與分析

1.6.1 測序數(shù)據(jù)的產(chǎn)出和質(zhì)量控制 使用高通量測序平臺對cDNA文庫進行測序，產(chǎn)出大量的高質(zhì)量原始數(shù)據(jù)（raw data），然后通過Perl 程序進行處理獲得清潔數(shù)據(jù)（clean reads），同時計算clean data 的Q20、Q30、GC含量和序列重復(fù)水平。

1.6.2 轉(zhuǎn)錄組數(shù)據(jù)和參考基因組序列比對 原始數(shù)據(jù)經(jīng)過處理后成為clean data，利用Hisat2 軟件將clean data 比對到參考基因組上，然后利用StringTie 把比對上的讀長進行組裝和定量。本研究選用的小體鱘參考基因組（https://www.ncbi.nlm.nih.gov/genome/?term=acipenser+）。

1.6.3 新基因的挖掘 將拼接好的序列與原有的基因組注釋信息比對，探尋原本未被注釋的轉(zhuǎn)錄區(qū)域，發(fā)掘該物種新的轉(zhuǎn)錄本和新基因，達到補充和完善基因組注釋信息的目的。

1.6.4 基因功能的注釋 運用BLAST 軟件將新基因與數(shù)據(jù)庫比對進行注釋，使用以下數(shù)據(jù)庫進行基因功能的注釋：Nr、Nt、Pfam、KOG/COG、Swiss Prot、KO、GO。

1.6.5 基因表達水平的定量 為了讓片段更真實的反映轉(zhuǎn)錄本表達水平，需要將樣品中的比對上基因組的讀長數(shù)目和轉(zhuǎn)錄本長度進行統(tǒng)一化處理。采用FPKM作為衡量轉(zhuǎn)錄本表達水平高低的指標。

1.6.6 差異表達分析 皮爾遜相關(guān)系數(shù)r 作為評價生物學(xué)重復(fù)相關(guān)性的評估指標，DESeq2 軟件分析樣品間的差異表達。參考Lau 等[9]將Fold Change≥1.5 且Pvalue＜0.01 作為篩選標準。通過對差異顯著性P值進行校正得到錯誤發(fā)生率（FDR）。FDR 作為差異表達基因篩選的關(guān)鍵指標。將3 個比較組的差異基因進行維恩圖的繪制。

1.6.7 差異表達基因功主要的GO 富集分析 GO 數(shù)據(jù)庫是一個結(jié)構(gòu)化的標準生物學(xué)注釋系統(tǒng)。GO 注釋系統(tǒng)包括3個主要分支：生物學(xué)過程（BP），分子功能（MF）和細胞組分（CC）。

1.6.8 差異表達基因KEGG 通路富集分析 以KEGG 數(shù)據(jù)庫中通路為單位，與整個基因組相比，差異表達基因顯著富集的通路。選取顯著性q值最小的前20個通路進行分析。

1.7 qRT-PCR驗證

從3 個密度組（高、中、低密度組），每個密度組中取3 個平行樣，共9 個樣本進行qRT-PCR 驗證。選取轉(zhuǎn)錄組文庫篩選獲得的差異倍數(shù)大的基因，包括與生長和應(yīng)激相關(guān)基因水通道蛋白-4（AQP4）、類絲氨酸蛋白酶抑制劑H1（SERPINH1-like）、磷脂酰肌醇3-激酶調(diào)節(jié)亞基（PIK3R1）、丙酮酸激酶（PKM）、轉(zhuǎn)化生長因子β-2（TGFB2）及與代謝相關(guān)的S-腺苷甲硫氨酸合成酶（MAT2a）、核苷二磷酸激酶2（NDKA2）、L-乳酸脫氫酶A（LDHA）、葡萄糖磷酸變位酶1（PGM1）進行實時定量PCR 檢測，根據(jù)轉(zhuǎn)錄組獲得序列使用primei Premier 5.0 軟件設(shè)計引物（表2）并由生工生物工程（上海）股份有限公司合成。熒光定量PCR 反應(yīng)程序為：95 ℃下預(yù)變性10 min，95 ℃下變性15 s，退火溫度60 ℃下反應(yīng)60 s，72 ℃下延伸30 s，共進行40 個循環(huán)數(shù)，熔解的反應(yīng)條件為65～95 ℃，讀板30 s 記錄熒光量。所有的qRT-PCR 重復(fù)3 次。

表2 實時熒光定量PCR引物Tab.2 Primers for real-time quantitative PCR

1.8 數(shù)據(jù)處理

以β-actin為內(nèi)參基因?qū)Φ玫降母鳂颖镜腃t 值做均一化處理，使用2-ΔΔCT方法進行相關(guān)基因定量，采用軟件SPSS 21.0進行顯著性差異分析，Excel 2019進行繪圖。

2 結(jié)果與分析

2.1 測序結(jié)果的產(chǎn)出及質(zhì)控

試驗完成3個密度組（高、中、低密度組），每個密度組取3個平行樣，共9個樣本進行轉(zhuǎn)錄組測序。此次轉(zhuǎn)錄組測序共得到60.70 Gb 的clean data，各樣品的clean data 均達到6.30 Gb，Q20 和Q30 堿基百分比在97.45%和93.42%以上。GC 含量所占比例平均為低密度組為50.47%、中密度組為50.59%、高密度組50.42%。Clean reads 平均為低密度組23.60×106，中密度組21.91×106，高密度組22.20×106。樣品的測序數(shù)據(jù)見統(tǒng)計表3。

表3 測序數(shù)據(jù)統(tǒng)計Tab.3 Statistical table of sequencing data

2.2 轉(zhuǎn)錄組數(shù)據(jù)與參考基因組比對

將本試驗轉(zhuǎn)錄組clean reads 與小體鱘參考基因組進行比對發(fā)現(xiàn)，各樣品的reads 與參考基因組的比對效率在72.64%～74.38%，比對到基因組唯一位置的reads比例在62.58～63.96%，比對到基因組多個位置的reads 在9.43%～10.42%，比對效果大于70%（表4）。結(jié)果表明，可以用小體鱘基因組作為參考基因組進行后續(xù)分析。

表4 測序數(shù)據(jù)與小體鱘參考基因組比對結(jié)果統(tǒng)計Tab.4 Comparison of sequencing data with the reference genome of small-body sturgeon

2.3 新基因的挖掘與分析

基于小體鱘基因組序列作為參考基因組，共發(fā)掘15 549 個新基因。將新基因與數(shù)據(jù)庫比對后進行注釋，最終得到各數(shù)據(jù)庫注釋的新基因數(shù)量統(tǒng)計見表5。結(jié)果發(fā)現(xiàn)，共12 636 個基因得到注釋，2 913 個基因未得到注釋，其中COG 數(shù)據(jù)庫注釋到2 264個新基因，GO 數(shù)據(jù)庫注釋到8 527個新基因，KEGG 數(shù)據(jù)庫注釋到8 086個新基因。

表5 新基因功能注釋結(jié)果統(tǒng)計Tab.5 New gene function annotation result statistics

2.4 差異基因統(tǒng)計分析

對注釋的基因進行差異表達分析發(fā)現(xiàn)，LD vs HD 比較組出現(xiàn)235 個差異基因（DEGs），其中108 個DEGs 基因表達上調(diào)，127 個下調(diào)；LD vs MD 比較組有205 個DEGs，其中111 個基因表達上調(diào)，94 個下調(diào)；MD vs HD比較組有274個DEGs，其中120個基因表達上調(diào)，154個下調(diào)（表6）。

表6 差異表達基因數(shù)目統(tǒng)計表Tab.6 Statistical table of number of differentially expressed genes

將每組差異基因按照要求進行維恩圖繪制發(fā)現(xiàn)，總計611 個DEGs，LD vs HD 和LD vs MD 比較組有33 個共有的DEGs；LD vs HD 和MD vs HD 比較組有40 個共有的DEGs；LD vs MD 和MD vs HD 比較組有31個共有的DEGs；只有1個DEG是3個比較組共同有的（圖1）。

圖1 差異表達基因維恩分析Fig.1 Venn diagram of differentially expressed genes

2.5 差異表達基因GO富集分析

對上述DEGs進行GO分類發(fā)現(xiàn)。LD vs HD比較組，共有148個DEGs獲得了GO注釋信息，其中包括62個上調(diào)和86個下調(diào)的DEGs，歸類到GO三大本體（BP、CC、MF）的59個子類別中。由圖2A可知，BP主要由代謝過程、細胞過程、單有機體過程、生物調(diào)控和對刺激的反應(yīng)組成；CC 主要由細胞、膜、膜部件、細胞組分和細胞器組成；MF 主要由催化活性和結(jié)合組成。在BP 中GTPase 活性的正向調(diào)、DNA 介導(dǎo)的轉(zhuǎn)位和DNA整合富集最顯著；在CC中，肌球蛋白復(fù)合物富集最顯著；在MF中，ATP結(jié)合、連接酶活性、鋅離子結(jié)合和肌動蛋白結(jié)合。

LD vs MD比較組，共有124個DEGs獲得了GO注釋信息，其中包括71個上調(diào)和53個下調(diào)的DEGs，由圖2B可知，代謝過程、細胞過程、單有機體過程和生物調(diào)控占據(jù)著BP大部分；CC主要由細胞、膜、膜部件、細胞組分和細胞器組成；MF主要由催化活性和結(jié)合這兩個子類別組成。在BP中，GTPase活性的正向調(diào)節(jié)、DNA整合、整合素介導(dǎo)的信號通路和甘油三酯生物合成過程富集最顯著；在CC中，肌球蛋白復(fù)合物富集最顯著；在MF中，ATP結(jié)、連接酶活性、鋅離子結(jié)合、肌動蛋白結(jié)合和GTPase激活劑活性富集最顯著。

MD vs HD 比較組，共有173 個DEGs 獲得了GO 注釋信息，其中包括78 個上調(diào)和95 個下調(diào)的DEGs，由圖2C可知，代謝過程、細胞過程、單有機體過程和生物調(diào)控均為BP主要組成；CC主要由細胞、膜、膜部件、細胞組分、細胞器和細胞器部分組成；MF 主要由催化活性和結(jié)合組成。在BP中，整合素介導(dǎo)的信號通路富集最顯著；在CC中，肌球蛋白復(fù)合物和細胞骨架富集最顯著；在MF中，ATP結(jié)合、連接酶活性、鋅離子結(jié)合、和肌動蛋白結(jié)合富集最顯著。

圖2 不同養(yǎng)殖密度組別間差異基因注釋信息Fig.2 Annotation information of differential genes between groups of different breeding densities

2.6 差異表達基因KEGG通路富集分析

對不同組別差別基因進行KEEG 富集分析發(fā)現(xiàn)，LD vs HD 組中，63個DEGs定位到62 個特定的KEGG 代謝途徑，主要顯著富集的通路為糖胺聚糖生物合成-硫酸乙酰肝素/肝素（3 個DEGs）、半胱氨酸與蛋氨酸代謝通路（3 個DEGs）（圖3A）；LD vs MD 組中，66 個DEGs 定位到107 個特定的KEGG 代謝途徑，主要富集的通路為丙氨酸、天冬氨酸和谷氨酸代謝（6 個DEGs）、吞噬體（10 個DEGs）和mTOR 信號通路（7 個DEGs）（圖3B）；MD vs HD 分組中，84 個DEGs 定位到89 個特定代謝途徑，主要富集在糖酵解/糖異生（10 個DEGs）、半乳糖代謝（7個DEGs）、戊糖磷酸途徑（6 個DEGs）、氨基酸的生物合成（8 個DEGs）這4個通路中（圖3C）。

圖3 不同密度組差異表達基因KEGG通路富集散點Fig.3 KEGG pathway rich hubs of differentially expressed genes in different density groups

2.7 差異基因中與生長及應(yīng)激相關(guān)基因的篩選與分析

選擇p 值小且差異倍數(shù)達的差異表達基因，并根據(jù)相關(guān)文獻報道進行篩選。從LD vs HD 分組中篩選獲得生長激素受體（GHR）、轉(zhuǎn)化生長因子β -2（TGB-beta2—like）和WNT1 誘導(dǎo)信號通路蛋白1（WISP1）、類絲氨酸蛋白酶抑制劑H1（SERPINH1-like）以及細胞外酪氨酸蛋白激酶PKDCC（Pkdcc）等與生長相關(guān)的DEGs及與應(yīng)激相關(guān)的熱休克蛋白（HSP70）等；從LD vs MD 組中篩選獲得類胰島素樣生長因子（IGF-1）、肌肉骨骼胚胎核蛋白1（MUSTN1）、整合素α-5（ITGA5）以及熱休克蛋白90（HSP90）；從MD vs HD 分組中篩選獲得α-1，6-甘露糖基轉(zhuǎn)移酶（ALG12）和轉(zhuǎn)化生長因子β-2（TGB-beta2—like）、7-脫氫膽固醇還原酶亞型X2（Dhcr7-X2）、成肌細胞測定蛋白1（MYOD1）、丙酮酸激酶11（PKM）、磷脂酰肌醇3-激酶調(diào)節(jié)亞基（PIK3R1）以及與應(yīng)激相關(guān)的DnaJ 同源亞家族C 第7 成員（DNAJC7）、水通道蛋白4（AQP4）、S-腺苷甲硫氨酸合成酶（MAT2a）、核苷二磷酸脫氫酶2（NDKA2）、L-乳酸脫氫酶A（LDHA）、葡萄糖磷酸變位酶1（PGM1）等（表7）。

表7 差異基因中生長和應(yīng)激相關(guān)基因的篩選Tab.7 Screening of growth and stress-related genes from differential genes

2.8 qRT-PCR驗證

選取9個篩選獲得的差異倍數(shù)大的基因進行qRT-PCR 驗證，結(jié)果發(fā)現(xiàn)，RNA-seq中與生長和應(yīng)激相關(guān)的水通道蛋白-4（AQP4）、類絲氨酸蛋白酶抑制劑H1（SERPINH1-like）、磷脂酰肌醇3-激酶調(diào)節(jié)亞基（PIK3R1）、丙酮酸激酶（PKM）、轉(zhuǎn)化生長因子β-2（TGFB2）及與代謝相關(guān)的S-腺苷甲硫氨酸合成酶（MAT2a）、核苷二磷酸激酶2（NDKA2）、L-乳酸脫氫酶A（LDHA）、葡萄糖磷酸變位酶1（PGM1）在三個密度組間的變化趨勢和qRT-PCR 結(jié)果相一致（圖4），說明轉(zhuǎn)錄組中篩選獲得的差異基因表達結(jié)果可靠，可以用于后續(xù)的基因功能研究。

圖4 qRT-PCR驗證RNA-seq結(jié)果Fig.4 RNA-Seq results were verified by qRT-PCR

3 討論

不同的密度會導(dǎo)致魚類個體生長產(chǎn)生差異，隨著養(yǎng)殖密度的增加，差異越明顯。研究認為高密度環(huán)境下的魚類，受到一系列復(fù)雜因素的作用，最終導(dǎo)致魚體生長速度減緩。本試驗對3個養(yǎng)殖密度，每個養(yǎng)殖密度3 個生物學(xué)重復(fù)，共9 個樣品進行RNA-seq，共得到60.70 Gb clean data，Q20 和Q30 堿基百分比在97.45%和93.42%以上，將轉(zhuǎn)錄數(shù)據(jù)與小體鱘基因組進行比對，比對效率在72.64%～74.38%（核對2.2 轉(zhuǎn)錄組數(shù)據(jù)與參考基因組比對）。根據(jù)轉(zhuǎn)錄組與參考基因組的比對，共發(fā)掘新基因15 549 個，將新基因與數(shù)據(jù)庫進行比對注釋，其中12 636 個基因得到注釋，2 913 個基因未得到注釋。重復(fù)組的設(shè)置是為了盡可能的消除生物學(xué)可變性的影響，對重復(fù)性進行評估發(fā)現(xiàn)，同一密度之間3個重復(fù)組的相關(guān)性均在0.907以上（結(jié)果未顯示），符合要求。

本試驗由3 個比較分組構(gòu)成，分別為：LD vs HD、LD vs MD、MD vs HD，對不同分組的基因表達進行定量，篩選出差異表達基因（DEGs），并進行注釋和功能富集分析。結(jié)果在LD vs HD 分組得到235 個DEGs，LD vs MD分組得到205個DEGs，MD vs HD分組得到274個DEGs。差異基因GO富集分析，可以發(fā)現(xiàn)各比較分組在GO二級分類富集的差異基因個數(shù)和類別均不同。其中，LD vs HD分組的差異基因主要富集在GTPase 活性的正向調(diào)節(jié)、DNA 介導(dǎo)的轉(zhuǎn)位、DNA 整合、肌球蛋白復(fù)合物、ATP 結(jié)合、連接酶活性、鋅離子結(jié)合和肌動蛋白結(jié)合，LD vs MD 分組的差異基因主要富集在GTPase活性的正向調(diào)節(jié)、DNA整合、整合素介導(dǎo)的信號通路和甘油三酯生物合成過程、肌球蛋白復(fù)合物、ATP 結(jié)合、連接酶活性、鋅離子結(jié)合、肌動蛋白結(jié)合和GTPase激活劑活性，MD vs HD 分組的差異基因主要富集在整合素介導(dǎo)的信號通路、肌球蛋白復(fù)合物、細胞骨架、ATP結(jié)合、連接酶活性、鋅離子結(jié)合、肌動蛋白結(jié)合。結(jié)果表明，養(yǎng)殖密度引起的GTPase 活性變化，有助于細胞分化、生長發(fā)育以及免疫應(yīng)答[10-11]，不同放養(yǎng)密度引起的肌肉運動活性變化主要體現(xiàn)在肌球蛋白復(fù)合物、肌動蛋白結(jié)合、ATP結(jié)合以及鋅離子結(jié)合[12]。

通過KEGG 分析發(fā)現(xiàn)，LD vs HD 比較組中，主要顯著富集的通路為糖胺聚糖生物合成-硫酸乙酰肝素/肝素、半胱氨酸與蛋氨酸代謝通路；LD vs MD 比較組中，主要富集的通路為丙氨酸、天冬氨酸和谷氨酸代謝、吞噬體和mTOR 信號通路；MD vs HD 比較組中，主要富集的通路為糖酵解/糖異生、半乳糖代謝、戊糖磷酸途徑、氨基酸的生物合成。在不同的分組中富集的通路不同，可能因為在不同密度下魚類的生長、生理狀況不同，導(dǎo)致影響程度不同。從這些富集通路可以發(fā)現(xiàn)，不同養(yǎng)殖密度顯著影響了機體基礎(chǔ)代謝，這可能是導(dǎo)致魚類個體出現(xiàn)差異的原因。

結(jié)合文獻報道，人工篩選獲得了27條與生長與應(yīng)激相關(guān)基因。其中，轉(zhuǎn)化生長因子β-2（TGFB2）具有調(diào)節(jié)細胞生長和分化的功能[13]；WNT1 誘導(dǎo)信號通路蛋白1（WISP1）是一種富含半胱氨酸的分泌型母細胞蛋白，調(diào)節(jié)細胞反應(yīng)，如細胞生長、分化和生存[14]；胰島素樣生長因子（IGF-1、IGF-2）作為GH/IGFs生長軸的主要成員，調(diào)節(jié)了大部分生長激素的促生長活性，如細胞生長和分化，DNA 和蛋白質(zhì)的合成記憶脂質(zhì)和碳水化合物的代謝等[15]；HSP70、HSP90作為評價應(yīng)激反應(yīng)的生物指標，同時作為分子伴侶能夠維持細胞和蛋白穩(wěn)態(tài)[16]；水通道蛋白-4（AQP4）介導(dǎo)水的跨膜轉(zhuǎn)運，是中樞神經(jīng)系統(tǒng)中含量最高的AQP[17]；丙酮酸激酶（PKM）促進平滑肌細胞表型轉(zhuǎn)化和新生內(nèi)膜增生[18]；S-腺苷甲硫氨酸（MAT2a）廣泛存在于所有組織和體液中，是生物體內(nèi)極其重要的代謝物質(zhì)[19]；NDKA2參與多種細胞活動，如細胞增殖、生長、黏附和分化等[20]；葡萄糖磷酸變位酶1（PGM1）是糖異生過程中的關(guān)鍵酶，參與糖原的分解和合成，在糖代謝過程中起重要的作用[21]；與生長相關(guān)的SERPINH1-like參與膠原蛋白的正確合成與分泌[22]；L-乳酸脫氫酶（LDHA）是一種關(guān)鍵的代謝酶，優(yōu)先催化丙酮酸轉(zhuǎn)化為乳酸，其活性的變化往往與外源環(huán)境的脅迫密切相關(guān)[23]；磷脂酰肌醇3-激酶調(diào)節(jié)亞基（PIK3R1）在細胞的生長和增殖中發(fā)揮重要作用[24]。

本試驗選取9 個篩選獲得的差異倍數(shù)大的基因進行qRT-PCR 驗證，其中與生長和應(yīng)激相關(guān)的水通道蛋白-4（AQP4）、丙酮酸激酶（PKM）、轉(zhuǎn)化生長因子β-2（TGFB2）和與代謝相關(guān)的S-腺苷甲硫氨酸合成酶（MAT2a）、核苷二磷酸激酶2（NDKA2）、葡萄糖磷酸變位酶1（PGM1）在MD 組中表達量最高，這表明MD 組對鱘的生長、應(yīng)激以及代謝更加有利。而與生長相關(guān)的類絲氨酸蛋白酶抑制劑H1（SERPINH1-like）和與代謝相關(guān)的L-乳酸脫氫酶A（LDHA）在LD 組表達量最高，隨著密度的增加其表達量逐漸降低，推測鱘養(yǎng)殖密度會對其生長和代謝過程產(chǎn)生一定程度的影響。在本試驗中PIK3R1在MD組的表達量最低，在HD 中的表達量最高，說明低密度更有利于鱘的正常發(fā)育生長。試驗選取的9 個基因的qRT-PCR 在3 個密度組之間的變化趨勢與RNA-seq 的趨勢一致，說明該轉(zhuǎn)錄組測序數(shù)據(jù)結(jié)果可靠。綜上，不同養(yǎng)殖密度對生長、代謝、應(yīng)激相關(guān)基因的表達有不同程度的影響，且在MD 組養(yǎng)殖對鱘的生長、代謝和應(yīng)激最有利。

研究對不同養(yǎng)殖密度下的西伯利亞雜交鱘魚幼魚進行了轉(zhuǎn)錄組測序，從生物信息學(xué)的角度揭示養(yǎng)殖密度對魚類生長影響的潛在機制，為后續(xù)研究提供了科學(xué)參考。

致謝：江西省現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系專項（JXARS-03）同時對本研究給予了資助，謹致謝意！