全基因組關聯(lián)分析揭示白羽肉雞孵化性狀的遺傳基礎

張高猛，丁紀強，劉昱宏，鄭麥青，文 杰，趙桂蘋，李慶賀

(中國農(nóng)業(yè)科學院北京畜牧獸醫(yī)研究所，北京 100193)

雞肉是我國第二大肉類生產(chǎn)和消費品，以孵化性狀為代表的繁殖性能越來越重要[1]，然而長期以來，在白羽肉雞繁殖性能的選育中，育種工作者更加注重母雞的繁殖效率選擇，導致公雞繁殖能力的選擇較少。而種公雞繁殖力的優(yōu)劣直接影響后代的生產(chǎn)性能和肉雞養(yǎng)殖的經(jīng)濟效益。因此，研究與公雞繁殖性狀相關的候選基因?qū)μ岣呷怆u生產(chǎn)的經(jīng)濟效益具有重要意義[2]，家禽育種中，傳統(tǒng)方法難以提高個體繁殖率，而數(shù)量性狀基因座(quantitative trait locus, QTL)潛在候選基因的研究對重要經(jīng)濟性狀影響很大。有研究利用GWAS方法分析了396只京海黃雞核心群母雞11個繁殖性狀相關的SNPs[3]，張濤等[4]利用GWAS方法研究了京海黃雞繁殖性狀相關的分子標記和候選基因，結果發(fā)現(xiàn)有7個SNPs與繁殖性狀相關。由此可見，標記輔助選擇(marker-assisted selection, MAS)可以對繁殖性狀進行改良，從而提高家禽養(yǎng)殖的經(jīng)濟效益。盡管動物的QTL研究長達20多年[5]，但對公雞受精率和受精蛋孵化率的QTL定位研究還比較少。伴隨著基因分型技術的不斷發(fā)展,可以更好的從基因組水平和轉(zhuǎn)錄組水平解析受精率和受精蛋孵化率的遺傳基礎[6]。本研究基于白羽肉雞A、B兩個資源群體，使用本團隊自主研發(fā)的“京芯一號”55K SNP芯片對A群體(556只)和B群體(398只)共954個個體進行基因分型，使用一般線性混合模型(LMM)對受精率和受精蛋孵化率進行全基因組關聯(lián)分析(GWAS)，篩選與孵化性狀相關的SNP位點。

隨著測序技術的不斷發(fā)展，繁殖性狀也被證明與基因表達水平有關[7]。有研究通過深度測序探究了卵巢、垂體和下丘腦miRNA與產(chǎn)蛋性能的關聯(lián)[8]，發(fā)現(xiàn)前列腺素D2和WNT信號轉(zhuǎn)導途徑相關的H-PGDS和WNT2在低精子活力表型公雞睪丸中的表達水平較低[9]。本研究對A群體G9世代8只公雞睪丸組織進行轉(zhuǎn)錄組測序，鑒定孵化性狀的轉(zhuǎn)錄水平差異，進一步解析孵化性狀的遺傳基礎，從而為后期的分子標記輔助選擇奠定基礎。

1 材料與方法

1.1 試驗動物

A、B資源群體采用籠養(yǎng)方式飼養(yǎng)，常規(guī)免疫，各世代營養(yǎng)水平保持不變，為減少其他因素對公雞孵化性狀的影響，每只公雞采精后，對12只母雞進行人工授精，收集10 d種蛋后，經(jīng)甲醛+高錳酸鉀3倍量(即每立方米空間使用甲醛42 mL，高錳酸鉀21 g)熏蒸20 min。第18天照蛋時統(tǒng)計12只母雞的受精率后取均值，受精率為入孵蛋數(shù)中受精蛋所占的比例。21 d出雛時統(tǒng)計12只母雞的孵化率后取均值，孵化率為出雛數(shù)占受精蛋數(shù)的比例。A資源群體5個世代分別有80、76、120、120、160只公雞統(tǒng)計受精率和孵化率，B資源群體3個世代分別有118、120、160只公雞統(tǒng)計受精率和孵化率。對每只公雞翅靜脈采血，用肝素鈉抗凝管保存于-20 ℃，用于基因組DNA的提取。A群體G9世代公雞400日齡時，根據(jù)受精率高低挑選8只公雞，試驗組和對照組各4只，解剖采樣上述睪丸組織，置于液氮凍存，用于RNA的提取。

1.2 遺傳力估計

使用R語言中的shapiro.test函數(shù)對受精率和孵化率表型進行正態(tài)檢驗，符合正態(tài)分布(P>0.05)。分別利用系譜信息和基因組信息構建親緣關系矩陣，使用Asremlv4.1軟件包分別計算受精率的遺傳力，Wald F法估計世代和批次的效應(P<0.01)，因此計算時將世代和批次作為固定效應進行矯正[10]，BLUP模型如下：

y=Xb+Zα+ε

其中,y是表型值向量，X和Z是固定效應和加性遺傳效應的關聯(lián)矩陣，b是固定效應向量，α是隨機加性遺傳效應向量，假設α～N(0，G/Aσ2α)，ε是隨機殘差向量，α～N(0，Iσ2e),G是全基因組標記構建的親緣關系矩陣(G矩陣)，A是基于系譜的血緣關系矩陣(A矩陣)。

1.3 SNP分型質(zhì)控

使用常規(guī)的酚-氯仿法提取血液中的基因組DNA，-80 ℃保存，將質(zhì)量合格的樣品送至北京康普森生物技術有限公司，使用“京芯一號”雞55K SNP芯片進行基因分型，使用plink(V1.9)軟件對芯片基因型數(shù)據(jù)進行質(zhì)控[11],剔除基因型缺失率大于10%的個體、最小等位基因頻率小于5%、樣本檢出率小于90%和分型錯誤的SNP。

1.4 全基因組關聯(lián)分析

采用GEMMA(V0.98.1)軟件(https://github.com/genetics-statistics/GEMMA/releases)中的單性狀混合線性模型對受精率、孵化率性狀進行GWAS分析，該模型包括SNP作為固定因子和個體的親緣關系作為隨機效應[12]，混合線性模型中加入批次作為固定效應，固定效應首先采用R(V3.6.0)軟件中model.matrix()函數(shù)轉(zhuǎn)換成0和1的設計矩陣，然后以類似于協(xié)變量的形式加入GEMMA軟件中。統(tǒng)計模型如下：

y=Wα+xβ+u+;u～MVNn(0,λτ-1K),～MVNn(0,τ-1In)

公式中，y代表表型值向量，W代表協(xié)變量(固定效應)的設計矩陣，包括第一列為1，α代表包括截距的相關系數(shù)向量，x代表標記基因型向量，β代表標記效應，u代表隨機效應向量，代表殘差向量，τ-1代表殘差的方差，λ代表兩個方差組分的比率，K代表由SNPs估計的中心親緣關系矩陣，In代表單位矩陣。MVNn代表n維多元正態(tài)分布。Wald檢驗用于篩選與性狀相關SNP的標準。使用plink(V1.90b)軟件中的參數(shù)-indep-pairwise 25 5 0.2來推斷獨立檢驗的有效SNP數(shù)量，最終推斷出有效獨立檢驗SNP為6 950個。因此全基因組水平顯著閾值和全基因組水平潛在顯著閾值分別為7.19×10-6(0.05/6 950) 和1.44×10-4(1/6 950)。

1.5 RNA提取、質(zhì)控、文庫構建和測序

采用氯仿法提取RNA，使用NanoDrop 2000微量分光光度計檢測RNA濃度，Agilent 2100 Bioanalyzer，Aglient RNA 6000 Nano Kit檢測RNA的樣品濃度和完整性，提取的RNA樣品濃度≥1 000 ng·μL-1，RIN值≥8且28S∶18S≥1.2?？俁NA樣品檢測合格后，根據(jù)不同來源mRNA的特性進行純化。真核生物mRNA 3′末端具有polyA尾結構，選用帶有Oligo(dT)的磁珠進行富集純化；原核生物mRNA不具有該特性，選用試劑盒去除rRNA以獲取更多有效信息。向純化得到的mRNA中加入Fragmentation Buffer使其片斷成為短片段。再以片斷化后的RNA為模板，利用隨機引物進行逆轉(zhuǎn)錄過程，實現(xiàn)cDNA第一鏈合成，并加入2nd Strand Marking Buffer、2nd Strand/End Repair Enzyme Mix合成cDNA第二鏈。后經(jīng)末端修復、加堿基A，加測序接頭，經(jīng)磁珠篩選回收目的片段，并進行PCR擴增，完成整個文庫構建。文庫構建完成后，先使用Qubit3.0進行初步定量，稀釋文庫至1 ng·μL-1，隨后使用Agilent 2100對文庫的insert size進行檢測，insert size符合預期后，使用Bio-RAD CFX96熒光定量PCR儀、Bio-RAD KIT iQ SYBR GRN進行q-PCR，對文庫的有效濃度進行準確定量(文庫有效濃度>10 mol·L-1)，以保證文庫質(zhì)量。質(zhì)量合格的文庫用Illumina平臺進行測序。為了保證后續(xù)分析數(shù)據(jù)的質(zhì)量，對原始序列進行過濾，去除接頭污染的Reads，低質(zhì)量的Reads(Reads中質(zhì)量值Q≤19的堿基占總堿基的50%以上)，去除含N比例大于5%的Reads后，得到高質(zhì)量的Clean Reads，再進行后續(xù)分析。

1.6 差異基因篩選及表達量分析

以A群體低受精率公雞為對照，以Fold change(FC)≥1.5、P-value<0.05且P-adj<0.05作為篩選標準篩選差異基因，兩組公雞睪丸組織中基因表達水平的差異以及差異的統(tǒng)計學顯著性通過火山圖顯示，F(xiàn)PKM[13](Fragments Per Kilobase of transcript per Million fragments mapped)是利用RNA-seq技術定量估計基因表達值的一個非常有效的工具，即每百萬fragments(測序核酸片段，PE中一條fragments對應兩條reads)中來自某一基因每千堿基長度的fragments數(shù)目，其同時考慮了測序深度和基因長度對fragments計數(shù)的影響，是目前最為常用的基因表達水平估算方法。采用FPKM作為衡量基因表達水平的指標，使用DESeq2軟件分析樣品間的差異表達。

1.7 實時熒光定量PCR驗證

cDNA反轉(zhuǎn)錄參照FastKing一步法除基因組cDNA第一鏈合成預混試劑試劑盒(天根，北京)操作說明，反轉(zhuǎn)錄產(chǎn)物于-20 ℃凍存?zhèn)溆?，RT-qPCR以β-actin為內(nèi)參基因，根據(jù)NCBI GenBank中雞THYN1、HMGCLL1、COA6基因和β-actin基因序列設計引物(表1)。擴增體系：2×ChamQ Universal SYBR qPCR Master Mix(Vazyme) 10.0 μL,上、下游引物各0.4 μL(10 μmol·L-1)，cDNA模板2.0 μL，7.2 μL水補足20 μL體系。擴增程序：95 ℃ 30 s；95 ℃ 10 s，60 ℃ 30 s，40個循環(huán)。采用7500型熒光定量PCR儀(ABI，美國)對THYN1、HMGCLL1、COA6基因和β-actin進行分析，每個樣品設置3個技術重復。采用比較Ct值法，即2-ΔΔCt法計算基因相對表達量，使用SAS 9.1.3軟件對每個基因的相對表達量進行雙尾非配對T檢驗，P<0.05表示差異顯著，P<0.01表示差異極顯著。

表1 RT-qPCR所用基因的引物信息

2 結 果

2.1 遺傳力估計

使用ASREML v4.1軟件包構建基于系譜信息的傳統(tǒng)動物模型(A矩陣)，基于基因型信息構建的GBLUP(G矩陣)，分別計算受精率的遺傳力。Wald F統(tǒng)計世代和批次的Wald顯著(P<0.01)?；贏群體5～9世代系譜和556個個體受精率表型值計算得遺傳力為0.21，基于G矩陣計算受精率遺傳力為0.14。

2.2 表型和基因型數(shù)據(jù)質(zhì)控

A群體556個個體受精率表型最大值和最小值分別為0.96和0.28、平均值為0.84、變異系數(shù)為9%；孵化率表型最大值最小值分別為1.00和0.68、平均值為0.88、變異系數(shù)為7.1%；B群體398個個體受精率表型最大值和最小值分別為0.95和0.44、平均值為0.79、變異系數(shù)為9.9%；孵化率最大值和最小值分別為0.98和0.60、平均值為0.85、變異系數(shù)為7.9%。質(zhì)控條件：剔除個體基因型缺失率大于10%的個體、最小等位基因頻率小于5%、樣本檢出率小于90%的SNP。質(zhì)控后，A群體保留32 459個SNPs，平均標記密度為31.7 kb/SNP(表2)；B群體保留38 150個SNPs用于后續(xù)分析，平均標記密度為27.0 kb/SNP(表3)。

表2 質(zhì)控后SNPs標記在A資源群體各染色體上的分布

表3 質(zhì)控后SNPs標記在B資源群體各染色體上的分布

2.3 主成分分析

PCA分析中，主成分1為橫坐標，主成分2為縱坐標，以此做主成分散點圖，由圖1可知，兩資源群公雞各個世代沒有聚攏在一起，說明群體之間出現(xiàn)分層。在以往的全基因組關聯(lián)分析當中，由于祖代或者父母代的遺傳差異導致的群體分層會導致分析結果出現(xiàn)假陽性[14]。利用全基因組標記的遺傳信息對樣本遺傳相關進行估計，比系譜信息更能真實的反映了個體間的相關度，可以通過估計個體間的遺傳相關來矯正群體結構分層對關聯(lián)分析結果的影響，因此在關聯(lián)分析中可以使用主成分1作為協(xié)變量，以校正群體結構對關聯(lián)分析的影響[15]。

a.A群體結構主成分分析圖；b.B群體結構主成分分析圖

2.4 全基因組關聯(lián)分析

受精率和孵化率的全基因組關聯(lián)分析結果(圖2)顯示，有9個位點與受精率顯著相關(P<1.01×10-4)，其中A群體556只個體的受精率篩選到4個顯著位點分別位于12號染色體的COPG1基因，11號染色體的ADAMTS18、CDH8基因，4號染色體的FABP2基因；B群體398只個體受精率篩選到4個顯著位點，分別位于20號染色體的SLCO4A1基因，9號染色體的ATP13A3基因，1號染色體的NELL2、VEZT基因；A、B兩群體基因型數(shù)據(jù)合并后，954只個體受精率篩選到1個顯著位點，位于5號染色體的IGHMBP2基因。B群體398只個體孵化率性狀篩選到3個顯著位點(P<4.57×10-8)，位于18號染色體的TMEM、SCO1、MYH1A基因。顯著SNPs信息見表4。

表4 受精率和孵化率達到5%基因組水平顯著的SNPs位點信息

a.A群體受精率GWAS分析QQ圖與曼哈頓圖；b.B群體受精率GWAS分析QQ圖與曼哈頓圖；c.B群體孵化率GWAS分析QQ圖與曼哈頓圖；d.A、B群體芯片數(shù)據(jù)合并后受精率GWAS分析QQ圖與曼哈頓圖

2.5 轉(zhuǎn)錄組測序樣本表型和測序分析結果統(tǒng)計

受精率高組4個個體受精率表型最大值和最小值分別為0.97和0.93、平均值為0.95；睪丸重最大值為47.34 g、最小值為24.08 g、平均值為36.3 g；受精率低組4個個體受精率表型最大值和最小值分別為0.85和0.79、平均值為0.83；睪丸重最大值為35.9 g、最小值為16.13 g、平均值為26.4 g。轉(zhuǎn)錄組測序分析結果見表5。

表5 8個樣本轉(zhuǎn)錄組測序分析結果統(tǒng)計

2.6 差異表達基因篩選和差異表達分析

以白羽肉雞受精率低組為對照組，以Fold change≥1.5、P-value<0.05且P-adj<0.05作為篩選標準，共篩選出差異表達基因17個，其中上調(diào)基因13個，下調(diào)基因4個。受精率高、低公雞睪丸組織中基因表達水平的差異通過FPKM來定量估計，對受精率高、低組的白羽肉雞差異基因進行表達量分析，發(fā)現(xiàn)低受精率組THYN1、TARBP1、ENSGALG00000036249、ENSGALG00000039805基因表達量較高，高受精率組HMGCLL1、COA6、ENSGALG00000033622基因表達量較高(圖3)。

a.差異表達基因火山圖；b.差異表達基因估計表達值圖。L.受精率低組；H.受精率高組，下同

2.7 高、低受精率公雞睪丸中HMGCLL1、COA6 mRNA的表達

從篩選出的差異表達基因中挑選HMGCLL1、COA6基因進行熒光定量驗證，結果顯示400日齡高受精率公雞睪丸中HMGCLL1 mRNA的表達是低受精率公雞的2.86倍，且差異顯著(P<0.05);并且高受精率公雞睪丸中COA6 mRNA的表達量是低受精率公雞的3倍，且差異顯著(P<0.05)(圖4)。

圖4 HMGCLL1和COA6 mRNA在高、低受精率公雞睪丸中表達

3 討 論

過去十年里，GWAS使用高通量的基因分型技術在全基因組范圍內(nèi)尋找SNPs，探討其與疾病以及復雜性狀的關系[16]。在人類GWAS的研究中，已經(jīng)鑒定出影響糖尿病、肥胖癥、乳腺癌在內(nèi)的幾十種疾病的基因位點[17]。畜禽GWAS的研究中，研究人員也將GWAS應用在雞[18-20]、豬[21-22]、牛[23-24]等重要經(jīng)濟性狀上，與生長性能[25]、免疫[26]、產(chǎn)蛋[27]性能顯著相關的SNPs位點和QTLs區(qū)域相繼被挖掘，而對肉雞孵化性狀的GWAS分析研究較少。

白羽肉雞在我國養(yǎng)禽業(yè)中占比高，其飼料轉(zhuǎn)化率高、蛋白質(zhì)含量高，是一種廣受市場歡迎的禽肉類。生長、繁殖性狀是肉雞養(yǎng)殖業(yè)兩大重要的經(jīng)濟性狀，由于生長性能和繁殖性能選育方向不同，導致肉雞在繁殖性狀中的選育還有很大空間。本研究采用系譜和基因組信息估計了白羽肉雞品系受精率性狀遺傳參數(shù)，受精率遺傳力為0.14～0.21，屬于低遺傳力性狀。本研究的ABLUP和GBLUP估計遺傳力的結果與純系蛋雞報道的結果(受精率遺傳力為0.16～0.33)基本一致。此外，利用全基因組關聯(lián)分析篩選與白羽肉雞受精率、孵化率相關的顯著SNPs位點，利用轉(zhuǎn)錄組測序篩選白羽肉雞受精率高、低組睪丸組織中的差異表達基因，利用生物信息學數(shù)據(jù)庫分析差異表達基因參與的生物學過程，挖掘白羽肉雞孵化性狀相關基因，揭示孵化性狀的遺傳基礎。全基因組關聯(lián)分析結果表明，A群體556只個體的受精率篩選到4個顯著基因，分別是位于12號染色體的COPG1基因、11號染色體的ADAMTS18、CDH8基因、4號染色體的FABP2基因；B群體398只個體受精率篩選到4個顯著基因，分別是位于20號染色體的SLCO4A1基因、9號染色體的ATP13A3基因、1號染色體的NELL2、VEZT基因；B群體398只個體孵化率篩選到3個顯著位點(P<4.57×10-8)，分別位于18號染色體的TMEM、SCO1、MYH1A基因； A、B兩群體芯片數(shù)據(jù)經(jīng)plink軟件合并后，954只個體受精率篩選到1個顯著位點，位于5號染色體的IGHMBP2基因。轉(zhuǎn)錄組測序篩選到HMGCLL1、COA6等差異基因。4號染色體上的FABP2基因與受精率顯著關聯(lián)，有研究發(fā)現(xiàn)血漿中FABP2水平與血清膽固醇、尿酸、肌酐水平成正相關、并與白蛋白、紅細胞和血紅蛋白水平呈負相關，這表明FABP2水平隨著糖尿病、腎病的嚴重程度而增加[28]，該基因可能是影響肉雞公雞孵化性狀的候選基因，同時鑒定出ADAMTS18、NELL2兩個基因與受精率存在顯著關聯(lián)。研究發(fā)現(xiàn)，ADAMTS18參與小鼠生殖道的發(fā)育，2周齡小鼠敲除ADAMTS18后，雄性ADAMTS18/-小鼠的表皮腺整個腺體顯著萎縮，7個月及以上的ADAMTS18/-小鼠的表皮腺體表面出現(xiàn)巨集觀腫脹，最后得出結論ADAMTS18缺乏導致雄性小鼠的預發(fā)性腺發(fā)育不良和纖維化[29]。有研究發(fā)現(xiàn)，睪丸分泌蛋白(NELL2)通過ROS1途徑發(fā)出信號，以調(diào)節(jié)附睪蛋白酶成熟和隨后的附睪蛋白酶(OVCH2)分泌，該蛋白酶處理ADAM3以獲得精子受精能力，最終得出結論，睪丸-附睪精子(NELL2-ROS1-OVCH2-ADAM3)信號傳導是雄性生育力所必需的[30]。很多研究都證明以上基因和雄性生殖有關，但本研究篩選出的NELL2基因與受精率的關系鮮見報道，有待進一步的研究。轉(zhuǎn)錄組測序篩選出的差異基因中，HMGCLL1、COA6調(diào)控雄性受精率的機制尚不清楚，還需進一步研究。THYN1基因?qū)儆诰€粒體黃體酶家族，催化α氧化過程和β氧化過程的脫氫步驟，這些步驟與脂肪酸β氧化有關。其可緩解脂肪沉積和脂質(zhì)代謝紊亂，從而顯著改善家禽的產(chǎn)卵率，證明THYN1與家禽的繁殖性能有關[31]。COL6A3基因會影響蛋白精氨酸甲基轉(zhuǎn)移酶7(Prmt7)在雄性生殖細胞中的表達，從而通過miRNA和目標基因調(diào)節(jié)精子增殖[32]。PARD3B基因與肥胖相關疾病有關[33]，肥胖和糖尿病對雄性的精液品質(zhì)有負面影響，并與睪丸激素水平低有關[34]。精液品質(zhì)會影響到孵化性狀[9]，因此這些基因可能是影響孵化性狀的候選基因。

4 結 論

本研究以白羽肉雞為試驗群體，基于系譜和基因型信息估計的受精率遺傳力分別為0.21和0.14。使用GWAS和RNA-seq技術分析了與白羽肉雞公雞孵化性狀相關的候選基因，GWAS篩選到12個SNPs位點與孵化性狀顯著相關，這些位點位于ADAMTS18、FABP2、NELL2等基因。RNA-seq技術鑒定到HMGCLL1、COA6、THYN1、COL6A3等差異表達基因，這些基因可能是影響白羽肉雞孵化性狀的候選基因。本試驗結果為白羽肉雞后續(xù)基因功能、品種改良等研究奠定了基礎。