白斑狗魚(yú)雌、雄基因組混池重測(cè)序研究

張俊杰，趙瑞陽(yáng)，蔣 麗，胡 瓊，李 菁，唐 露，李勝忠

(1.新疆農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)學(xué)院，烏魯木齊 830052；2.中國(guó)水產(chǎn)科學(xué)研究院，北京 100141)

白斑狗魚(yú)(EsoxLucius)屬鮭形目(Salmoniformes)狗魚(yú)科(Esocidae)狗魚(yú)屬(Esox)，自然分布于北美洲及歐亞大陸北部淡水流域，在我國(guó)僅見(jiàn)于新疆北部的額爾齊斯河流域[1,2]。白斑狗魚(yú)作為自然水體中的主要頂級(jí)捕食者，其生長(zhǎng)迅速，體型較大，具有重要經(jīng)濟(jì)價(jià)值[3]。自人工繁殖取得成功以來(lái)，白斑狗魚(yú)已經(jīng)成為新疆重要的特色養(yǎng)殖魚(yú)類(lèi)，并逐漸向全國(guó)推廣，取得了顯著的經(jīng)濟(jì)效益。

研究表明，白斑狗魚(yú)雌魚(yú)的生長(zhǎng)速度明顯快于雄魚(yú)，而且性成熟時(shí)間更晚，壽命更長(zhǎng)，體型更大[4-7]。Luczynski等[8]通過(guò)性逆轉(zhuǎn)與雌核發(fā)育結(jié)合，確定了白斑狗魚(yú)屬于XX/XY(雌性同型)的染色體決定類(lèi)型，但近年來(lái)對(duì)白斑狗魚(yú)的研究主要集中在生物學(xué)特性、繁殖生物學(xué)、生態(tài)學(xué)等方面，其性別鑒定只能在性成熟后通過(guò)觀察和擠壓生殖孔等方法進(jìn)行[9-13]?；蚪M重測(cè)序技術(shù)能夠在個(gè)體或群體樣本的全基因組水平上篩選出相關(guān)性狀的基因序列差異，已在多種魚(yú)類(lèi)中獲得性別特異分子標(biāo)記[14-17]。本實(shí)驗(yàn)通過(guò)對(duì)雌、雄白斑狗魚(yú)進(jìn)行混池重測(cè)序(BSA-seq)，探究白斑狗魚(yú)雌、雄性別的基因組差異，以期篩選出雌、雄特異分子標(biāo)記，實(shí)現(xiàn)對(duì)白斑狗魚(yú)的性別控制技術(shù)及全雌魚(yú)單性生產(chǎn)。

1 材料與方法

1.1 實(shí)驗(yàn)材料

實(shí)驗(yàn)所用白斑狗魚(yú)購(gòu)自烏魯木齊北園春市場(chǎng)，體重500～800 g，解剖后通過(guò)性腺觀察進(jìn)行生理性別鑒定，然后剪取鰭條，在無(wú)水乙醇中保存，用于后續(xù)基因組DNA的提取。

1.2 實(shí)驗(yàn)試劑

實(shí)驗(yàn)所用試劑主要有：組織裂解液[1 mL Tris-HCl(1 mol/L，pH 8.0)，1 mL EDTA(0.5 mol/L，pH 8.0)，20 mL NaCl(0.5 mol/L)，10 mL SDS(10%)]、乙酸銨(7.5 mol/L)和PCR反應(yīng)預(yù)混液(天根生化科技北京有限公司)。所用引物均由生工生物工程(上海)股份有限公司合成。

1.3 實(shí)驗(yàn)方法

1.3.1 白斑狗魚(yú)DNA制備

選取鑒定性別后的雌、雄白斑狗魚(yú)各20尾，使用醋酸銨法提取基因組DNA，經(jīng)全波長(zhǎng)酶標(biāo)儀(InfiniteM200，瑞士TECAN公司)測(cè)定濃度和A260 nm/A280 nm值后，送至上海凌恩生物科技有限公司進(jìn)行測(cè)序。用同樣方法提取和測(cè)定雌、雄各24尾白斑狗魚(yú)的DNA，并稀釋至50 ng/μL，用于性別特異標(biāo)記的PCR驗(yàn)證。

1.3.2 混池重測(cè)序與數(shù)據(jù)質(zhì)控

將20尾雌性和20尾雄性白斑狗魚(yú)的DNA分別混合，構(gòu)建雌、雄2個(gè)基因組DNA樣品池，將其序列分別隨機(jī)打斷，電泳回收長(zhǎng)度為300～500 bp的DNA片段，兩端加上接頭后進(jìn)行cluster制備，構(gòu)建重測(cè)序文庫(kù)，在Illumina Hiseq 4000平臺(tái)進(jìn)行測(cè)序(上海凌恩生物科技有限公司)。將重測(cè)序獲得的原始測(cè)序數(shù)據(jù)(reads)，通過(guò)Trimmomatic(https://github.com/timflutre/trimmomatic)軟件進(jìn)行過(guò)濾和質(zhì)量檢測(cè)，得到待分析測(cè)序數(shù)據(jù)(clean reads)。

1.3.3 性別標(biāo)記篩選關(guān)聯(lián)分析與差異片段篩選

使用bwa(https://github.com/lh3/bwa)和samtools(https://github.com/samtools/samtools)將重測(cè)序獲得的clean reads與NCBI中的白斑狗魚(yú)參考基因組序列(NCBI登錄號(hào)：GCA_000721915.3)進(jìn)行比對(duì)和排序。將clean reads和參考基因組使用bwa進(jìn)行比對(duì)，再用samtools將比對(duì)后的基因組文件進(jìn)行格式轉(zhuǎn)換和排序，接著使用GATK(https://github.com/broadinstitute/gatk)軟件包對(duì)上一步的結(jié)果文件進(jìn)行SNP和Indel檢測(cè)，通過(guò)MarkDuplicates進(jìn)行去重以及HaplotypeCaller分析，得到SNP和Indel的數(shù)據(jù)集，將比對(duì)質(zhì)量值大于20的SNP位點(diǎn)作為最終SNP位點(diǎn)，并對(duì)得到的每個(gè)SNP和Indel位點(diǎn)進(jìn)行注釋。使用SNP-Index算法進(jìn)行標(biāo)記與性別之間的關(guān)聯(lián)分析，篩選性別特異SNP位點(diǎn)，并對(duì)性別相關(guān)侯選區(qū)域進(jìn)行選擇，在這些候選區(qū)域中進(jìn)行候選基因的注釋?zhuān)S后在STRING(https://string-db.org)網(wǎng)站對(duì)這些基因進(jìn)行功能互作分析。

將比對(duì)到參考基因組上的雌、雄基因組序列進(jìn)行相互比對(duì)，找出差異序列作為性別特異候選序列。同時(shí)將不能比對(duì)到白斑狗魚(yú)參考基因組上的雌、雄基因組序列使用SOAPdenovo(http://soap.genomics.org.cn)進(jìn)行拼接，然后將拼接的雌、雄序列進(jìn)行相互比對(duì)，排除共有序列后，剩余序列也作為該性別特異候選序列，最后逐一將這些特異序列使用BLAST進(jìn)行比對(duì)。

1.3.4 性別特異候選序列的驗(yàn)證

隨機(jī)挑取部分雌、雄候選特異序列，使用primer 5.0設(shè)計(jì)PCR引物，以雌、雄個(gè)體基因組DNA為模板，用PCR擴(kuò)增進(jìn)行兩個(gè)階段的驗(yàn)證。首先用4尾雌魚(yú)和4尾雄魚(yú)的DNA進(jìn)行驗(yàn)證驗(yàn)，再根據(jù)驗(yàn)證結(jié)果，用24尾雌魚(yú)和24尾雄魚(yú)的DNA進(jìn)行驗(yàn)證。PCR反應(yīng)體系如下：DNA模板1 μL，PCR反應(yīng)預(yù)混液6.25 μL，正反引物各0.5 μL(10 μmol/L)，加ddH2O使體積達(dá)到12.5 μL。PCR反應(yīng)程序?yàn)椋?5 ℃預(yù)變性5 min；95 ℃變性30 s，57 ℃退火30 s，72 ℃延伸1 min，35個(gè)循環(huán)；72 ℃再延伸5 min，最后4 ℃保存。用1.5%瓊脂糖凝膠電泳檢測(cè)PCR擴(kuò)增產(chǎn)物。

2 結(jié)果

2.1 混池重測(cè)序數(shù)據(jù)質(zhì)量評(píng)估

雌、雄樣本混池后進(jìn)行重測(cè)序得到的原始測(cè)序讀長(zhǎng)數(shù)目分別為304 730 496條和305 121 064條；堿基數(shù)分別為45 709 574 400 bp和45 768 159 600 bp；堿基測(cè)序錯(cuò)誤率(Q30)分別為90.66%和89.39%；測(cè)序深度50 X。對(duì)reads進(jìn)行過(guò)濾后得到高質(zhì)量的雌、雄混池?cái)?shù)據(jù)，clean reads數(shù)目分別為282 400 688條和277 391 000條；堿基數(shù)分別為41 394 508 708 bp和40 501 950 971 bp；Q30分別提升至94.44%和93.98%，統(tǒng)計(jì)結(jié)果如表1所示。

表1 白斑狗魚(yú)混池重測(cè)序數(shù)據(jù)質(zhì)量統(tǒng)計(jì)Tab.1 Quality statistics of mixed pool resequencing data of E.lucius

2.2 性別特異標(biāo)記篩選

通過(guò)表2可以看到，雌性池SNP位點(diǎn)總計(jì)2 076 647個(gè)，純合SNP位點(diǎn)11 167個(gè)，占0.54%；雌性池Indel位點(diǎn)總計(jì)105 810個(gè)，缺失突變有54 234個(gè)，比插入突變多2 658個(gè)，雌性特異候選序列有524條。雄性池SNP位點(diǎn)總計(jì)2 076 942個(gè)，純合SNP有12 573個(gè)，占0.61%；雄性池Indel位點(diǎn)共計(jì)647 466個(gè)，缺失突變有337 940個(gè)，比插入突變多28 414個(gè)，雄性特異候選序列500條。篩選得到了10個(gè)比對(duì)質(zhì)量值大于20的性別特異SNP位點(diǎn)，其中轉(zhuǎn)換位點(diǎn)7個(gè)，顛換位點(diǎn)3個(gè)。同時(shí)雌性SNP位點(diǎn)的堿基信息和雄性SNP位點(diǎn)堿基信息具有顯著差異，具體情況如表3所示。

表2 白斑狗魚(yú)雌、雄混池SNP和Indel位點(diǎn)以及性別特異候選序列Tab.2 SNP and Indel loci and sex specific candidate sequences in the female and male mixed pool of E.lucius

表3 重測(cè)序獲得的白斑狗魚(yú)性別特異SNP位點(diǎn)Tab.3 Sex specific SNP loci of E.lucius obtained by resequencing

2.3 性別特異候選基因分析

使用SNP-Index算法獲得Δ(SNP-Index)的分布和局部加權(quán)回歸(locally weighted regression，LOESS)擬合曲線(圖1)。經(jīng)bootstrap方法篩選得到26個(gè)差異顯著(P<0.01)的Δ(SNP-Index)區(qū)間作為性別特異候選區(qū)域，SNP數(shù)目介于89～946。26個(gè)候選區(qū)域中存在不同數(shù)目的基因，基因總數(shù)為101，詳細(xì)信息見(jiàn)表4。經(jīng)STRING網(wǎng)站分析發(fā)現(xiàn)有27個(gè)基因存在部分相互作用，這些基因分別在生物體結(jié)構(gòu)發(fā)育、細(xì)胞增殖和能量代謝等方面發(fā)揮作用，其中Pdhx和Prkaa1基因調(diào)控主要進(jìn)行的信號(hào)通路(圖2)。

2.4 性別特異片段的PCR實(shí)驗(yàn)驗(yàn)證

將雌、雄混池?cái)?shù)據(jù)比對(duì)到參考基因組后篩選得到雌性特異序列35條、雄性特異序列14條，未比對(duì)到參考基因組數(shù)據(jù)(unmapping reads)經(jīng)組裝后篩選得到雌性特異序列489條、雄性特異序列486條。在1 024條性別特異片段中隨機(jī)選取了276條可設(shè)計(jì)引物的片段進(jìn)行PCR引物設(shè)計(jì)，其中雌性引物116對(duì)，雄性160對(duì)。通過(guò)對(duì)4尾雌魚(yú)和4尾雄魚(yú)的驗(yàn)證后，篩選出20對(duì)條帶清晰且具有性別偏向性的引物(圖3)。然后又通過(guò)對(duì)24尾雌魚(yú)和24尾雄魚(yú)的驗(yàn)證，最終得到1對(duì)引物(編號(hào)為2043138)(表5)的擴(kuò)增產(chǎn)物出現(xiàn)性別差異，雌性陽(yáng)性比率為4/24，雄性陽(yáng)性比率為16/24，具有顯著的雄性偏向性(圖4)。將白斑狗魚(yú)參考基因組和該序列進(jìn)行BLAST比對(duì)，結(jié)果表明該序列不存在于白斑狗魚(yú)參考基因組，組裝出來(lái)的序列長(zhǎng)度為773 bp(圖5)，其中PCR擴(kuò)增產(chǎn)物長(zhǎng)度為384 bp(圖4)。

圖1 白斑狗魚(yú)染色體SNP標(biāo)記的Δ(SNP-Index)分布和LOESS擬合曲線Fig.1 Δ(SNP-Index)distribution and LOESS fitting curve of SNP loci on the chromosome of E.lucius

表4 白斑狗魚(yú)染色體性別特異候選區(qū)域信息統(tǒng)計(jì)Tab.4 Data statistics of sex special chromosome candidate region of E.lucius

續(xù)表

注：“染色體區(qū)域”中信息依次表示為“染色體編號(hào)：起始?jí)A基數(shù)-結(jié)束堿基數(shù)”

圖2 白斑狗魚(yú)27個(gè)性別特異候選基因的相互作用Fig.2 Interaction of 27 sex special candidate genes of E.lucius藍(lán)綠色線-相互關(guān)系來(lái)源于數(shù)據(jù)庫(kù)；紫色線-相互關(guān)系來(lái)源于實(shí)驗(yàn)測(cè)定；綠色線-預(yù)測(cè)為基因鄰接關(guān)系；紅色線-預(yù)測(cè)為基因融合；藍(lán)色線-預(yù)測(cè)為基因同現(xiàn)；黃色線-文本挖掘；黑色線-共表達(dá)；灰色線-同源蛋白質(zhì)

圖3 3個(gè)性別特異片段在4尾雌性和4尾雄性白斑狗魚(yú)的PCR驗(yàn)證結(jié)果Fig.3 PCR result for 3 sex special fragments in 4 female and 4 male E.lucius性別特異片段編號(hào)從左向右依次為2019708、2022414、2043138

圖4 性別特異片段編號(hào)2043138在24尾雌性和24尾雄性白斑狗魚(yú)的PCR擴(kuò)增結(jié)果Fig.4 PCR result for sex special fragment 2043138 in 24 female and 24 male E.lucius

圖5 白斑狗魚(yú)性別特異片段編號(hào)2043138完整序列Fig.5 Total sequence for sex special fragment 2043138 of E.lucius選中黑色部分為白斑狗魚(yú)性別特異片段編號(hào)2043138的PCR產(chǎn)物區(qū)域。

表5 白斑狗魚(yú)性別特異片段(編號(hào)2043138)引物設(shè)計(jì)Tab.5 Primer for sex special fragment 2043138 of E.lucius

3 討論

魚(yú)類(lèi)性別相關(guān)研究不僅可以加強(qiáng)魚(yú)類(lèi)性別發(fā)育機(jī)制的基礎(chǔ)研究，對(duì)于水產(chǎn)養(yǎng)殖業(yè)來(lái)說(shuō)，進(jìn)行魚(yú)類(lèi)性別控制也具有重要的應(yīng)用價(jià)值[18]。但魚(yú)類(lèi)是低等脊椎動(dòng)物，性染色體分化程度較低，多數(shù)魚(yú)類(lèi)性染色體間差異非常小，通過(guò)核型分析往往很難鑒定魚(yú)類(lèi)性別。由于魚(yú)類(lèi)性別決定方式非常復(fù)雜，而且環(huán)境因素也常作為重要影響因子，和遺傳因子共同作用影響性別的形成過(guò)程，甚至出現(xiàn)生理性別與遺傳性別不一致的情況[19]。通過(guò)篩選魚(yú)類(lèi)基因組中性別特有DNA片段或位點(diǎn)，開(kāi)發(fā)出性別分子標(biāo)記進(jìn)行魚(yú)類(lèi)性別鑒定成為一種簡(jiǎn)便快速的方法。隨著基因組測(cè)序技術(shù)的發(fā)展，開(kāi)發(fā)出混合群體分離分析(又稱(chēng)混池重測(cè)序，Bulked Segregant Analysis，BSA)技術(shù)，為魚(yú)類(lèi)性別研究提供了一種有效方法[20]。林曉煜等[14]結(jié)合混池重測(cè)序和個(gè)體重測(cè)序，篩選出一個(gè)大黃魚(yú)(Larimichthyscrocea)雄性特異SNP位點(diǎn)，Lin等[16]通過(guò)混池重測(cè)序和個(gè)體重測(cè)序，在大黃魚(yú)基因組鑒定出雄性特異的15 bp缺失片段。

本研究利用BSA方法得到雌、雄混池?cái)?shù)據(jù)，篩選得到26個(gè)性別特異候選區(qū)域，包含101個(gè)性別相關(guān)的基因，有27個(gè)基因之間存在不同程度的相互作用。其中Pdhx(pyruvate dehydrogenase complex component X，丙酮酸脫氫酶復(fù)合物X)和Prkaa1(AMP-activated protein kinase catalytic subunit alpha-1，AMP激活蛋白激酶α1基因)與其他基因互作聯(lián)系較多，Pdhx基因[21]催化丙酮酸不可逆脫羧作用生成乙酰輔酶a，同時(shí)形成NADH，而且丙酮酸的還原氧化還與增強(qiáng)細(xì)胞增殖能力有關(guān)。Prkaa1基因[22]參與機(jī)體的能量代謝過(guò)程，在能量平衡的調(diào)節(jié)過(guò)程中起重要作用，且Prkaa1多態(tài)性與肌內(nèi)脂肪含量顯著相關(guān)。而在白斑狗魚(yú)的雌、雄個(gè)體之間存在性別二型性，個(gè)體生長(zhǎng)速率與能量代謝和脂肪含量之間可能存在性別差異，所以推測(cè)Pdhx基因和Prkaa1基因通過(guò)調(diào)控細(xì)胞增殖和能量代謝的過(guò)程與個(gè)體性別發(fā)育存在一定的關(guān)聯(lián)。

經(jīng)過(guò)雌、雄混池序列分析，以及PCR擴(kuò)增的實(shí)驗(yàn)驗(yàn)證，獲得1條源于unmapping reads具有雄性偏向性的序列，全長(zhǎng)773 bp?；蚪M比對(duì)產(chǎn)生unmapping reads的主要原因是由于基因組結(jié)構(gòu)變異、堿基數(shù)目大于組裝軟件所允許的錯(cuò)配數(shù)目和gap數(shù)目、測(cè)序誤差、樣本污染，同時(shí)非模式生物的基因組還可能受到核參考基因組質(zhì)量差、不完全共生體、細(xì)胞器基因組的影響，由此unmapping reads仍可能有新的信息發(fā)現(xiàn)，例如它們可能提供有關(guān)病原體、共生體或參考基因組中缺失的序列/基因的信息[23-25]。將該序列經(jīng)BLAST后未能獲得與其他物種數(shù)據(jù)庫(kù)中相匹配的序列，也未能在Pan等[26]上傳的白斑狗魚(yú)參考基因組數(shù)據(jù)庫(kù)中搜索得到該序列甚至是相似序列，僅和鯽(Carassiusauratus)、恐龍魚(yú)(Erpetoichthyscalabaricus)、金線鲃(Sinocyclocheilusrhinocerous)的uncharacterized protein K02 A2.6-like具有61.11%、63.24%、65.75%的相似度，和南亞野鯪(Labeorohita)的voltage-dependent L-type calcium channel subunit alpha-1D protein具有63.01%的相似度。同時(shí)該序列在雌、雄各15尾白斑狗魚(yú)個(gè)體重測(cè)序數(shù)據(jù)中經(jīng)BLAST比對(duì)后，結(jié)果顯示該序列似乎并不存在于某一特定染色體，而是對(duì)應(yīng)于不同個(gè)體的不同染色體(結(jié)果另文發(fā)表)，其生物學(xué)功能還需要進(jìn)一步研究。

本研究通過(guò)對(duì)白斑狗魚(yú)雌、雄混池基因組的測(cè)序和篩選分析，對(duì)unmapping reads組裝后篩選獲得了1條雄性偏向性序列，可為以后進(jìn)一步對(duì)白斑狗魚(yú)的性別發(fā)育機(jī)制研究提供一定理論基礎(chǔ)。