不同養殖密度虹鱒幼魚生長和肌肉轉錄組比較分析

朱敏杰,曾 起,程 超,肖 敏,李愛平,簡少卿,趙大顯

( 1.南昌大學 生命科學學院,江西省水產動物資源與利用重點實驗室, 江西 南昌 330031;2.井岡山農業科技園,江西 吉安 343000;3.井岡山市竹佬總冷水魚養殖專業合作社,江西 吉安 343000 )

虹鱒(Oncorhynchusmykiss)屬硬骨魚綱鮭科太平洋鮭屬,是我國重要的冷水性經濟魚類,2020年國內鱒魚產量達3.78×104t[1]。隨著市場需求量增大和養殖產業發展,高密度集約化人工養殖迅速興起,而這種養殖模式也導致養殖過程中疾病的發生,嚴重制約了該產業的健康發展。

在集約化水產養殖中,養殖密度被認為是影響魚類生長性能和養殖產量的重要因素之一[2]。隨著養殖密度增大,魚體生長率、存活率和食物利用率都有所下降[3-5];養殖密度的增大,除降低養殖魚體的生長性能外,還會引起魚體氧化應激反應[4]、皮質醇水平變化[6]以及魚體擁擠脅迫相關基因表達水平發生改變[7]等;而低放養密度由于空間利用效率低下,也可能導致更高的生產成本和利潤率下降。因此,適宜的養殖密度對水產動物健康養殖至關重要。

由于魚類在生態和經濟方面的重要性,利用功能基因組學通過基因調控來研究魚類和非生物環境之間的相互作用關系越來越受到廣泛關注[8-12]。隨著轉錄組測序技術在水產研究領域的廣泛應用,其也被針對性地用于探討魚體受養殖密度誘導的應激反應的調節機制[13-14]。近年來,基于二代高通量測序平臺的RNA-seq技術,可提供大量的序列和豐度信息,逐漸替代了微陣列技術,成為轉錄組分析更好的方法[15]。與傳統的基因芯片技術相比,RNA-seq技術在注釋編碼和非編碼基因上有較大優勢,它不需要知道目標物種的基因信息,就可以對多數生物體的轉錄水平進行分析。母尹楠[16]利用RNA-seq技術研究了感染嗜水氣單胞菌(Aeromonashydrophila)的大黃魚(Larimichthyscrocea)轉錄組,分析得到了1996個差異表達基因(DEGs),包括1133個上調表達基因和863個下調表達基因,其中727個差異表達基因顯著上調,489個差異表達基因顯著下調,這些基因分布在Toll-like通路、JAK-STAT通路、MAPK通路以及參與T細胞受體信號通路中;Schunter等[17]采用RNA-seq技術研究德氏三鰭鳚(Tripterygiondelaisi)交配策略中雄性二態性的分子特征,結果發現,兩種雄性表型間差異表達基因數目比雌雄表型間差異表達基因還多;RNA-seq技術還被用于識別大西洋鮭(Salmosalar)肌肉組織中的miRNA,作為毒理學應激的潛在生物標志物[18]。因此,RNA-seq技術在魚類差異表達基因的篩選和候選基因的發掘上發揮重要功能。

筆者以3種養殖密度條件下的虹鱒幼魚為試驗對象,通過比較分析其生長指標變化,并采用RNA-seq高通量測序技術進行測序和分析,篩選獲得養殖密度和生長應激響應之間的關聯基因信息,為進一步開展虹鱒健康養殖和基因組學研究提供參考數據,也為進一步深入挖掘和開發利用虹鱒功能基因提供科學依據。

1 材料與方法

1.1 試驗設計

養殖試驗在井岡山市竹佬總冷水魚養殖專業合作社進行,所有的虹鱒均由同一批受精卵孵化而來,采用室內流水養殖模式進行試驗。試驗用魚均來自同一個養殖池(記為來源池),經人工捕撈、篩選出規格一致的虹鱒幼魚,放養在規格一致、底面積為1 m2的9個養殖網箱中。試驗水溫17～21 ℃,pH 6.5～7.5,出水口水流速度為1 kg/s,池內水深0.4 m。

試驗用虹鱒幼魚初始平均體質量(4.66±0.26)g,初始平均體長(5.98±0.16)cm。試驗設置3個養殖密度:6.99、4.66、2.33 kg/m3,比例為3∶2∶1。按養殖密度設為低密度組、中密度組、高密度組,每組設置3個平行。低密度組每池放養200尾,中密度組每池放養400尾,高密度組每池放養600尾。試驗魚體健康,規格基本一致。試驗自2019年8月6日開始,每隔8 d從每個網箱中取30尾幼魚進行體長和體質量的測量,養殖32 d后從每個網箱中取3尾虹鱒幼魚肌肉組織混樣作為1個樣品,3個密度組每組3個平行,共9個樣本進行轉錄組測序。

1.2 RNA提取

各組樣品使用TRIzol?試劑(Invitrogen)并按照說明書操作提取總RNA,使用DNase I (TaKaRa)去除基因組DNA,使用Nanodrop ND-2000分光光度計(美國賽默飛)、Aglient 2100 分析儀器對總RNA的純度、濃度和完整性進行檢測。RIN值>7的RNA用于下游試驗。

1.3 mRNA-seq文庫構建和Illumina測序

使用NEBNext UltraTM RNA Library Prep Kit(NEB,美國)試劑盒生成測序文庫。利用帶有Oligo(dT)磁珠純化mRNA,隨機六聚體引物和M-MuLV逆轉錄酶合成第一鏈cDNA,DNA聚合酶I和RNase H合成第二鏈cDNA,AMPure XP系統(美國貝克曼庫爾特有限公司)對文庫片段進行純化;然后使用USER酶(NEB,美國) 與大小合適的接頭連接cDNA,37 ℃孵育15 min,95 ℃孵育5 min后進行PCR擴增并對產物進行純化,利用Agilent Bioanalyzer 2100系統對文庫質量進行評價,合格后使用Illumina-HiSeq 4000平臺進行測序。

1.4 數據處理與分析

1.4.1 生長數據處理與分析

以每隔8 d取樣稱量獲得的體長、體質量作為生長指標,按下式計算生長性能指標。

wWGR=(m2-m1)/m1×100%

RSG=(lnm2-lnm1)/t×100%

CF=100m/L3

式中,wWGR為質量增加率(%),RSG為特定生長率(%/d),CF為肥滿度,m1為初始體質量(g),m2為終末體質量(g),t為養殖時間(d),m為體質量(g),L為體長(cm)。

試驗數據均以平均值±標準誤表示,采用Excel和SPSS 21.0軟件進行單因素方差分析。

1.4.2 測序數據產出和質量控制

基于高通量測序平臺IlluminaHiseq 4000對cDNA文庫進行測序,通過Perl腳本對產出的原始數據進行處理,同時計算處理后的有效數據、Q20、Q30、GC含量和序列重復水平。

1.4.3 轉錄組數據和參考基因組序列比對

利用Hisat2軟件 (http://ccb.jhu.edu/software/hisat2/index.shtml,參數:-p 24-dta-x-1-2-S) 對參考基因組進行比對,再利用StringTie軟件(http://ccb.jhu.edu/software/stringtie/,參數:-p 24-G-o-l)進行組裝和定量。

1.4.4 新基因的挖掘

基于參考基因組序列,將拼接好的序列與原有的基因組注釋信息比對,尋找原來未被注釋的轉錄區,發掘該物種性的轉錄本和新基因,補充和完善基因組注釋信息。使用BLAST軟件 (http://blast.ncbi.nlm.nih.gov/Blast.cgi,參數:-query-db-out-outfmt 6-max-target-seqs 1) 將新基因與數據庫比對進行注釋。

1.4.5 基因功能的注釋

基因功能的注釋基于Nr、Nt、Pfam、KOG/COG、Swiss Prot、KEGG、GO數據庫。

1.4.6 差異表達分析

采用DESeq2軟件(http://www.bioconductor.org/packages/release/bioc/html/DESeq.html)分析樣品組間的差異表達。將差異倍數≥1.5且P<0.05作為篩選標準,對每組差異基因進行維恩圖的繪制。

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

GO數據庫是一個結構化的標準生物學注釋系統。GO注釋系統是一個有向無環圖,包括三個主要分支:生物學過程、分子功能和細胞組分。

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

根據差異表達基因進行KEGG富集分析,并選取q值最小的前20個通路進行分析。

1.4.9 生長與應激相關基因的篩選

根據轉錄組測序結果,結合相關文獻報道,篩選與生長和應激相關的基因。

2 結 果

2.1 養殖密度對虹鱒幼魚生長的影響

養殖32 d后,低密度組終末體質量、質量增加率、特定生長率和肥滿度均大于中密度組和高密度組,但3個組別之間上述4個生長指標數據間差異不顯著(P>0.05)(表1)。

表1 不同養殖密度中虹鱒幼魚生長參數統計

2.2 測序結果的產出及質控

本試驗完成3個密度組(低密度組、中密度組、高密度組),每個密度組3個平行,共9個樣本的轉錄組測序,共得到數據67.41 Gb,各樣品的數據均達到5.99 Gb以上。Q20和Q30堿基百分比分別在98.07%和94.72%以上,各密度組GC堿基占總堿基平均百分比為:低密度組51.07%、中密度組51.33%、高密度組51.08%。有效數據總數(取各密度組3個平行組的有效數據平均值)為低密度組2.694×107、中密度組2.478×107、高密度組2.327×107(表2)。

表2 測序數據統計

2.3 轉錄組數據與參考基因組比對

將本試驗轉錄組獲得的有效數據與虹鱒參考基因組進行比對發現,各樣品的序列與參考基因組的比對效率為91.67%～92.77%,比對效果大于70%。另外,單次匹配的比對率為79.26%～80.06%,多次匹配的比對率為11.82%～13.37%(表3)。

表3 樣品測序數據與虹鱒參考基因組的序列比對結果統計

2.4 新基因的分析

基于虹鱒基因組序列作為參考基因組,使用StringTie軟件對Mapped Reads進行拼接,并與原有的基因組注釋信息進行比較,尋找原來未被注釋的轉錄區,發掘虹鱒的新轉錄本和新基因,從而補充和完善原有的基因組注釋信息。過濾掉編碼的肽鏈過短(少于50個氨基酸殘基)或只包含單個外顯子的序列,共發掘5501個新基因,其中4121個基因得到注釋,1380個基因未得到注釋(表4)。

表4 新基因功能注釋結果統計

2.5 差異表達分析

使用DESeq2軟件對虹鱒養殖高、中、低密度組轉錄組兩兩比較進行差異表達分析(FDR<0.05)發現:低密度組與高密度組相比出現2760個差異表達基因,其中1279個差異表達基因上調,1481個差異表達基因下調;低密度組與中密度組相比有654個差異表達基因,其中307個上調,347個下調;中密度組與高密度組相比有2316個差異表達基因,其中1226個上調,1090個下調(表5)。

表5 差異表達基因數目統計

將每組差異基因繪制維恩圖發現,共4030個非冗余差異表達基因中的38個是3個比較組之間共有的(圖1)。

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

2.6 差異基因GO富集分析

對上述差異表達基因進行GO分類發現,低密度組與高密度組比較共有1943個差異表達基因獲得了GO注釋,包括823個上調和1120個下調的差異表達基因。歸類到GO 3大本體的59個子類別中,其中:生物學過程主要由代謝過程(741個)、細胞過程(950個)、單有機體過程(710個)和生物調控(559個)組成;細胞組分主要由細胞(920個)、膜(530個)、膜部件(463個)、細胞組分(909個)和細胞器(649個)組成;分子功能主要由催化活性(706個)和結合(1067個)組成(圖2a)。

低密度組與中密度組比較共有438個具有顯著性差異表達基因獲得了GO注釋信息,其中包括191個上調和247個下調的差異表達基因。代謝過程(131個)、細胞過程(172個)、單有機體過程(160個)和生物調控(157個)占據著生物學過程涉及的大部分差異表達基因;細胞組分主要由細胞(194個)、膜(151個)、膜部件(137個)、細胞組分(188個)和細胞器(137個)組成,分子功能主要由催化活性(135個)和結合(232個)這兩個子類別組成(圖2b)。

中密度組與高密度組進行比較共有1602個具有顯著性差異表達基因獲得了GO注釋信息,其中包括811個上調和791個下調的差異表達基因。生物學過程中有611個差異表達基因參與代謝過程,755個差異表達基因參與細胞過程,586個差異表達基因參與單有機體過程,426個差異表達基因參與生物調控;細胞組分主要由細胞(680個)、膜(480個)、膜部件(439個)、細胞組分(667個)和細胞器(440個)組成,分子功能主要由催化活性(593個)和結合(844個)組成(圖2c)。

通過進一步分析發現,低密度組與高密度組以及中密度組與高密度組的差異表達基因都主要富集在DNA整合、DNA介導的轉位(生物學過程),褶皺(細胞組分)和細胞骨架的結構組成(分子功能),低密度組與中密度組的差異表達基因主要富集在DNA整合、DNA介導的轉位(生物學過程),褶皺(細胞組分)和DNA結合、氧化還原酶活性作用于成對供體結合或重組(分子功能)過程(圖2)。

2.7 差異基因KEGG富集通路分析

通過KEGG富集通路分析發現,低密度組與高密度組進行比較,920個差異表達基因定位到212個特定的KEGG代謝途徑。低密度組與中密度組相比較,196個差異表達基因定位到112個特定的KEGG代謝途徑;中密度組與高密度組進行對比,743個差異表達基因定位到189個特定代謝途徑。

選取q值最小的20個通路進行分析發現:低密度組與高密度組相比,差異表達基因顯著富集到RNA轉運、剪切體、DNA復制和真核生物中核糖體合成4個通路(圖3a);低密度組與中密度組相比,差異表達基因顯著富集到緊密連接通路(圖3b);中密度組與高密度組相比,差異表達基因顯著富集到剪切體,甘氨酸、絲氨酸和蘇氨酸代謝,氨基酸的生物合成和真核生物中核糖體合成4條信號途徑(圖3c)。

2.8 生長與應激相關基因的篩選

選擇P值小且差異倍數大的差異表達基因,并根據相關文獻報道進行篩選。從低密度組-高密度組中篩選獲得絲氨酸/蘇氨酸蛋白磷酸酶pp1-β催化亞基、胰島素受體2、組蛋白結合蛋白4、無翅型MMTV整合位點家族成員2BA、肌肉骨骼-胚胎核蛋白1b、V型質子ATP酶亞基e1等與生長相關的差異表達基因及顆粒蛋白前體、整合素alpha-V、轉換B細胞復合物亞基、組織蛋白酶、熱休克蛋白等與生長相關的差異表達基因;從低密度組-中密度組分組中篩選獲得肌肉生長抑制素、無翅型MMTV整合位點家族成員2BA等與生長相關的差異表達基因及與應激相關的視黃醛脫氫酶2等差異表達基因;從中密度組-高密度組組中篩選獲得生長相關的差異表達基因,包括骨形態發生蛋白2、骨形態發生蛋白5、轉化生長因子β-2前蛋白、肌肉生長抑制素B、組蛋白結合蛋白4(表6)。

表6 篩選獲得的與生長和應激相關的差異表達基因

3 討 論

3.1 養殖密度對虹鱒幼魚生長的影響

養殖密度作為水產健康養殖的一個重要因子,對魚類生長性能產生重要影響。目前,關于養殖密度影響魚類生長的機制還存在著爭議。普遍認為在高密度環境生活的魚類中,受到多種因素的影響,這些影響因素單獨或者疊加地影響魚類個體的生長,使得養殖密度對生長機制研究變得更為復雜。有些魚類隨著養殖密度的增加,其體質量增長(率)均不斷降低,如地圖魚(Astronotusocellatus)[19]、大菱鲆(Scophthalmusmaximus)[20]、施氏鱘(Acipenserschrenckii)[21];有些魚類高密度養殖條件下生長速度和飼料利用率更好,如大西洋白姑魚(Argyrosomusregius)[22]等;有些魚類養殖密度超過臨界范圍,其生長和存活不受密度的影響,如達氏鰉(Husodauricus)[23]、瓦氏海馬(Hippocampuswhitei)[24]稚魚和塞內加爾鰨(Soleasenegalensis)[25]等。筆者發現,隨著養殖密度的增加,虹鱒幼魚體質量增長率等生長指標均有不同程度的降低,雖然這種影響未呈現顯著性差異(P>0.05),但從轉錄組分析結果可以發現,一些與生長和應激相關基因表達模式在不同養殖組之間發生了顯著性變化,同時經過初步分析還發現,高密度養殖條件下虹鱒幼魚抗氧化酶活性與低、中密度組之間存在顯著差異(數據未列出),因此,綜合分析可以推測,隨著養殖時間的延長,高密度養殖條件下虹鱒幼魚的生長受到負面影響。

3.2 不同養殖密度條件下虹鱒幼魚轉錄組比較分析

為進一步闡明養殖密度對虹鱒幼魚生長性能影響的分子機制,筆者通過RNA-seq技術對3個不同養殖密度虹鱒幼魚肌肉組織進行了高通量轉錄組測序,對數據進行過濾和質控后,通過組裝、拼接和轉錄本功能注釋,得到67.41 Gb 的轉錄組有效數據信息,與其他淡水魚類如江鱈(Lotalota)[26]、黃尾鲴(Xenocyprisdavidi)[27]、達氏鰉[28]、花斑裸鯉(Gymnocypriseckloni)[29]、烏鱧(Channaargus)和斑鱧(C.maculatus)[30]一樣都處于較高水平。這些數據的積累為今后開展虹鱒基因組學、免疫學以及健康養殖提供了參考資料。

從轉錄組質量來看,大部分堿基質量打分能達到或超過Q30。通常Q30值在80%以上則可認為測序質量非常可靠,本試驗9個樣品中Q30值最低為94.72%,最高為95.62%,均大于90%,表明筆者利用RNA-seq技術構建的轉錄組數據庫準確可靠,可用于后續虹鱒基因克隆及功能研究。

通過過濾掉編碼的肽鏈過短(少于50個氨基酸殘基)或只包含單個外顯子的序列,筆者共發掘了5501個新基因,其中4121個基因得到注釋,1380個基因未得到注釋。這些基因的發現進一步補充和完善了虹鱒原有的基因組數據,對這些新基因的深入研究將有助于闡明養殖密度與虹鱒幼魚響應的分子機制。

3.3 生長與應激相關基因在虹鱒幼魚生長中的功能

對轉錄組進行GO和KEGG分析,并結合相關文獻報道,人工篩選獲得了19條與生長和應激相關的基因序列,其中:肌肉生長抑制素和肌肉生長抑制素B在魚類早期發育、肌肉生長和免疫系統中發揮重要作用[31-33];無翅型MMTV整合位點家族成員2BA調控魚類卵巢分化和發育[34];草魚(Ctenopharyngodonidella)進食高淀粉飼料后,絲氨酸/蘇氨酸蛋白磷酸酶pp1-β催化亞基可以使過量的碳水化合物轉換為糖原[35];胰島素受體2廣泛分布于各個組織,在生長發育過程中有重要作用[36];肌肉骨骼-胚胎核蛋白1b是一種小型的核蛋白,參與了肌肉和骨骼的發育和再生[37];骨形態發生蛋白2和骨形態發生蛋白5參與了幾個重要的胚胎過程和成脂前體細胞的誘導分化[38],骨形態發生蛋白5在1齡的團頭魴的肌間骨中廣泛高表達[39];轉化生長因子β-2前蛋白可作為間充質祖細胞和巨噬細胞(如破骨巨噬細胞[40]破骨細胞[41]、軟骨細胞[42])的形態發生誘導因子[43];顆粒蛋白前體具有營養和抗炎活性[44];整合素alpha-V影響著身體和內臟的左右不對稱[45];組織蛋白酶B參與了細菌感染和接種過程中的宿主免疫應答[46];熱休克蛋白60可以在先天免疫系統細胞中引發強效的促炎反應,被認為是細胞受到壓力或受損的危險信號[47];轉換B細胞復合物亞基70作為外周肌動蛋白籠與吞噬體連接支架有著關鍵作用[48]。這些基因在不同養殖密度條件表達模式發現顯著變化,說明它們在虹鱒幼魚應對不同的養殖密度中發揮重要功能。這些基因的發現為進一步深入研究養殖密度對虹鱒生長性能影響的分子機制提供了基礎數據,也為確定生產上適宜的養殖密度提供了參考指標。

4 結 論

本試驗結果顯示,隨著養殖密度的增加,虹鱒幼魚質量增加率等生長指標呈現降低的趨勢。通過RNA-seq技術共獲得了3個養殖密度條件下虹鱒幼魚肌肉組織67.41 Gb轉錄組有效數據,通過序列比對共發掘5501個新基因,其中4121個基因得到注釋,1380個基因未得到注釋。差異表達基因分析發現:低密度組與高密度組以及中密度組與高密度組的差異表達基因都主要富集在DNA整合、DNA介導的轉位(生物學過程)、褶皺(細胞組分)和細胞骨架的結構組成(分子功能),低密度組與高密度組進行比較,920個差異表達基因定位到212個特定的KEGG代謝途徑,中密度組與高密度組進行對比,743個差異表達基因定位到189個特定代謝途徑;低密度組與中密度組的差異表達基因主要富集在DNA整合、DNA介導的轉位(生物學過程),褶皺(細胞組分)和DNA結合、氧化還原酶活性作用于成對供體結合或重組(分子功能)過程,196個差異表達基因定位到112個特定的KEGG代謝途徑。篩選獲得了19條與生長和應激相關的基因序列。此試驗結果可為開展鮭科魚類組學研究提供了基因資源庫,也可為今后虹鱒功能基因的挖掘和健康養殖發提供理論依據。