陳小明 李佳凱 王志勇 蔡明夷 韓 芳 劉賢德

(集美大學(xué)水產(chǎn)學(xué)院, 農(nóng)業(yè)部東海海水健康養(yǎng)殖重點實驗室, 廈門 361021)

基于簡化基因組測序的大黃魚耐高溫性狀全基因組關(guān)聯(lián)分析

陳小明 李佳凱 王志勇 蔡明夷 韓 芳 劉賢德

(集美大學(xué)水產(chǎn)學(xué)院, 農(nóng)業(yè)部東海海水健康養(yǎng)殖重點實驗室, 廈門 361021)

利用Illumina HiSeqTM2500測序平臺, 對通過高溫脅迫實驗篩選得到的20尾耐高溫和20尾不耐高溫的大黃魚(Larimichthys crocea)進行了簡化基因組測序(SLAF-seq), 每個樣本的平均測序深度達到10.26×, 共獲得419211個高質(zhì)量的群體單核苷酸多態(tài)性(SNP)位點 。利用TASSEL軟件的混合線性模型(MLM)進行全基因組關(guān)聯(lián)分析(GWAS), 共篩選到38個與大黃魚耐高溫性狀顯著相關(guān)的SNP位點(P<2.39E–08)。利用BLAST程序定位每個SNP位點在大黃魚基因組中的位置, 并分析其周圍的功能基因。結(jié)果在38個SNPs附近共找到26個已知的功能基因, 這些基因主要與細胞轉(zhuǎn)錄、代謝、免疫等功能相關(guān)。研究結(jié)果可為下一步大黃魚耐高溫分子機制解析及耐高溫品種的選育提供參考。

大黃魚; 高溫脅迫; 簡化基因組測序; 單核苷酸多態(tài)性; 全基因組關(guān)聯(lián)分析

大黃魚(Larimichthys crocea)是我國重要的海洋經(jīng)濟魚類, 在自然海區(qū)分布于30—60 m水深, 適應(yīng)溫度在10—32℃, 最適生長溫度在18—25℃[1,2]。當前, 大黃魚的養(yǎng)殖模式仍以淺海網(wǎng)箱養(yǎng)殖為主,網(wǎng)箱深度為4—6 m。由于水深較淺, 夏季大黃魚處在(或接近)其可耐受高溫的時間較長, 持續(xù)高溫會導(dǎo)致大黃魚生長減緩、抗病力下降, 加上病原生物感染, 常常引發(fā)大黃魚大量發(fā)病死亡。因此, 開展大黃魚耐高溫選育的研究, 對提高養(yǎng)殖大黃魚的度夏成活率具有重要的參考意義。

傳統(tǒng)的育種是依據(jù)性狀的表型進行選擇[3],對于耐高溫品種選育, 常規(guī)的做法是先通過高溫脅迫實驗, 篩選出耐高溫的親魚, 用于進一步選育[4],但這種高溫脅迫實驗, 常常會造成親魚的死亡, 即使沒有死亡, 其生殖細胞也會受到損傷, 對繁殖下一代產(chǎn)生不良影響。隨著分子生物學(xué)和基因組學(xué)的發(fā)展, 魚類的育種研究也正在經(jīng)歷由傳統(tǒng)的選擇育種、雜交育種到基于基因組信息的分子育種的轉(zhuǎn)變[5,6]。基于基因組的分子育種核心內(nèi)容是對基因型和表型多態(tài)性的研究, 其中單核苷酸多態(tài)性(Single Nucleotide Polymorphism, SNP)因其分布廣、可實現(xiàn)高通量檢測等優(yōu)點, 逐漸成為主流分子標記[7—9]。如果能夠通過前期實驗, 篩選出與高溫相關(guān)的SNP標記, 借助分子標記對大黃魚進行耐高溫選育, 則將可以有效提高選擇的準確性, 加速育種進程。

全基因組關(guān)聯(lián)分析(Genome-wide association study, GWAS)旨在從全基因組范圍內(nèi)篩選與性狀關(guān)聯(lián)的SNPs[10]。GWAS分析的成本主要有2個: 第一, 用于分析的樣本數(shù)量；第二, 基因分型的費用。數(shù)量越多, 其分型成本、測定性狀的成本也越高。為減少成本, 只測序極端樣品的GWAS方法被開發(fā)出來了, 如BSR-seq (RNA-seq based BSA)、XP-GWAS (Extreme-phenotype genome-wide association study)等。研究證明XP-GWAS可有效降低基因分型的工作量, 能夠進行低成本高效益的SNP篩選[11]。SNP分型技術(shù)也逐漸由早期階段的中、低通量, 發(fā)展到高通量的基因芯片、重測序技術(shù)[12]。其中SLAF-seq (Specific-locus amplified fragment sequencing)就是一種便宜、高效的SNP發(fā)掘與分型高通量方法[13]。

綜合考慮實驗成本和效率, 本研究采用SLAF-seq技術(shù)結(jié)合XP-GWAS策略在全基因組范圍內(nèi)尋找與大黃魚耐高溫性狀相關(guān)的SNP位點, 旨在尋找新的遺傳標記, 為大黃魚耐高溫性狀的改良提供基礎(chǔ)資料。

1 材料與方法

1.1 實驗材料

本實驗于2014年5—6月在福建省寧德市金鈴水產(chǎn)科技有限公司育苗場進行。所用實驗材料為700尾2月齡大黃魚, 平均體重0.98 g, 平均體長4.87 cm,其親本為寧德三都澳海區(qū)的普通養(yǎng)殖群體(88♀× 50♂)。實驗魚放在室內(nèi)長方形水泥池(2 m×1 m×1 m)中暫養(yǎng)7d后開始升溫實驗, 暫養(yǎng)期間水溫(25.0± 0.3)℃, 正常喂食, 每天換水一次。

1.2 實驗方法

高溫脅迫實驗大黃魚暫養(yǎng)7d后, 參照Diegane等[14]的實驗方法每天升高1℃進行緩慢升溫,當魚開始死亡時, 停止升溫, 維持出現(xiàn)死亡時的溫度1d后, 再按每天升高0.5℃進行升溫。在實驗期間, 連續(xù)充氣, 正常投喂餌料, 每天換水1次, 每2h測量水溫1次, 記錄大黃魚的死亡情況, 包括每條魚的死亡時間和死亡溫度。至大黃魚全部死亡時, 停止加熱升溫。在上述大黃魚幼魚群體熱耐受性實驗中, 采集最先死亡的大黃魚20尾(定義為不耐高溫個體S)與最后死亡的大黃魚20尾(定義為耐高溫個體R), 去除內(nèi)臟取全魚, 固定于95%乙醇中, –20℃保存?zhèn)溆谩?/p>

大黃魚基因組DNA的提取與檢測利用常規(guī)的苯酚鯰氯仿鯰異戊醇法抽提DNA, 1.0%瓊脂糖凝膠電泳檢測基因組DNA的完整性, 紫外分光光度計測定OD260與OD280, 確定DNA的質(zhì)量并將濃度調(diào)至50—100 ng/μL, 確保符合后續(xù)測序分型的要求,–20℃保存?zhèn)溆谩?/p>

測序分型本研究利用大黃魚基因組(Gen-Bank登錄號: GCA_000742935.1)作為參考基因組,采用酶切預(yù)測軟件(北京百邁客生物科技有限公司,北京)對其進行電子酶切預(yù)測。根據(jù)選定的最適酶切方案, 對各樣品基因組DNA分別進行酶切。對得到的酶切片段進行5′端修復(fù)和磷酸化處理, 3′端加A, 連接solexa接頭, 用瓊脂糖凝膠電泳進行片段大小選擇, 通過PCR擴增增大文庫量, 文庫質(zhì)檢合格后使用Illumina HiSeqTM2500測序平臺測序。

統(tǒng)計分析利用接頭序列對測序的原始數(shù)據(jù)進行分類, 得到各個樣品的測序讀段。去除質(zhì)量得分<20的測序讀段, 利用SOAP2[15]軟件將合格的測序讀段比對到大黃魚參考基因組上, 將兩端都比對到參考基因組上的測序讀段用來定義SLAF標簽。選取樣本平均深度在2以上的組來定義SLAF標記。然后根據(jù)SLAF標記, 對40個樣本進行SNP檢測。

使用Plink (v1.07)[16]對得到的SNP進行質(zhì)量控制, 棄去最小檢出率低于80%和最小等位基因頻率低于5%的SNP。使用TASSEL[17]軟件中的混合線性模型進行SNPs和表型性狀的全基因組關(guān)聯(lián)分析,使用模型的公式如下: y=Xα+Qβ+e, 其中, y為表型值向量, X為SNP基因型矩陣, Q為admixture[18]生成的群體結(jié)構(gòu)矩陣(K=1); α為基因型效應(yīng)向量, β為固定效應(yīng)向量, e為隨機誤差向量。采用Bonferroni方法對關(guān)聯(lián)P值進行校正, 得到全基因組顯著性閾值。本研究SNP標記的數(shù)目為419211個, 因此Bonferroni校正的全基因組顯著性閾值為2.39E–08 (0.01/419211)。使用R (https://www.R-project.org)軟件生成關(guān)聯(lián)分析結(jié)果P值的Quantile-Quantile (QQ)圖和曼哈頓圖。

2 結(jié)果

2.1 耐高溫和不耐高溫大黃魚的篩選

通過高溫脅迫實驗觀察到大黃魚開始出現(xiàn)死亡時的溫度為33℃, 對高溫脅迫實驗中大黃魚死亡情況按照發(fā)生死亡時的時間進行歸納整理, 最先死亡的20尾大黃魚(定義為不耐高溫個體S)和最后死亡的20尾大黃魚(定義為耐高溫個體R)的具體死亡時間見表 1。

2.2 測序結(jié)果分析

測序結(jié)果和SLAF標簽統(tǒng)計利用大黃魚基因組為參考基因組進行電子酶切預(yù)測, 最終確定使用RsaⅠ+HaeⅢ酶切組合, 酶切片段長度在264—294 bp的序列定義為SLAF標簽, SLAF標簽在基因組上基本均勻分布, 位于重復(fù)序列區(qū)的SLAF標簽比例為2.52%。

對40個樣品進行測序, 最終共得到98620785條測序讀段, 測序平均Q30為84.71%, 平均GC含量為43.10%。使用SOAP2軟件將各個樣品數(shù)據(jù)比對到大黃魚參考基因組上, 兩端都比對到基因組上為可靠的測序讀段, 以比對到基因組唯一位置的兩端序列來做SLAF標簽開發(fā)。本研究共計開發(fā)146560個SLAF標簽, 每個樣品的平均測序深度為 10.26×。

SNP篩選根據(jù)開發(fā)得到的146560個SLAF標記統(tǒng)計SNP信息, 共得到2423837個群體SNP。根據(jù)MAF>0.05和完整度>0.8進行篩選, 共得到419211個高質(zhì)量群體SNP可用于后續(xù)的GWAS分析。

表 1 20尾耐高溫(R)和20尾不耐高溫(S)大黃魚的死亡時間Tab.1 The information of death time of 20-tails thermal-tolerance (R) and 20-tails thermal-sensitive (S) large yellow croaker

群體結(jié)構(gòu)分析通過admixture軟件分析大黃魚的群體結(jié)構(gòu), 根據(jù)交叉驗證錯誤率的谷值確定最優(yōu)分群數(shù)(K), 如圖 1所示, 當K=1時, 交叉驗證錯誤率最低, 說明本實驗分析的樣品應(yīng)該是來自于同一個原始的祖先, 也就是樣品不存在明顯的分群,適合做后續(xù)的關(guān)聯(lián)分析。繪制了混合線性模型下群體的QQ-plot圖(圖 2)。觀測值與期望值前期基本相符, 在后期翹起, 說明選用的分析模型合適, 關(guān)聯(lián)結(jié)果可靠。

圖 1 每個K值對應(yīng)的交叉驗證誤差Fig.1 The CV value of each K value

圖 2 耐熱性狀的Quantile-Quantile圖Fig.2 QQ-plot for thermal tolerance

圖 3 大黃魚耐熱性狀全基因組關(guān)聯(lián)分析的曼哈頓圖Fig.3 Manhttan plot of GWAS of thermal tolerance in large yellow croaker

全基因組關(guān)聯(lián)分析采用混合線性模型進行全基因組關(guān)聯(lián)分析, 大黃魚耐高溫性狀關(guān)聯(lián)分析的曼哈頓圖如圖 3所示, 根據(jù)關(guān)聯(lián)分析的結(jié)果, 共篩選得到38個與大黃魚耐高溫性狀顯著相關(guān)的SNP位點(P<2.39E–08), 其中有14個SNPs位點位于scaffold scf716上, 有10個SNPs位點位于scaffold scf2上, 有10個SNPs位點位于scaffold scf466上, 有成簇分布的特點。取SNP位點上下游各1000 bp的序列, 采用BLAST程序與大黃魚基因組(Gen-Bank登錄號: GCA_000972845.1)進行比對, E-value的閾值為1e–3, 獲得每個SNP位點周圍的基因,除去3個未注釋到基因以及蛋白功能重復(fù)的, 共發(fā)現(xiàn)26個已知功能的基因, 這些基因主要與細胞轉(zhuǎn)錄、代謝、免疫等功能相關(guān)(表 2)。

表 2 與耐熱性狀顯著關(guān)聯(lián)的SNPs位點Tab.2 SNPs significantly associated with thermal tolerance trait

3 討論

對于GWAS分析, 其成本主要來自于2個方面,樣本數(shù)量和SNP分型方法。樣本數(shù)量可以通過僅分析極端表型個體的策略來減少[11], 而SNP分型方法有很多, 有低通量(CAPS、一代測序)、中等通量(SNaPshot、質(zhì)譜法)以及高通量(基因芯片、重測序)等方法[12]。目前大黃魚尚無商業(yè)用的SNP檢測芯片, 且芯片方法也只能檢測已知的SNP變異。重測序方法費用比較高, 但簡化基因組測序可以有效降低這個費用。目前簡化基因組技術(shù)有多種方法,如RAD、GBS、2b-RAD、dd-RAD、SLAF等[13,19],各種方法均有其優(yōu)勢和劣勢。SLAF-seq技術(shù)可以在減少測序成本的同時, 獲得覆蓋度較高的、能夠在群體間進行分型的標記, 為非模式生物全基因組范圍內(nèi)的SNP標記發(fā)掘提供了高效的方法[13]。目前, 研究人員已利用該方法在玉米、黃瓜上進行全基因組關(guān)聯(lián)分析研究[20,21], 在玉米上發(fā)現(xiàn)1個新的玉米花序分生組織突變體、在黃瓜上找到140個與白粉病抗性相關(guān)的SLAF標簽, 2個熱點區(qū)域?；诔杀竞托实目紤], 本研究采用了前述的研究方案。

通過GWAS分析, 本研究共篩選到38個與大黃魚耐高溫性狀顯著相關(guān)的SNP標記(P<2.39E–08),這些SNP位點主要分布在scf2、scf466和scf716上,表現(xiàn)出成簇分布的特點, 這點在其他物種基因組關(guān)聯(lián)分析中也有類似的結(jié)果[22,23], 可能與高溫相關(guān)的這些SNP位點主要集中在部分染色體上。本研究發(fā)現(xiàn)的與高溫相關(guān)的SNP位點涉及基因的功能主要有細胞轉(zhuǎn)錄、代謝、免疫等, 可能高溫對大黃魚來說是一種刺激, 引起其各種生理反應(yīng), 因而基因轉(zhuǎn)錄、復(fù)制以及免疫應(yīng)答變得非?；钴S。遺憾的是, 在本次篩查出SNP位點周圍并沒有找到直接與高溫相關(guān)的基因, 分析其原因可能與我們所用的簡化基因組測序有關(guān), 導(dǎo)致一些相關(guān)的SNP位點沒有篩查出來, 因而對應(yīng)的基因也沒有發(fā)掘出來。因此,下一步還需要用重測序的方法在大群體中進行進一步的篩查。

[1]Xue B, He Y N, Guo Y M, et al.Research of ecological culture system of Larimichthys crocea [J].Modern Agricultural Sciences and Technology, 2014, (16): 244—249 [薛彬, 何依娜, 郭遠明, 等.大黃魚生態(tài)養(yǎng)殖系統(tǒng)研究.現(xiàn)代農(nóng)業(yè)科技, 2014, (16): 244—249]

[2]Li J K, Wang Z Y, Liu X D, et al.Effects of high temperature on serum biochemical indices of large yellow croaker Larimichthys crocea [J].Marine Science Bulletin, 2015, 34(4): 457—462 [李佳凱, 王志勇, 劉賢德, 等.高溫對大黃魚(Larimichthys crocea)幼魚血清生化指標的影響.海洋通報, 2015, 34(4): 457—462]

[3]Hulata G.Genetic manipulations in aquaculture: a review of stock improvement by classical and modern technologies [J].Genetica, 2001, 111(1—3): 155—173

[4]Ma A J, Huang Z H, Wang X A, et al.The selective breeding of thermal tolerance family and appraisal of performance in turbot Scophthalmus maximus [J].Oceanologia et Limnologia Sinica, 2012, 43(4): 797—804 [馬愛軍, 黃智慧, 王新安, 等.大菱鲆(Scophthalmus maximus)耐高溫品系選育及耐溫性能評估.海洋與湖沼, 2012, 43(4): 797—804]

[5]Xu K, Duan W, Xiao J, et al.Development and application of biological technologies in fish genetic breeding [J].Science China Life Sciences, 2015, 58(2): 187—201

[6]Yue G H.Recent advances of genome mapping and marker-assisted selection in aquaculture [J].Fish and Fisheries, 2014, 15(3): 376—396

[7]Vignal A, Milan D, SanCristobal M, et al.A review on SNP and other types of molecular markers and their use in animal genetics [J].Genetics Selection Evolution, 2002, 34(3): 275—306

[8]Liu Z J, Cordes J F.DNA marker technologies and their applications in aquaculture genetics [J].Aquaculture, 2004, 238(1): 1—37

[9]Quan Y C, Ma D M, Bai J J, et al.SNPs identification in RNA-seq data of largemouth bass (Micropterus salmoides) fed on formulated feed and association analysis with growth trait [J].Acta Hydrobiologica Sinica, 2016, 40(6): 1128—1134 [全迎春, 馬冬梅, 白俊杰, 等.大口黑鱸轉(zhuǎn)錄組SNPs篩選及其與生長的關(guān)聯(lián)分析.水生生物學(xué)報, 2016, 40(6): 1128—1134]

[10]Korte A, Farlow A.The advantages and limitations of trait analysis with GWAS: a review [J].Plant Methods, 2013, 9(1): 29

[11]Yang J, Jiang H, Yeh C T, et al.Extreme-phenotype genome-wide association study (XP-GWAS): a method for identifying trait-associated variants by sequencing pools of individuals selected from a diversity panel [J].The Plant Journal, 2015, 84(3): 587—596

[12]Zhao Q Y, Li X, Zhou D G, et al.SNP genotyping methods for crops in post-genomic era [J].Molecular Plant Breeding, 2010, 8(1): 125—133 [趙瓊一, 李信, 周德貴,等.后基因組時代下作物的SNP分型方法.分子植物育種, 2010, 8(1): 125—133]

[13]Sun X, Liu D, Zhang X, et al.SLAF-seq: An efficient method of large-scale de novo SNP discovery and genotyping using high throughput sequencing [J].PLoS One, 2013, 8(3): e58700, doi: 10.1371/journal.pone.0058700

[14]Diegane N D, Chen Y Y, Lin Y H, et al.The immune response of tilapia Oreochromis mossambicus and its susceptibility to Streptococcus iniae under stress in low andhigh temperatures [J].Fish & Shellfish Immunology, 2007, 22(6): 686—694

[15]Li R, Yu C, Li Y, et al.SOAP2: an improved ultrafast tool for short read alignment [J].Bioinformatics, 2009, 25(15): 1966—1967

[16]Purcell S, Neale B, Todd-Brown K, et al.PLINK: a tool set for whole-genome association and population-based linkage analyses [J].The American Journal of Human Genetics, 2007, 81(3): 559—575

[17]Bradbury P J, Zhang Z, Kroon D E, et al.TASSEL: software for association mapping of complex traits in diverse samples [J].Bioinformatics, 2007, 23(19): 2633—2635

[18]Alexander D H, Novembre J, Lange K.Fast model-based estimation of ancestry in unrelated individuals [J].Genome Research, 2009, 19(9): 1655—1664

[19]Davey J W, Hohenlohe P A, Etter P D, et al.Genomewide genetic marker discovery and genotyping using next-generation sequencing [J].Nature Reviews Genetics, 2011, 12(7): 499—510

[20]Xia C, Chen L, Rong T, et al.Identification of a new maize inflorescence meristem mutant and association analysis using SLAF-seq method [J].Euphytica, 2015, 202(1): 35—44

[21]Zhang P, Zhu Y, Wang L, et al.Mining candidate genes associated with powdery mildew resistance in cucumber via super-BSA by specific length amplified fragment (SLAF) sequencing [J].BMC Genomics, 2015, 16: 1058, doi: 10.1186/s12864-015-2041-z

[22]Li H, Peng Z, Yang X, et al.Genome-wide association study dissects the genetic architecture of oil biosynthesis in maize kernels [J].Nature Genetics, 2013, 45(1): 43—50

[23]Huang X, Wei X, Sang T, et al.Genome-wide association studies of 14 agronomic traits in rice landraces [J].Nature Genetics, 2010, 42(11): 961—967

GENOME-WIDE ASSOCIATION STUDY OF THERMAL TOLERANCE IN LARGE YELLOW CROAKER LARIMICHTHYS CROCEA BASED ON SLAF-SEQ TECHNOLOGY

CHEN Xiao-Ming, LI Jia-Kai, WANG Zhi-Yong, CAI Ming-Yi, HAN Fang and LIU Xian-De
(Key Laboratory of Mariculture for the East China Sea, Ministry of Agriculture; Fisheries College, Jimei University, Xiamen 361021, China)

Twenty thermal-tolerant and twenty thermal-sensitive individuals of Larimichthys crocea were sequenced using specific-locus amplified fragment (SLAF-seq) technology based on Illumina HiSeqTM2500 platform.419211 SNPs were identified with an average read depth of 10.26× for each sample.Thirty-eight SNPs (P<2.39E–08) significantly related with thermal tolerance trait were identified according to association analysis.The SNP locations in large yellow croaker genome were identified using BLAST program, and functional genes around SNP were annotated.Twenty-six genes with known functions were discovered around 38 SNPs, which mainly regulate cell transcription, metabolism and immunity.These results provide basic information to analyze thermal-tolerant molecular mechanism and develop thermal-tolerant lines of Larimichthys crocea in the future.

Larimichthys crocea; High temperature; SLAF-seq; SNP marker; GWAS

Q344+.1

A

1000-3207(2017)04-0735-06

10.7541/2017.91

2016-06-22;

2016-11-05

國家自然科學(xué)基金(31172397和31402339); 福建省高等學(xué)校新世紀優(yōu)秀人才支持計劃(JA14167)資助 [Supported by the National Natural Science Foundation of China (31172397, 31402339); the New Century Excellent Talents of Fujian Province University (JA14167)]

陳小明(1991—), 男, 江西南康人; 碩士研究生; 主要從事水產(chǎn)生物遺傳育種研究。E-mail: x_m_chen@163.com

劉賢德, 教授; 主要從事水產(chǎn)生物遺傳育種研究。E-mail: xdliu@jmu.edu.cn