基于簡化基因組測序的大黃魚耐高溫性狀全基因組關聯分析

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

(集美大學水產學院, 農業部東海海水健康養殖重點實驗室, 廈門 361021)

基于簡化基因組測序的大黃魚耐高溫性狀全基因組關聯分析

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

(集美大學水產學院, 農業部東海海水健康養殖重點實驗室, 廈門 361021)

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

大黃魚; 高溫脅迫; 簡化基因組測序; 單核苷酸多態性; 全基因組關聯分析

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

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

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

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

1 材料與方法

1.1 實驗材料

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

1.2 實驗方法

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

大黃魚基因組DNA的提取與檢測利用常規的苯酚鯰氯仿鯰異戊醇法抽提DNA, 1.0%瓊脂糖凝膠電泳檢測基因組DNA的完整性, 紫外分光光度計測定OD260與OD280, 確定DNA的質量并將濃度調至50—100 ng/μL, 確保符合后續測序分型的要求,–20℃保存備用。

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

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

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

2 結果

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

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

2.2 測序結果分析

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

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

SNP篩選根據開發得到的146560個SLAF標記統計SNP信息, 共得到2423837個群體SNP。根據MAF>0.05和完整度>0.8進行篩選, 共得到419211個高質量群體SNP可用于后續的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

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

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

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

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

全基因組關聯分析采用混合線性模型進行全基因組關聯分析, 大黃魚耐高溫性狀關聯分析的曼哈頓圖如圖 3所示, 根據關聯分析的結果, 共篩選得到38個與大黃魚耐高溫性狀顯著相關的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個未注釋到基因以及蛋白功能重復的, 共發現26個已知功能的基因, 這些基因主要與細胞轉錄、代謝、免疫等功能相關(表 2)。

表 2 與耐熱性狀顯著關聯的SNPs位點Tab.2 SNPs significantly associated with thermal tolerance trait

3 討論

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

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

[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 [薛彬, 何依娜, 郭遠明, 等.大黃魚生態養殖系統研究.現代農業科技, 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 [全迎春, 馬冬梅, 白俊杰, 等.大口黑鱸轉錄組SNPs篩選及其與生長的關聯分析.水生生物學報, 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

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

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

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