敖漢細(xì)毛羊羊毛經(jīng)濟(jì)性狀的全基因組關(guān)聯(lián)分析

摘 要： 旨在利用全基因組關(guān)聯(lián)分析探尋敖漢細(xì)毛羊羊毛性狀新的分子標(biāo)記和候選基因。本研究采集1～2周歲的健康敖漢細(xì)毛羊耳組織與羊毛作為試驗(yàn)素材，其中，母羊248只，公羊81只，總計329只。羊毛進(jìn)行性狀測定（包括纖維直徑、自然長度、伸直長度、伸直率），并對表型數(shù)據(jù)進(jìn)行描述性統(tǒng)計和相關(guān)性分析。利用綿羊40K液相SNP芯片對全部個體進(jìn)行基因分型。使用Plink 1.07軟件對芯片數(shù)據(jù)進(jìn)行質(zhì)控，使用GCTA軟件和PopLDdecay軟件對質(zhì)控數(shù)據(jù)進(jìn)行群體結(jié)構(gòu)分析。利用GMEMA混合線性模型對4種羊毛性狀進(jìn)行了全基因組關(guān)聯(lián)分析（genome-wide association study，GWAS），利用在線軟件對候選基因進(jìn)行GO和KEGG富集分析。質(zhì)控后得到329只個體的30 079個SNPs位點(diǎn)用于后續(xù)分析。通過GWAS分析篩選出4個在全基因組上顯著相關(guān)的SNPs位點(diǎn)可能影響羊毛經(jīng)濟(jì)性狀，分別位于1號、6號及8號染色體上。篩選出9個在染色體水平上顯著相關(guān)的SNPs位點(diǎn)可能對羊毛性狀具有潛在意義，分別位于3、5、8、11、18、21、22、25號染色體，尋找到39個可能影響羊毛性狀的候選基因。本研究結(jié)果為后續(xù)探究敖漢細(xì)毛羊羊毛性狀的遺傳機(jī)制及分子育種標(biāo)記開發(fā)提供重要參考。

關(guān)鍵詞： 敖漢細(xì)毛羊；羊毛性狀；全基因組關(guān)聯(lián)分析；候選基因

中圖分類號：S826.2

文獻(xiàn)標(biāo)志碼：A

文章編號：0366-6964（2024）10-4346-14

收稿日期：2024-04-15

基金項目：山東省自然科學(xué)基金面上項目（ZR2020MC164；ZR2022MC094）；山東省研究生教育教學(xué)改革研究資助項目（SDYJSJGC2023061）

作者簡介：林曉坤（1999-），男，山東煙臺人，碩士，主要從事動物遺傳育種與繁殖研究，E-mail：243252756@qq.com

*通信作者：趙金山，主要從事羊遺傳育種與生態(tài)畜牧業(yè)研究，E-mail： 201501005@qau.edu.cn；李和剛，主要從事動物遺傳育種研究，E-mail： 201701018@qau.edu.cn

Genome-Wide Association Study of Wool Economic Traits in Aohan Fine Wool Sheep

LIN" Xiaokun， DU" Mengmeng， ZHOU" Lisheng， HUANG" Zhengang， WANG" Di， ZHOU" Donghui， CAO

Xinxin， HE" Jianning， ZHAO" Jinshan*， LI" Hegang*

（College of Animal Science and Technology， Qingdao Agricultural University， Qingdao 266109， China）

Abstract： The study aimed to identify new molecular markers and candidate genes for wool traits in Aohan fine wool sheep using genome-wide association study. The ear tissue and wool of 1-2 years old healthy Aohan fine wool sheep were collected as experimental materials， 248 ewes and 81 rams， totaling 329. The wool traits were detected （including fiber diameter， natural length， elongation length， and elongation ratio）， and the phenotypic data were subjected to descriptive statistics and correlation analysis. All individuals were genotyped using the 40K liquid phase SNP chip for sheep. The chip data were quality controlled using PLINK 1.07 software， and the quality-controlled data were analyzed for population structure using GCTA software and PopLDdecay software. The GMEMA mixed linear model was used to perform a genome-wide association study （GWAS） on 4 wool traits， and online software was used to conduct GO and KEGG enrichment analyses on candidate genes. After quality control， 329 individuals with 30 079 SNP loci were used for subsequent analysis. Through GWAS analysis， 4 SNP loci were identified as significantly associated with wool economic traits across the genome， located on chromosomes 1， 6， and 8. Additionally， 9 SNP loci were identified as significantly associated with wool traits at the chromosome level， located on chromosomes 3， 5， 8， 11， 18， 21， 22， and 25， and 39 candidate genes were identified as potentially affecting wool traits. The results of this study provide important references for further exploring the genetic mechanisms and molecular breeding markers of wool traits in Baohan fine wool sheep.

Key words： Aohan fine wool sheep; wool traits; genome-wide association study; candidate gene

*Corresponding authors：" ZHAO Jinshan， E-mail： 201501005@qau.edu.cn; LI Hegang， E-mail： 201701018@qau.edu.cn

敖漢細(xì)毛羊是由蒙古羊與高加索細(xì)毛羊、斯達(dá)夫細(xì)毛羊雜交培育而成的毛肉兼用品種，主要分布在內(nèi)蒙古赤峰地區(qū)［1］。敖漢細(xì)毛羊具有適應(yīng)能力強(qiáng)、抗病力強(qiáng)等特點(diǎn)，在羊毛細(xì)度和長度方面表現(xiàn)突出，是一種毛肉兼用細(xì)毛羊品種［2-3］。羊毛性狀決定了羊毛制品的品質(zhì)，對綿羊生產(chǎn)具有重要的經(jīng)濟(jì)意義。對于細(xì)毛羊育種者來說，尋找與羊毛性狀相關(guān)的遺傳標(biāo)記可以提高羊毛品質(zhì)，從而提升產(chǎn)業(yè)價值。

隨著商業(yè)化SNP芯片和基因分型技術(shù)的發(fā)展，全基因組關(guān)聯(lián)分析（genome-wide association study，GWAS）已經(jīng)成為檢測畜禽復(fù)雜數(shù)量性狀候選基因的關(guān)鍵方法［4］。例如，有研究利用OvineHD BeadChip在美利奴羊和美利奴雜交綿羊中進(jìn)行了單性狀GWAS、多性狀GWAS和基因組預(yù)測的方法得到了FOXI3［5］、BMPER［6］等基因并發(fā)現(xiàn)它們與毛囊發(fā)育和羊毛發(fā)生有關(guān)［7］。隨著GWAS分析在細(xì)毛羊中的不斷應(yīng)用，也得到了許多與羊毛性狀相關(guān)的基因，例如PTPN3基因與羊毛產(chǎn)量有關(guān)［8］，SLIT3和ZNF280B基因與羊毛纖維直徑有關(guān)［9］，KAP6-1［10］、KAP22-1［11］、FST［12］、IGF-1［13］等基因均對羊毛不同性狀有影響。

雖然細(xì)毛羊羊毛性狀的相關(guān)基因研究較多，但羊毛性狀作為受多基因調(diào)控的復(fù)雜數(shù)量性狀，目前缺乏主效基因和SNP位點(diǎn)，給分子標(biāo)記輔助選擇帶來較大困難。并且毛囊生長發(fā)育受到多種基因和多個信號通路的參與，共同調(diào)控毛囊的周期性發(fā)育過程，研究毛囊生長發(fā)育的規(guī)律和結(jié)構(gòu)特征對羊毛生長發(fā)育至關(guān)重要。因此，本研究通過對敖漢細(xì)毛羊羊毛經(jīng)濟(jì)性狀的免疫相關(guān)候選基因關(guān)聯(lián)分析和全基因組關(guān)聯(lián)分析，探究羊毛生長和毛囊發(fā)育的分子機(jī)制，為細(xì)毛羊分子育種提供重要參考。

1 材料與方法

1.1 試驗(yàn)動物

本研究試驗(yàn)動物由內(nèi)蒙古自治區(qū)赤峰市敖漢旗細(xì)毛羊保種場提供，采集1～2周歲的敖漢細(xì)毛羊耳組織與羊毛樣本，母羊248只，公羊81只，總計329只敖漢細(xì)毛羊樣本。

1.2 耳組織采集與基因組DNA提取

通過耳缺鉗在外耳廓邊緣采集敖漢細(xì)毛羊耳組織，組織大小約為1 cm（長1 cm，厚度0.5 cm），將組織保存在無水乙醇中，利用苯酚-氯仿法提取基因組DNA［14］。提取后，通過微量分光光度計來評估DNA樣品的質(zhì)量和濃度，并進(jìn)行凝膠電泳檢測。將全部DNA樣本-20℃保存用于后續(xù)綿羊40K液相芯片測序。

1.3 羊毛收集與性狀測定

在對應(yīng)羊只的肩胛后緣一掌、體側(cè)靠左中線稍上處采集毛樣（羊毛采集標(biāo)準(zhǔn)，ZB B45005-88），用剪刀緊貼皮膚剪下，帶回實(shí)驗(yàn)室進(jìn)行測量。

測定方法：1）利用顯微鏡測定纖維直徑。隨機(jī)選取羊毛纖維放在蓋玻片上，測量8根毛樣數(shù)據(jù)并記錄，并計算平均值用于后續(xù)全基因組關(guān)聯(lián)分析。2）測量羊毛的自然長度。將毛樣放在黑絨板上，沿毛叢平行方向用小鋼尺精確測取長度，精確到0.1 cm。測量8根毛樣數(shù)據(jù)并記錄，計算平均值用于后續(xù)全基因組關(guān)聯(lián)分析。3）測量羊毛的伸直長度。將毛樣放在黑絨板上，使用鑷子由毛叢根部抽出纖維，直至彎曲消失為止，并準(zhǔn)確測量長度至0.1 cm。測量8根毛樣數(shù)據(jù)并記錄，計算平均值用于后續(xù)全基因組關(guān)聯(lián)分析。4）計算伸直率的公式為：（伸直長度-自然長度）/自然長度［15］。

對得到的表型數(shù)據(jù)進(jìn)行描述性統(tǒng)計分析，剔除異常數(shù)據(jù)，用R語言性狀相關(guān)性分析判斷表型數(shù)據(jù)是否符合正態(tài)分布，對符合正態(tài)分布的表型數(shù)據(jù)用于后續(xù)的GWAS分析。

1.4 SNP分型和質(zhì)量控制

綿羊40K液相芯片利用來源于全世界共587只綿羊及其野生近緣種的重測序數(shù)據(jù)，其中中國羊434只，有更多中國地方綿羊品種特有的中高頻SNP位點(diǎn)，測序深度高，分型準(zhǔn)確，位點(diǎn)在全基因組上均勻分布。本研究采用的敖漢細(xì)毛羊是特有品種，使用該芯片性價比高，并能找到更多SNPs位點(diǎn)。用Plink 1.07［16］軟件對芯片數(shù)據(jù)進(jìn)行質(zhì)控，質(zhì)控標(biāo)準(zhǔn)為：剔除檢出率小于95%的樣本，并過濾質(zhì)量不合格的SNPs：包括檢出率低于90%，最小等位基因頻率低于5%，以及Hardy-Weinberg平衡檢驗(yàn)Plt;1×10-6［17］。

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

通過基于個體基因組SNPs差異程度進(jìn)行成分分析（principal component analysis，PCA），將個體按主成分聚類成不同亞群。本試驗(yàn)使用GCTA軟件對經(jīng)過濾后的SNPs數(shù)據(jù)進(jìn)行PCA分析，生成.eigenval和.eigenvec文件，并在R中繪制PCA圖，以主成分1和主成分2作為x軸和y軸；使用GCTA［18］軟件構(gòu)建基因組G矩陣，可以更真實(shí)地反映個體之間的親緣關(guān)系［19］；使用PopLDdecay（https：∥github.com/BGI-shenzhen/PopLDdecay）軟件對質(zhì)控后的SNPs數(shù)據(jù)進(jìn)行LD衰減分析［20］。比較LD衰減（到0.1）距離和標(biāo)記間的平均距離，以判斷標(biāo)記是否對全基因組具有足夠的覆蓋度。

1.6 全基因組關(guān)聯(lián)分析

全基因組關(guān)聯(lián)分析（GWAS）以博瑞迪綿羊40K液相芯片測序數(shù)據(jù)和329只敖漢細(xì)毛羊羊毛表型數(shù)據(jù)為研究對象。芯片測序數(shù)據(jù)經(jīng)質(zhì)控后得到30 079個位點(diǎn)，使用GEMMA（v0.93）的混合線性模型［21-22］對4種羊毛性狀進(jìn)行GWAS分析，該模型考慮了種群分層和樣本結(jié)構(gòu)［23］，

y=Wα+xβ+u+ε，

u～MVNn（0，λτ-1K），ε～MVNn（0，τ-1In）

其中，y是n個個體性狀（纖維直徑、自然長度、伸直長度、伸直率）的n向量；W=（w1， w2，…，wc）是協(xié)變量（固定效應(yīng)）的n×c矩陣，其中一列為1；α是包含截距的相應(yīng)系數(shù)的c向量；x是標(biāo)記基因型的n載體；β為標(biāo)記物的效應(yīng)大??；u是隨機(jī)效應(yīng)的n向量；ε是誤差的n向量；τ-1為殘差方差；λ為兩個方差分量之比；K是已知的n×n關(guān)聯(lián)矩陣，In是n×n單位矩陣。MVNn表示n維多元正態(tài)分布。

關(guān)聯(lián)矩陣K的計算公式如下：

K=1P∑Pt=1xi-1nx-ixi-1nx-iT

其中，xi為第i列，表示第i個SNP的基因型;x-i為樣本均值;1n為1的n×1向量。GEMMA對每個SNP檢驗(yàn)備選假設(shè)H1，β≠0和零假設(shè)H0，β=0為了糾正多重假設(shè)檢驗(yàn)。用R語言［24］繪制曼哈頓圖和QQ圖將結(jié)果可視化，利用GEC［25］軟件確定閾值線。

1.7 GO和KEGG分析

利用在線軟件Metascape（https：∥metascape.org/gp/index.html）和DAVID線上網(wǎng)站（https：∥david.ncifcrf.gov/home.jsp）對候選基因進(jìn)行GO和KEGG富集分析，尋找候選基因共同參與的信號通路，并結(jié)合候選基因國內(nèi)外相關(guān)研究綜合分析候選基因的生物學(xué)功能。

2 結(jié) 果

2.1 質(zhì)量控制和性狀數(shù)據(jù)分析

綿羊40K液相芯片測序得到的40 156個SNPs位點(diǎn)經(jīng)過數(shù)據(jù)質(zhì)量控制后，得到了329只個體的30 079個SNPs位點(diǎn)。描述性統(tǒng)計分析是對一組數(shù)據(jù)進(jìn)行特征性分析，旨在詳細(xì)描述樣本的測量特征以及所代表總體的特征。這為構(gòu)建統(tǒng)計模型奠定了基礎(chǔ)，并通過觀察表型數(shù)據(jù)的值來初步分析羊毛表現(xiàn)形式信息，從而描繪群體整體特點(diǎn)。本研究敖漢細(xì)毛羊的纖維直徑、自然長度、伸直長度和伸直率4個羊毛性狀的描述性統(tǒng)計結(jié)果如表1所示。此外，通過對羊毛纖維直徑和其他羊毛性狀進(jìn)行相關(guān)性分析發(fā)現(xiàn)，羊毛纖維直徑與羊毛自然長度、羊毛伸直長度、伸直率存在較弱的相關(guān)性（圖1）。

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

2.2.1 主成分分析

將敖漢細(xì)毛羊群體進(jìn)行主成分分析發(fā)現(xiàn)樣本群體比較分散，根據(jù)第一主成分和第二主成分可看出采集的敖漢細(xì)毛羊出現(xiàn)了群體分層現(xiàn)象，如圖2所示。這提示在下一步的全基因組關(guān)聯(lián)研究中需要考慮該結(jié)果。因此，將主成分前10個結(jié)果作為協(xié)變量，用于后續(xù)全基因組關(guān)聯(lián)分析，來校正群體結(jié)構(gòu)。

2.2.2 親緣關(guān)系矩陣分析

通過對329只敖漢細(xì)毛羊的親緣關(guān)系分析得到了親緣矩陣的熱圖，如圖3所示，樣本為同一品種且親緣關(guān)系較近。

2.2.3 LD衰減分析

敖漢細(xì)毛羊群體LD衰減分析結(jié)果如圖4所示，通過LD衰減速度判斷群體多樣性的差異，一般野生群體的LD衰減速度快于馴化群體，LD衰減距離（0.1）和標(biāo)記間平均距離的比較判斷GWAS關(guān)聯(lián)分析中標(biāo)記是否足夠，從圖4中LD衰減距離（0.1）和標(biāo)記間平均距離可以判斷用于GWAS分析的標(biāo)記足夠。

2.3 全基因組關(guān)聯(lián)分析

2.3.1 羊毛纖維直徑全基因組關(guān)聯(lián)分析

羊毛纖維直徑是影響羊毛產(chǎn)品關(guān)鍵因素之一，也是評估細(xì)毛羊經(jīng)濟(jì)價值的重要指標(biāo)。在羊毛纖維直徑的全基因組關(guān)聯(lián)分析結(jié)果中（圖5a），發(fā)現(xiàn)有1個SNPs效應(yīng)位點(diǎn)達(dá)到了全基因組顯著水平，位于1號染色體（rs204762729，P=2.71×10-5），根據(jù)綿羊基因組（GCF_000298735.2_Oar_v4.0）上的基因注釋得到了TTC14和PEX5L兩個基因。此外，Sahana等［26］定義在染色體范圍顯著性水平的SNPs位點(diǎn)作為染色體范圍上顯著性相關(guān)的標(biāo)準(zhǔn)，并且這種顯著性相關(guān)的SNPs位點(diǎn)對應(yīng)的P值應(yīng)小于10-4。因此，得到了3個SNPs位點(diǎn)在染色體水平顯著相關(guān)，分別位于21（rs6974831，P=6.25×10-5）、22（rs25513088，P=8.52×10-5）、25（rs39039643，P=7.88×10-5）號染色體上（表2）。

2.3.2 自然長度性狀全基因組關(guān)聯(lián)分析

羊毛產(chǎn)品中自然長度的重要性僅次于纖維直徑，在相同品質(zhì)下，較長的自然長度意味著更好的羊毛品質(zhì)。在全基因組關(guān)聯(lián)分析中未發(fā)現(xiàn)與自然長度顯著相關(guān)的SNPs位點(diǎn)（圖5b）。但發(fā)現(xiàn)了2個在染色體水平上達(dá)到顯著相關(guān)的SNPs位點(diǎn)，分別位于3（rs55043579，P=5.38×10-5）和11（rs51282011，P=5.29×10-5）號染色體（表2）。

2.3.3 伸直長度性狀全基因組關(guān)聯(lián)分析

伸直長度是評估羊毛品質(zhì)的另一個重要標(biāo)準(zhǔn)，依據(jù)全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)了2個與全基因組顯著相關(guān)的SNPs位點(diǎn)（圖5c），均位于8號染色體（rs30842589，P=7.05×10-7；rs30837398，P=1.40×10-5），如表2所示，根據(jù)綿羊基因組（GCF_000298735.2_Oar_v4.0）上的注釋基因信息，未得到注釋基因。

2.3.4 伸直率性狀全基因組關(guān)聯(lián)分析

在伸直率的全基因組關(guān)聯(lián)分析中，顯著相關(guān)的SNPs位點(diǎn)" 有3個（圖5d），分別在5（rs102267048，P=7.98×10-5）、6（rs103863192，P=3.68×10-5）、18號（rs41285606，P=6.97×10-5）染色體上，如表2所示，且在顯著位點(diǎn)未找到注釋基因。

2.3.5 GO分析和KEGG分析結(jié)果

為進(jìn)一步探究羊毛經(jīng)濟(jì)性狀的相關(guān)分子標(biāo)記以及功能基因，結(jié)合綿羊基因組（GCF_000298735.2_Oar_v4.0）和NCBI上的公共信息，在顯著位點(diǎn)前后500 kb范圍附近找到了TTC14、PEX5L、SORCS3、RAB38、LAMP3、MCF2L2、CCSER2、MCCC1等39個候選基因（表3），并進(jìn)行了功能富集來驗(yàn)證這些候選基因的功能，富集結(jié)果顯示共有21個GO項達(dá)到了具有統(tǒng)計學(xué)意義的水平（Plt;0.05），如圖6所示。這些GO條目主要涉及外肽酶活性、抗原產(chǎn)生與運(yùn)輸?shù)冗^程，其中外肽酶活性的GO條目（GO：0008238）是富集結(jié)果中最顯著的。此外，使用DAVID線上網(wǎng)站（https：∥david.ncifcrf.gov/home.jsp）對候選基因的信號通路查找結(jié)果顯示，MSX1基因參與毛囊發(fā)育相關(guān)通路；ATG5、RPTOR共同參與自噬相關(guān)通路，如圖7所示。

3 討 論

羊毛作為細(xì)毛羊的一個重要性狀，在以往的研究中通過高密度基因分型和GWAS對其他品種的細(xì)毛羊檢測了羊毛經(jīng)濟(jì)性狀的遺傳標(biāo)記和候選基因。在目前的研究結(jié)果中，基因USP13、NLGN1、MTR、DKK1、KIF16B、FGF5等都與羊毛性狀有關(guān)［27-31］。本研究中，綿羊40K液相芯片測序位點(diǎn)也包括對USP13、NLGN1兩個基因的注釋，但在研究結(jié)果中并未關(guān)聯(lián)到，我們推測由于羊毛性狀是受遺傳和環(huán)境因素控制的復(fù)雜數(shù)量性狀，不同的細(xì)毛羊品種關(guān)聯(lián)結(jié)果會有差異。

關(guān)于羊毛性狀分析結(jié)果，先前研究報道的其他細(xì)毛羊品種羊毛性狀描述性統(tǒng)計結(jié)果［32］與本研究中敖漢細(xì)毛羊羊毛性狀的描述性統(tǒng)計結(jié)果基本一致，這就說明本研究所測得羊毛性狀的觀測值基本符合敖漢細(xì)毛羊的羊毛性狀標(biāo)準(zhǔn)（國家標(biāo)準(zhǔn)，標(biāo)準(zhǔn)號：GB/T 25242-2010），因此，本研究測定的表型數(shù)據(jù)是比較可靠的。此外，本研究還對4個羊毛性狀表型進(jìn)行了相關(guān)性分析，性狀相關(guān)性結(jié)果顯示，自然長度與伸直長度、伸直率有強(qiáng)相關(guān)性，這就提示了我們在尋找與羊毛自然長度相關(guān)的候選基因時可以在伸直長度和伸直率兩個性狀的關(guān)聯(lián)結(jié)果中篩選。

在GWAS分析中，模型的選擇對于分析結(jié)果至關(guān)重要，關(guān)系到結(jié)果的準(zhǔn)確性。我們認(rèn)為群體分層和個體親屬關(guān)系是影響假陽性關(guān)聯(lián)的主要因素。統(tǒng)計學(xué)家和遺傳學(xué)家已經(jīng)開發(fā)了幾種方法來克服這兩個影響因素［33-34］。因此，本研究使用GEMMA軟件，采用了全基因組高效混合模型關(guān)聯(lián)（GEMMA）分析，該模型可將群體分層和個體親緣關(guān)系兩因素包含在內(nèi)進(jìn)行GWAS分析。本研究的QQ圖表明，大多數(shù)位點(diǎn)的實(shí)際值和理論值是一致的，這反過來表明方法在分層群體控制方面是有效的。羊毛性狀值通常受性別影響。因此，我們還將性別作為固定因素，以減少性別對特定表型的影響。

根據(jù)以往研究表明，由于Bonferroni方法確定閾值有較強(qiáng)的限制性。可能會導(dǎo)致一些有用的位點(diǎn)被篩除，出現(xiàn)假陰性的結(jié)果［35］。根據(jù)我們的研究結(jié)果，按照Bonferroni方法確定的閾值線為-lg（0.05/30079）=5.78，會導(dǎo)致GWAS分析結(jié)果信號變少，一些羊毛性狀的基因位點(diǎn)被篩除，因此，本試驗(yàn)利用GEC軟件確定閾值線，避免關(guān)聯(lián)分析中出現(xiàn)假陰性結(jié)果。

本研究根據(jù)GWAS分析結(jié)果發(fā)現(xiàn)了4個SNPs位點(diǎn)與羊毛纖維直徑性狀相關(guān)；2個SNPs與自然長度相關(guān)；4個SNPs與伸直長度相關(guān)；3個SNPs與伸直率相關(guān)，并根據(jù)綿羊基因組在這些相關(guān)SNPs區(qū)域附近得到了TTC14、PEX5L、RAB38等39個候選基因（表4）。關(guān)于TTC14基因，Wang等［36］發(fā)現(xiàn)在肺腺癌組織中TTC14基因的信使RNA水平增加，還有研究表明TTC14基因參與外胚層分化超級通路，而外胚層分化中有一部分組織會分化形成皮膚等器官［37］。以往研究表明，皮膚是羊毛發(fā)生和生長的場所，由此，我們推測TTC14對羊毛性狀具有潛在影響。PEX5L基因參與蛋白質(zhì)輸入過氧化物酶體基質(zhì)、對接和cAMP介導(dǎo)的信號傳導(dǎo)的調(diào)節(jié)，在維持蛋白質(zhì)位置和調(diào)節(jié)膜電位的上游或內(nèi)部起作用。位于胞質(zhì)溶膠中，是受體復(fù)合物的一部分［38］。參與醚脂質(zhì)生物合成通路［39］，而醚脂類代表了一類新的、有前途的化合物，其特征是能夠選擇性地調(diào)節(jié)某些膜蛋白的活性，通過原始的作用機(jī)制具有抗增殖特性［40］，因此推測其在機(jī)體免疫方面具有延緩癌細(xì)胞等相關(guān)細(xì)胞增殖的作用。此外，還有研究表明，SORCS3是一種神經(jīng)元特異性跨膜蛋白，參與細(xì)胞內(nèi)囊泡和質(zhì)膜之間蛋白質(zhì)的運(yùn)輸，與多種神經(jīng)精神疾病和行為表型有關(guān)［41］，具有抑制腫瘤作用［42］。RAB38是一種GTP酶蛋白，與調(diào)節(jié)腫瘤細(xì)胞增殖和存活有關(guān)［43］。LAMP3在干燥綜合征（一種慢性自身免疫疾?。┱{(diào)控機(jī)制中發(fā)揮作用，具有誘導(dǎo)細(xì)胞凋亡的作用［44］。在小鼠中LAMP3在調(diào)節(jié)肺表面活性物質(zhì)穩(wěn)態(tài)和正常肺功能方面起著關(guān)鍵作用［45］。還有研究表明，LAMP3通?？赡芘c腫瘤細(xì)胞遷移到周圍淋巴管和免疫功能有關(guān)［46-47］。B3GNT5過表達(dá)和糖基化與乳腺癌中增強(qiáng)的癌癥干細(xì)胞特性關(guān)聯(lián)［48］。許多研究調(diào)查了RPTOR在致癌中的作用［49-50］以及RPTOR甲基化與癌癥之間的關(guān)聯(lián)［51-52］。

有研究表明，ATG5作為一種自噬因子，可作為惡性間皮瘤早期檢測的生物標(biāo)志物［53］。PRDM1基因（編碼BLIMP1轉(zhuǎn)錄因子）參與彌漫性大B細(xì)胞淋巴瘤［54-55］。在腫瘤中PRDM1/BLIMP1過表達(dá)是調(diào)節(jié)腫瘤生長的一把雙刃劍。PRDM1/BLIMP1過表達(dá)抑制細(xì)胞內(nèi)源性細(xì)胞生長，同時通過上調(diào)PD-L1和抑制CD8T細(xì)胞抗腫瘤免疫反應(yīng)來促進(jìn)腫瘤細(xì)胞免疫逃避［56］。此外，還有研究表明，ATG5和CRYBG1兩個基因是腫瘤抑制基因，可以阻止腫瘤細(xì)胞增殖［57］。

研究發(fā)現(xiàn)，MSX1能夠抑制Wnt/β-catenin信號通路，并且調(diào)節(jié)Wnt/β-catenin信號通路的能力，對于MSX1抑制膠質(zhì)母細(xì)胞瘤細(xì)胞遷移和侵襲至關(guān)重要［58］。ERAP2是一種參與免疫抗原呈遞的氨肽酶，參與抗原呈遞過程，其作為全長蛋白的表達(dá)在幾種免疫介導(dǎo)的疾病中具有保護(hù)作用［59］。LNPEP是一種同源基因，是一種專門處理細(xì)胞外抗原以呈遞給T細(xì)胞的氨肽酶［60］。有研究發(fā)現(xiàn)，ERAP1和ERAP2很可能都通過基因復(fù)制起源于LNPEP［61］。ZBTB49作為miR-21-3p的靶點(diǎn)基因在TGF-β和Wnt信號通路中發(fā)揮作用［62］。Tu等［63］確定了NSG1在食管鱗狀細(xì)胞癌中的致癌作用，并顯示NSG1高表達(dá)與食管鱗狀細(xì)胞癌患者預(yù)后不良之間存在相關(guān)性。NUBPL可通過誘導(dǎo)上皮-間充質(zhì)轉(zhuǎn)化（癌細(xì)胞轉(zhuǎn)移的主要機(jī)制）和激活細(xì)胞外調(diào)節(jié)蛋白激酶對結(jié)直腸癌產(chǎn)生重要影響［64］。

綜上所述，并結(jié)合本研究基因富集的結(jié)果，我們假設(shè)上述所有與免疫相關(guān)的基因是通過影響毛囊免疫赦免相關(guān)機(jī)制調(diào)控，影響了毛囊發(fā)育和羊毛的生長，進(jìn)而改變了羊毛性狀，預(yù)測結(jié)果如圖8所示，根據(jù)研究結(jié)果和相關(guān)基因研究發(fā)現(xiàn)PEX5L、ERAP2、LNPEP、SORCS3、ATG5、PRDM1、CRYBG1、LAMP3、MCF2L2、B3GNT5、PRTOR、NUBPL、RAB38等13個基因都與癌癥腫瘤相關(guān)。腫瘤癌癥作為免疫逃逸也不會引起免疫系統(tǒng)產(chǎn)生免疫反應(yīng)，這一現(xiàn)象與免疫赦免機(jī)制類似，所以我們推測這13個基因?qū)γ颐庖呱饷饩哂袧撛诘恼{(diào)控作用。其中，PEX5L、ERAP2、LNPEP、SORCS3、ATG5、PRDM1、CRYBG1等7個基因的表達(dá)對腫瘤癌癥具有抑制作用，包括抑制癌癥腫瘤細(xì)胞增殖，轉(zhuǎn)移等，我們推測其對毛囊免疫赦免有潛在負(fù)向調(diào)控。RAB38、LAMP3、MCF2L2、B3GNT5、PRTOR、NUBPL等6個基因表達(dá)均具有致癌作用，誘導(dǎo)癌癥細(xì)胞增殖分化，推測其對毛囊免疫赦免正向調(diào)控，此外，在本研究結(jié)果中還發(fā)現(xiàn)了可以直接參與毛囊發(fā)育和羊毛發(fā)生的基因如MSX1（參與毛囊發(fā)育：細(xì)胞分化和Wnt通路）、TTC14（外胚層分化）、ZBTB49（參與Wnt通路）。還有一些基因并未直接參與免疫相關(guān)機(jī)制，但在細(xì)胞分化、膜運(yùn)輸、能量代謝等一些信號通路中發(fā)揮作用，如STX18［65］、SUCLG1［66］、MCCC1［66］等，推測其可能間接影響毛囊免疫赦免（包括信號傳導(dǎo)、物質(zhì)合成與運(yùn)輸?shù)确矫妫?，上述發(fā)現(xiàn)可為我們后續(xù)研究羊毛性狀主效基因以及尋找遺傳標(biāo)記提供重要參考。

4 結(jié) 論

GWAS分析在全基因組上得到了4個SNPs位點(diǎn)可能影響羊毛經(jīng)濟(jì)性狀，分別位于1號、6號及8號染色體上，9個在染色體水平上顯著相關(guān)的SNPs位點(diǎn)可能對羊毛性狀具有潛在意義，分別位于3、5、8、11、18、21、22、25號染色體。依據(jù)GWAS分析結(jié)果尋找到39個可能潛在影響羊毛性狀的候選基因。本研究可為后續(xù)探究羊毛性狀的遺傳機(jī)制和尋找羊毛相關(guān)分子標(biāo)記以及候選基因的研究提供重要遺傳參考。

參考文獻(xiàn)（References）：

［1］ 李 莉，榮威恒，白俊艷，等.影響敖漢細(xì)毛羊早期主要性狀的非遺傳因素分析［J］.畜牧與飼料科學(xué)，2006，27（2）：20-23.

LI L，RONG W H，BAI J Y，et al.Analysis of non-genetic factors influencing early main traits in Aohan merino sheep［J］.Animal Husbandry and Feed Science，2006，27（2）：20-23.（in Chinese）

［2］ 梅 花，榮威恒.敖漢細(xì)毛羊主要數(shù)量性狀遺傳力的估計［J］.中國草食動物，2011，31（4）：42-43.

MEI H，RONG W H.Estimation of heritability of main quantitative traits in Aohan fine wool sheep［J］.China Herbivores， 2011，31（4）： 42-43.（in Chinese）

［3］ 柳 楠，王春亮，賀建寧，等.敖漢細(xì)毛羊不同部位皮膚毛囊發(fā)育及形態(tài)結(jié)構(gòu)研究［J］.中國畜牧雜志，2015，51（17）：1-5.

LIU N，WANG C L，HE J N，et al.Study on skin hair follicle development and morphology of different body parts of Aohan fine wool sheep［J］.Chinese Journal of Animal Science，2015，51（17）：1-5.（in Chinese）

［4］ WITTENBURG D，BONK S，DOSCHORIS M，et al.Design of experiments for fine-mapping quantitative trait loci in livestock populations［J］.BMC Genet，2020，21（1）：66.

［5］ SHIROKOVA V，BIGGS L C，JUSSILA M，et al.Foxi3 deficiency compromises hair follicle stem cell specification and activation［J］.Stem Cells，2016，34（7）：1896-1908.

［6］ BOTCHKAREV V A，SHAROV A A.BMP signaling in the control of skin development and hair follicle growth［J］. Differentiation， 2004，72（9/10）：512-526.

［7］ BOLORMAA S，SWAN A A，BROWN D J，et al.Multiple-trait QTL mapping and genomic prediction for wool traits in sheep［J］. Genet Sel Evol， 2017，49（1）：62.

［8］ WANG C，YUAN Z H，HU R X，et al.Association of SNPs within PTPN3 gene with wool production and growth traits in a dual-purpose sheep population［J］.Anim Biotechnol，2023，34（4）：1429-1435.

［9］ YUE L，LU Z K，GUO T T，et al.Association of SLIT3 and ZNF280B gene polymorphisms with wool fiber diameter［J］.Animals， 2023，13（22）：3552.

［10］ SALLAM A M，GAD-ALLAH A A，AL-BITAR E M.Association analysis of the ovine KAP6-1 gene and wool traits in Barki sheep［J］.Anim Biotechnol，2021，32（6）：733-739.

［11］ SALLAM A M，GAD-ALLAH A A，ALBETAR E M.Genetic variation in the ovine KAP22-1 gene and its effect on wool traits in Egyptian sheep［J］.Arch Anim Breed，2022，65（3）：293-300.

［12］ MA G W，CHU Y K，ZHANG W J，et al.Polymorphisms of FST gene and their association with wool quality traits in Chinese Merino sheep［J］.PLoS One，2017，12（4）：e0174868.

［13］ DARWISH H R，EL-SHORBAGY H M，ABOU-EISHA A M，et al.New polymorphism in the 5′ flanking region of IGF-1 gene and its association with wool traits in Egyptian Barki sheep［J］.J Genet Eng Biotechnol，2017，15（2）：437-441.

［14］ K?CHL S，NIEDERST?TTER H，PARSON W.DNA extraction and quantitation of forensic samples using the phenol-chloroform method and real-time PCR［J］.Methods Mol Biol，2005，297：13-30.

［15］ 劉書東.中國美利奴羊（新疆型）毛品質(zhì)性狀全基因組關(guān)聯(lián)分析［D］.石河子：石河子大學(xué)，2017.

LIU S D.Genome-wide association study for wool production traits in Chinese Merino sheep （Xinjiang type）［D］.Shihezi：Shihezi University，2017.（in Chinese）

［16］ CHANG C C，CHOW C C，TELLIER L C A M，et al.Second-generation PLINK：rising to the challenge of larger and richer datasets［J］.GigaScience，2015，4（1）：7.

［17］ TUERSUNTUOHETI M，ZHANG J H，ZHOU W，et al.Exploring the growth trait molecular markers in two sheep breeds based on Genome-wide association analysis［J］.PLoS One，2023，18（3）：e0283383.

［18］ YANG J，LEE S H，GODDARD M E，et al.GCTA：a tool for genome-wide complex trait analysis［J］.Am J Hum Genet，2011， 88（1）：76-82.

［19］ LI L，LI Y F，MA Q，et al.Analysis of family structure and paternity test of tan sheep in Yanchi Area，China［J］.Animals， 2022，12（22）：3099.

［20］ ZHANG C，DONG S S，XU J Y，et al.PopLDdecay：a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files［J］.Bioinformatics，2019，35（10）：1786-1788.

［21］ YU J M，PRESSOIR G，BRIGGS W H，et al.A unified mixed-model method for association mapping that accounts for multiple levels of relatedness［J］.Nat Genet，2006，38（2）：203-208.

［22］ ZHOU X，STEPHENS M.Genome-wide efficient mixed-model analysis for association studies［J］.Nat Genet，2012，44（7）： 821-824.

［23］ PARK M N，CHOI J A，LEE K T，et al.Genome-wide association study of chicken plumage pigmentation［J］.Asian Australas J Anim Sci，2013，26（11）：1523-1528.

［24］ TURNER S D.Qqman：an R package for visualizing GWAS results using Q-Q and manhattan plots［J］. J Open Source Softw， 2018， 3（25）：731.

［25］ LI M X，YEUNG J M Y，CHERNY S S，et al.Evaluating the effective numbers of independent tests and significant p-value thresholds in commercial genotyping arrays and public imputation reference datasets［J］.Hum Genet，2012，131（5）：747-756.

［26］ SAHANA G，GULDBRANDTSEN B，BENDIXEN C，et al.Genome-wide association mapping for female fertility traits in Danish and Swedish Holstein cattle［J］.Anim Genet，2010，41（6）：579-588.

［27］ WANG Z P，ZHANG H，YANG H，et al.Genome-wide association study for wool production traits in a Chinese Merino sheep population［J］.PLoS One，2014，9（9）：e107101.

［28］ RONG E G，YANG H，ZHANG Z W，et al.Association of methionine synthase gene polymorphisms with wool production and quality traits in Chinese Merino population［J］.J Anim Sci，2015，93（10）：4601-4609.

［29］ MU F，RONG E G，JING Y，et al.Structural characterization and association of ovine Dickkopf-1 gene with wool production and quality traits in Chinese Merino［J］.Genes （Basel），2017，8（12）：400.

［30］ LI W R，LIU C X，ZHANG X M，et al.CRISPR/Cas9-mediated loss of FGF5 function increases wool staple length in sheep［J］.FEBS J，2017，284（17）：2764-2773.

［31］ ZHAO H C，GUO T T，LU Z K，et al.Genome-wide association studies detects candidate genes for wool traits by re-sequencing in Chinese fine-wool sheep［J］.BMC Genomics，2021，22：127.

［32］ BECKER G M，WOODS J L，SCHAUER C S，et al.Genetic association of wool quality characteristics in United States Rambouillet sheep［J］.Front Genet，2023，13：1081175.

［33］ MARCHINI J，CARDON L R，PHILLIPS M S，et al.The effects of human population structure on large genetic association studies［J］.Nat Genet，2004，36（5）：512-517.

［34］ VANRADEN P M.Efficient methods to compute genomic predictions［J］.J Dairy Sci，2008，91（11）：4414-4423.

［35］ ARMSTRONG R A.When to use the Bonferroni correction［J］.Ophthalmic Physiologic Optic，2014，34（5）：502-508.

［36］ WANG J，SHIVAKUMAR S，BARKER K，et al.Comparative study of autoantibody responses between lung adenocarcinoma and benign pulmonary nodules［J］.J Thorac Oncol，2016，11（3）：334-345.

［37］ SALOMONIS N，WILLIGHAGEN E，ROUDBARI Z，et al.Ectoderm differentiation （WP2858）［EB/OL］.（2021-05-12）.https：∥ pubchem.ncbi.nlm.nih.gov/pathway/WikiPathways：WP2858.

［38］ GHOSH M，DENKERT N，REUTER M，et al.Dynamics of the translocation pore of the human peroxisomal protein import machinery［J］.Biol Chem，2023，404（2/3）：169-178.

［39］ LIPIDS C，SLENTER D，WILLIGHAGEN E，et al.Ether lipid biosynthesis （WP5275）［EB/OL］.（2024-01-30）.https：∥www. wikipathways. org/pathways/WP5275.html.

［40］ GOMES M A G B，BAUDUIN A，LE ROUX C，et al.Synthesis of ether lipids：natural compounds and analogues［J］.Beilstein J Org Chem，2023，19（1）：1299-1369.

［41］ KAMRAN M，LAIGHNEACH A，BIBI F，et al.Independent associated SNPs at SORCS3 and its protein interactors for multiple brain-related disorders and traits［J］.Genes，2023，14（2）：482.

［42］ ZHANG Y Q，LI Y，F(xiàn)AN Y H，et al.SorCS3 promotes the internalization of p75NTR to inhibit GBM progression［J］.Cell Death Dis，2022，13（4）：313.

［43］ BIANCHETTI E，BATES S J，NGUYEN T T T，et al.RAB38 facilitates energy metabolism and counteracts cell death in glioblastoma cells［J］.Cells，2021，10（7）：1643.

［44］ TANAKA T，WARNER B M，MICHAEL D G，et al.LAMP3 inhibits autophagy and contributes to cell death by lysosomal membrane permeabilization［J］.Autophagy，2022，18（7）：1629-1647.

［45］ LUNDING L P，KRAUSE D，STICHTENOTH G，et al.LAMP3 deficiency affects surfactant homeostasis in mice［J］.PLoS Genet，2021，17（6）：e1009619.

［46］ KANAO H，ENOMOTO T，KIMURA T，et al.Overexpression of LAMP3/TSC403/DC-LAMP promotes metastasis in uterine cervical cancer［J］.Cancer Res，2005，65（19）：8640-8645.

［47］ HARALAMBIEVA I H，OBERG A L，DHIMAN N，et al.High-dimensional gene expression profiling studies in high and low responders to primary smallpox vaccination［J］.J Infect Dis，2012，206（10）：1512-1520.

［48］ MIAO Z R，CAO Q H，LIAO R C，et al.Elevated transcription and glycosylation of B3GNT5 promotes breast cancer aggressiveness［J］.J Exp Clin Cancer Res，2022，41（1）：169.

［49］ SHARLOW E R，LEIMGRUBER S，LIRA A，et al.A small molecule screen exposes mTOR signaling pathway involvement in radiation-induced apoptosis［J］.ACS Chem Biol，2016，11（5）：1428-1437.

［50］ SAXTON R A，SABATINI D M.mTOR signaling in growth，metabolism，and disease［J］.Cell，2017，168（6）：960-976.

［51］ SHIRKAVAND A，BOROUJENI Z N，ALEYASIN S A.Examination of methylation changes of VIM，CXCR4，DOK7，and SPDEF genes in peripheral blood DNA in breast cancer patients［J］.Indian J Cancer，2018，55（4）：366-371.

［52］ SONG L L，SHEN L J，LI H，et al.Age at natural menopause and hypertension among middle-aged and older Chinese women［J］.J Hypertens，2018，36（3）：594-600.

［53］ TOMASETTI M，MONACO F，STROGOVETS O，et al.ATG5 as biomarker for early detection of malignant mesothelioma［J］. BMC Res Notes，2023，16（1）：61.

［54］ PASQUALUCCI L，COMPAGNO M，HOULDSWORTH J，et al.Inactivation of the PRDM1/BLIMP1 gene in diffuse large B cell lymphoma［J］.J Exp Med，2006，203（2）：311-317.

［55］ LIU Y Y，LEBOEUF C，SHI J Y，et al.Rituximab plus CHOP （R-CHOP） overcomes PRDM1-associated resistance to chemotherapy in patients with diffuse large B-cell lymphoma［J］.Blood，2007，110（1）：339-344.

［56］ LI Q，ZHANG L R，YOU W H，et al.PRDM1/BLIMP1 induces cancer immune evasion by modulating the USP22-SPI1-PD-L1 axis in hepatocellular carcinoma cells［J］.Nat Commun，2022，13（1）：7677.

［57］ HUANG Y，DE REYNIèS A，DE LEVAL L，et al.Gene expression profiling identifies emerging oncogenic pathways operating in extranodal NK/T-cell lymphoma，"" nasal type［J］.Blood，2010，115（6）：1226-1237.

［58］ TAO H Q，GUO L，CHEN L F，et al.MSX1 inhibits cell migration and invasion through regulating the Wnt/β-catenin pathway in glioblastoma［J］.Tumor Biol，2016，37（1）：1097-1104.

［59］ MATTORRE B，TEDESCHI V，PALDINO G，et al.The emerging multifunctional roles of ERAP1，ERAP2 and IRAP between antigen processing and renin-angiotensin system modulation［J］.Front Immunol，2022，13：1002375.

［60］ MARUSINA A I，JI-XU A，LE S T，et al.Cell-specific and variant-linked alterations in expression of ERAP1，ERAP2，and LNPEP aminopeptidases in psoriasis［J］.J Invest Dermatol，2023，143（7）：1157-1167.e10.

［61］ PALADINI F，F(xiàn)IORILLO M T，TEDESCHI V，et al.The multifaceted nature of aminopeptidases ERAP1，ERAP2，and LNPEP：from evolution to disease［J］.Front Immunol，2020，11：1576.

［62］ ONG J，TIMENS W，RAJENDRAN V，et al.Identification of transforming growth factor-beta-regulated microRNAs and the microRNA-targetomes in primary lung fibroblasts［J］.PLoS One，2017，12（9）：e0183815.

［63］ TU M S，YIN X Q，ZHUANG W Z，et al.NSG1 promotes glycolytic metabolism to enhance Esophageal squamous cell carcinoma EMT process by upregulating TGF-β［J］.Cell Death Discov，2023，9（1）：391.

［64］ WANG Y H，WU N，SUN D L，et al.NUBPL， a novel metastasis-related gene，promotes colorectal carcinoma cell motility by inducing epithelial-mesenchymal transition［J］.Cancer Sci，2017，108（6）：1169-1176.

［65］ GUILLEMYN B，DE SAFFEL H，BEK J W，et al.Syntaxin 18 defects in human and zebrafish unravel key roles in early cartilage and bone development［J］.J Bone Miner Res，2023，38（11）：1718-1730.

［66］ HANSPERS K，WILLIGHAGEN E，WEITZ E.Amino acid metabolism （WP3925）［EB/OL］.（2021-05-16）. https：∥pubchem. ncbi.nlm.nih.gov/pathway/WikiPathways：WP3925.

（編輯 郭云雁）