基于簡化基因組測序的紅羅非魚低溫體色變異全基因組關(guān)聯(lián)分析

徐鴻飛，陳 詔*，黃彩林，袁宗偉，趙何勇，李 華，楊賓蘭，周大顏，蘇換換，朱華平

(1.廣西壯族自治區(qū)水產(chǎn)引育種中心，廣西 南寧 530031；2.中國水產(chǎn)科學(xué)研究院珠江水產(chǎn)研究所/農(nóng)業(yè)農(nóng)村部熱帶亞熱帶水產(chǎn)資源利用與養(yǎng)殖重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510380)

【研究意義】紅羅非魚(Red tilapia，Oreochromisspp.)隸屬于鱸形目麗魚科，由羅非魚屬間雜交選育獲得，其體色是由色素細(xì)胞分泌的色素及遺傳、環(huán)境、營養(yǎng)和生理等多種因素共同作用的結(jié)果。紅羅非魚肉質(zhì)細(xì)嫩、營養(yǎng)豐富[1-2]，體色艷麗，無肌間刺且腹腔無黑膜，近年來深受消費(fèi)者青睞，在我國和東南亞等國家的養(yǎng)殖面積不斷擴(kuò)大[3]。但在遺傳育種和越冬過程中，紅羅非魚的膚色會發(fā)生變異，進(jìn)而制約了其養(yǎng)殖產(chǎn)業(yè)的良性發(fā)展[4]。簡化基因組測序(RAD-Seq)因具有酶切DNA片段大小一致、文庫構(gòu)建簡單快速、標(biāo)簽密度易于調(diào)節(jié)、成本低廉和準(zhǔn)確性高等優(yōu)點(diǎn)已得到廣泛應(yīng)用[5]；而全基因組關(guān)聯(lián)分析(GWAS)在揭示復(fù)雜性狀的數(shù)量遺傳變異規(guī)律、解析關(guān)鍵基因的遺傳調(diào)控機(jī)制及推動分子輔助育種等方面具有廣闊前景[6-9]。因此，基于RAD-Seq和GWAS技術(shù)分析紅羅非魚低溫體色變異的分子機(jī)制，篩選出與體色變異相關(guān)的分子標(biāo)記，對擴(kuò)大我國紅羅非魚養(yǎng)殖規(guī)模及提升其養(yǎng)殖效益具有重要意義。【前人研究進(jìn)展】傳統(tǒng)的遺傳育種主要根據(jù)不同表型性狀直接進(jìn)行選擇[10]，是一個(gè)耗時(shí)費(fèi)力的繁雜過程。隨著分子生物學(xué)和基因組學(xué)的發(fā)展，魚類遺傳育種研究正經(jīng)歷由傳統(tǒng)選擇育種、雜交育種到基于基因組信息分子育種的轉(zhuǎn)變[11-12]。Vignal等[13]、Liu和Cordes[14]、全迎春等[15]研究認(rèn)為，基于基因組分子育種的核心內(nèi)容是進(jìn)行基因型和表型多態(tài)性研究，而單核苷酸多態(tài)性(SNP)因具有分布廣、可實(shí)現(xiàn)高通量檢測等優(yōu)點(diǎn)已逐漸成為主流分子標(biāo)記。Tsai等[16]基于高密度SNP序列開展了大西洋鮭(Salmosalar)幼魚生長性狀的GWAS和基因組預(yù)測分析。RAD-Seq主要包括復(fù)雜度降低的多態(tài)序列(CRoPS)測序[17]、限制性酶切位點(diǎn)相關(guān)的DNA(RAD)測序[18]和基因分型測序(GBS)[19]，其中RAD分子標(biāo)記測序技術(shù)已在多個(gè)物種中得到廣泛應(yīng)用，通過RAD測序能在大多數(shù)生物中獲得成千上萬的SNP分子標(biāo)記[20-21]。Houston等[22]利用RAD-Seq技術(shù)對傳染性胰臟壞死病病毒(IPN)感染抗性和敏感的大西洋鮭親本及14個(gè)子代進(jìn)行基因分型，并構(gòu)建了遺傳連鎖圖譜。Baird等[23]利用RAD-Seq技術(shù)對3種鏈胞菌屬(Neurosporacrassa)進(jìn)行基因分型，發(fā)現(xiàn)13000多個(gè)SNPs位點(diǎn)。Carmichael等[24]研究顯示，利用RAD-Seq技術(shù)檢測到鮭魚虱(Lepeophtheirussalmonis)不同性別群體中存在Lsa101901標(biāo)記位點(diǎn)，該位點(diǎn)在雌性中為雜合子，在雄性中為純合子。GWAS則在動植物復(fù)雜性狀研究中發(fā)揮著重要作用，可從全基因組范圍內(nèi)篩選出與目標(biāo)性狀相關(guān)聯(lián)的位點(diǎn)[25]，目前已在大西洋鮭(SalmosalarL.)[6]、虹鱒(Oncorhynchusmykiss)[7]、鯉(Cyprinuscarpio)[8]和鯰(Ictalurusspp.)[9]等多種水生動物中得到應(yīng)用。【本研究切入點(diǎn)】目前，從分子水平上揭示紅羅非魚體色變異機(jī)理的研究主要集中于轉(zhuǎn)錄組[3]和紅斑數(shù)量性狀基因座位點(diǎn)繪制[26]，而有關(guān)紅羅非魚在低溫后體色性狀發(fā)生變異的分子標(biāo)記篩選鮮見研究報(bào)道?！緮M解決的關(guān)鍵問題】在RAD-Seq的基礎(chǔ)上，結(jié)合GWAS策略在全基因組范圍內(nèi)尋找與紅羅非魚低溫處理后體色變異性狀相關(guān)的SNP位點(diǎn)，旨在篩選出與紅羅非魚低溫體色變異相關(guān)的遺傳標(biāo)記，為紅羅非魚體色變異的改良與優(yōu)良品種選育提供基礎(chǔ)研究資料。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)所用200尾紅羅非魚由廣西壯族自治區(qū)水產(chǎn)引育種中心桂虹羅非魚良種場提供，平均體重50.11 g/尾，其親本為廣西壯族自治區(qū)水產(chǎn)引育種中心桂虹羅非魚良種場的普通生產(chǎn)群體。試驗(yàn)紅羅非魚在室內(nèi)長方形水泥池(7.6 m×2.3 m×0.8 m)暫養(yǎng)7 d后開始降溫試驗(yàn)，暫養(yǎng)期間水溫保持25 ℃，正常喂食，每2 d換水1次。暫養(yǎng)7結(jié)束后，以每天降低1 ℃的速率進(jìn)行緩慢降溫，降至15 ℃時(shí)維持7 d。7 d后選取體色變化明顯和體色不變的紅羅非魚各30尾，利用100 mg/L MS-222(美國Sigma MO公司)進(jìn)行麻醉處理，刮去鱗片后取其皮膚組織投入液氮中速凍，-80 ℃保存?zhèn)溆谩?/p>

1.2 基因組DNA提取與檢測

根據(jù)傳統(tǒng)的苯酚∶氯仿∶異戊醇法對紅羅非魚皮膚組織樣品進(jìn)行基因組DNA提取，并利用酶標(biāo)儀和Qubit Fluorometer檢測基因組DNA濃度，采用1.0 % 瓊脂糖凝膠電泳檢測基因組DNA的完整性，用微量分光光度計(jì)測定其質(zhì)量，確定基因組DNA純度，并置于-20 ℃冰箱保存?zhèn)溆谩?/p>

1.3 RAD文庫構(gòu)建及測序

RAD文庫的構(gòu)建及測序委托深圳華大基因股份有限公司完成。操作步驟如下：取1 μg基因組DNA，使用EcoR I(NEB)進(jìn)行酶切。對得到的酶切片段進(jìn)行5′端修復(fù)和磷酸化處理，3′端加poly(A)，連接P1 Index接頭，用瓊脂糖凝膠電泳篩選300～500 bp目的DNA，參照QIAquick PCR Purification Kit試劑盒說明進(jìn)行純化回收，對回收產(chǎn)物進(jìn)行末端修復(fù)、末端加poly(A)及P2接頭連接，然后對連接產(chǎn)物進(jìn)行PCR擴(kuò)增。PCR產(chǎn)物經(jīng)瓊脂糖凝膠電泳后，切取目的片段，參照QIAquick Gel Extraction Kit(QIAGEN)試劑盒說明進(jìn)行純化回收，適量Elution Buffer溶解，完成文庫構(gòu)建。構(gòu)建好的RAD文庫使用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System進(jìn)行質(zhì)量和產(chǎn)量檢測，合格文庫在Illumina HiseqTM4000平臺上進(jìn)行測序。

1.4 數(shù)據(jù)過濾與SNP分型

對測序原始數(shù)據(jù)(Raw data)中的低質(zhì)量序列和接頭序列進(jìn)行剪切和過濾得到有效數(shù)據(jù)(Clean data)。以BWA為比對工具，將測序數(shù)據(jù)與參考基因組(https://www.ncbi.nlm.nih.gov/assembly/GCA_001858045.3)進(jìn)行在線比對分析；然后應(yīng)用GATK進(jìn)行SNP位點(diǎn)篩選和基因分型，并采用BGI自主研發(fā)的軟件對其進(jìn)行注釋和統(tǒng)計(jì)分析。對多個(gè)樣品的SNP位點(diǎn)篩選條件為：同一SNP位點(diǎn)上SNP缺失率<50 % ；同一SNP位點(diǎn)上每個(gè)樣品支持該位點(diǎn)的Reads數(shù)≥5；同一SNP位點(diǎn)的堿基類型數(shù)≥2。

1.5 GWAS分析

應(yīng)用一般線性模型(GLM)和混合線性模型(MLM)算法進(jìn)行GWAS分析，并用R語言繪制Q-Q圖(Quantile-Quantile Plot)和曼哈頓圖。其中，曼哈頓圖中的-log(P)閾值選取α=0.05/SNP的個(gè)數(shù)，圖中的虛線為y=-lg(α)；Q-Q圖主要用來估計(jì)數(shù)量性狀觀測值與預(yù)測值間的差異，獲得的數(shù)量性狀數(shù)據(jù)均為正態(tài)分布數(shù)據(jù)。Q-Q圖的X和Y軸分別代表SNP標(biāo)記位點(diǎn)的預(yù)測-lg(P)值和觀察-lg(P)值。

2 結(jié)果與分析

2.1 測序數(shù)據(jù)質(zhì)量及比對結(jié)果

通過Illumina HiSeqTM4000平臺測序得到的Raw data經(jīng)Barcode拆分、過濾獲得Clean data。在所有Reads中質(zhì)量值大于20的堿基占總Reads長度的比例均大于97.44 % ，質(zhì)量值大于30的堿基占總Reads長度的比例均大于94.07 % ，因此可判斷測序數(shù)據(jù)質(zhì)量值高，可用于后續(xù)的信息分析。以尼羅羅非魚(Oreochromisniloticus)的基因組序列作為參考基因組，應(yīng)用BWA軟件將所有測序數(shù)據(jù)與參考基因組進(jìn)行比對，結(jié)果發(fā)現(xiàn)所有樣本的比對率均介于99.34 %～99.45 % ，說明測序數(shù)據(jù)質(zhì)量可靠，可進(jìn)行后續(xù)數(shù)據(jù)分析和試驗(yàn)操作。

2.2 SNP分型結(jié)果

通過BWA與參考基因組進(jìn)行比對分析后，利用GATK軟件進(jìn)行SNP位點(diǎn)篩選和基因分型，結(jié)果共得到1823 228個(gè)SNPs位點(diǎn)。根據(jù)這些SNPs位點(diǎn)的完整度進(jìn)一步篩選，得到46 889個(gè)高質(zhì)量SNPs位點(diǎn)可用于GWAS。

2.3 群體結(jié)構(gòu)分析結(jié)果

通過Admixture軟件分析紅羅非魚的群體結(jié)構(gòu)，根據(jù)交叉驗(yàn)證錯(cuò)誤率的估值確定最優(yōu)分群數(shù)(K)，結(jié)果(圖1)表明，K均在0.6以上，說明所分析的樣品來自同一個(gè)原始祖先，即樣品不存在明顯的分群，適合進(jìn)行后續(xù)的關(guān)聯(lián)分析。從混合線性模型下群體的Q-Q圖(圖2)可看出，實(shí)際觀測值與預(yù)測值在前期基本相符，在后期產(chǎn)生分離，實(shí)際觀測值曲線逐漸翹起，說明所選用的分析模型正確，關(guān)聯(lián)分析結(jié)果可靠。

2.4 GWAS分析結(jié)果

采用混合線性模型對紅羅非魚低溫處理后體色的變異性狀進(jìn)行GWAS分析，并繪制曼哈頓圖。從圖3可看出，GWAS分析共篩選出24個(gè)與紅羅非魚低溫體色變異性狀顯著相關(guān)的SNPs位點(diǎn)(P<5.15E-08)，其中有15個(gè)SNPs位點(diǎn)位于1號染色體上，有9個(gè)SNPs位點(diǎn)位于12號染色體上，且具有成簇分布的特點(diǎn)。取SNPs位點(diǎn)上、下游各1000 bp序列，采用BLAST程序與尼羅羅非魚基因組進(jìn)行在線比對，獲得每個(gè)SNP位點(diǎn)周圍的基因，結(jié)果顯示除5個(gè)未注釋到基因及蛋白功能重復(fù)的基因外，共發(fā)現(xiàn)19個(gè)已知功能的基因，這些基因主要與信號轉(zhuǎn)導(dǎo)、合成和免疫等功能相關(guān)(表1)。

表1 與紅羅非魚低溫體色性狀顯著關(guān)聯(lián)SNPs位點(diǎn)基因的注釋情況Table 1 Gene annotation of SNPs significantly associated with body color trait in red tilapia after low temperature treament

3 討 論

SNP作為一種新興的分子標(biāo)記，在草魚(Ctenopharyngodonidella)[27]、尖吻鱸(Latescalcari)[28]、大口黑鱸(Micropterussalmoides)[29]和牙鲆(Paralichthysolivaceu)[30]等多種魚類中得到廣泛應(yīng)用，是一種有效的分子輔助育種手段，也是優(yōu)良品種選育的基礎(chǔ)。本研究通過對降溫條件下體色變化明顯與體色不變各30尾紅羅非魚進(jìn)行RAD-Seq，共篩選到24個(gè)與紅羅非魚低溫體色變異性狀顯著相關(guān)的SNPs位點(diǎn)，這些SNP位點(diǎn)主要分布在1號染色體和12號染色體上(在其他染色體上均未發(fā)現(xiàn))，且表現(xiàn)出成簇分布的特點(diǎn)，與Li等[31]、Huang等[32]對其他物種基因組的關(guān)聯(lián)分析結(jié)果相似，可能與低溫體色變異相關(guān)的SNP位點(diǎn)主要集中在部分染色體上有關(guān)。

魚類體色主要通過色素細(xì)胞及其所含色素而形成，且受神經(jīng)系統(tǒng)和內(nèi)分泌系統(tǒng)調(diào)節(jié)，能隨著環(huán)境的改變而產(chǎn)生適應(yīng)性變化。Li等[26]研究發(fā)現(xiàn)，馬來西來紅羅非魚遺傳性紅斑形成相關(guān)的微衛(wèi)星主要位于5號染色體和15號染色體上，雖然在這2條染色體上發(fā)現(xiàn)2個(gè)與色素合成相關(guān)的基因(SLC45A2和ASIP)，但這2個(gè)基因未在黑羅非魚和紅羅非魚的皮膚顯著表達(dá)，說明自身色素合成機(jī)制不一定能直接影響紅羅非魚體色變異，而外界環(huán)境刺激也可能對其體色變化產(chǎn)生影響。本研究通過對低溫條件下與體色變異性狀相關(guān)SNP位點(diǎn)附近的基因進(jìn)行GWAS分析，除5個(gè)SNPs位點(diǎn)未注釋到的基因外，共發(fā)現(xiàn)19個(gè)已知功能的基因，這些基因主要與信號轉(zhuǎn)導(dǎo)、生物合成和免疫等功能相關(guān)，與Zhu等[4]在不同體色馬來西亞紅羅非魚皮膚組織中的轉(zhuǎn)錄組測序結(jié)果較一致。此外，所注釋到的基因均與信號轉(zhuǎn)導(dǎo)、物質(zhì)合成及免疫等相關(guān)，究其原因可能是溫度降低對紅羅非魚是一種強(qiáng)烈的外界刺激。紅羅非魚在應(yīng)對外界溫度降低刺激時(shí)引起各種生理反應(yīng)，導(dǎo)致其需調(diào)動大量的基因來應(yīng)答外界刺激，以減少刺激對自身的影響，使與信號轉(zhuǎn)導(dǎo)、合成及免疫應(yīng)答等功能相關(guān)的基因表現(xiàn)較活躍。本研究在所注釋到的基因中未找到直接與紅羅非魚體色變異相關(guān)的基因，可能與用于RAD-Seq的紅羅非魚群體數(shù)量較少有關(guān)，因此，下一步還需在更大的紅羅非魚群體中進(jìn)行篩查。

4 結(jié) 論

給予RAD-Seq和GWAS分析共篩選出24個(gè)與紅羅非魚低溫體色變異性狀顯著相關(guān)的SNPs位點(diǎn)，并注釋到19個(gè)已知的功能基因，且均與信號轉(zhuǎn)導(dǎo)、生物合成及免疫相關(guān)，可作為后續(xù)開展紅羅非魚體色相關(guān)分子輔助育種研究的參考依據(jù)。