基于簡化基因組測序的紅羅非魚低溫體色變異全基因組關聯分析

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 基因組DNA提取與檢測

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

1.3 RAD文庫構建及測序

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

1.4 數據過濾與SNP分型

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

1.5 GWAS分析

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

2 結果與分析

2.1 測序數據質量及比對結果

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

2.2 SNP分型結果

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

2.3 群體結構分析結果

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

2.4 GWAS分析結果

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

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

3 討 論

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

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

4 結 論

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