基于轉(zhuǎn)錄組測(cè)序的奶牛產(chǎn)奶性狀重要功能基因挖掘

陳一丹 張昱 楊潔 張勤，2 姜力

（1. 中國(guó)農(nóng)業(yè)大學(xué)動(dòng)物科學(xué)技術(shù)學(xué)院 畜禽育種國(guó)家工程實(shí)驗(yàn)室 農(nóng)業(yè)農(nóng)村部動(dòng)物遺傳育種與繁殖重點(diǎn)實(shí)驗(yàn)室，北京 100193；2. 山東農(nóng)業(yè)大學(xué)動(dòng)物科技學(xué)院，泰安 271018）

牛奶營(yíng)養(yǎng)豐富，是人類(lèi)良好的蛋白質(zhì)和鈣的來(lái)源，同時(shí)也是生產(chǎn)多種加工食品的原料，對(duì)人類(lèi)生活具有重要意義［1］。產(chǎn)奶性狀是奶牛最重要的生產(chǎn)性狀，挖掘影響奶牛產(chǎn)奶性狀的功能基因或遺傳變異一直是奶牛分子育種中的研究熱點(diǎn)。如何準(zhǔn)確的確定基因型與表型之間的關(guān)系是當(dāng)前畜禽遺傳研究的重要挑戰(zhàn)之一。而基因組變異與基因表達(dá)以及表型關(guān)系的整合分析是揭示基因型與表型關(guān)系的有效途徑［2］。

近年來(lái)，快速發(fā)展的二代測(cè)序技術(shù)為揭示復(fù)雜性狀遺傳基礎(chǔ)以及表型變異的機(jī)制提供了重要工具。RNA-seq技術(shù)不但可以檢測(cè)細(xì)胞或組織中所有基因的表達(dá)，而且為鑒定轉(zhuǎn)錄本中的遺傳變異提供了新的機(jī)會(huì)［3］，為深入挖掘奶牛產(chǎn)奶性狀功能基因和重要突變位點(diǎn)提供了便利條件。眾所周知，轉(zhuǎn)錄本中的遺傳變異會(huì)對(duì)基因表達(dá)和基因的轉(zhuǎn)錄后調(diào)控具有調(diào)控作用。因此，本研究以中國(guó)荷斯坦奶牛為研究對(duì)象，選擇高產(chǎn)奶牛和低產(chǎn)奶牛泌乳期的血液組織進(jìn)行轉(zhuǎn)錄組分析，通過(guò)檢測(cè)差異表達(dá)基因和轉(zhuǎn)錄本中的遺傳變異位點(diǎn)，并結(jié)合生物信息學(xué)分析，篩選與產(chǎn)奶性狀相關(guān)的重要功能基因和遺傳變異。本研究旨在進(jìn)一步挖掘和鑒定影響奶牛產(chǎn)奶性狀的遺傳標(biāo)記位點(diǎn)，為提高我國(guó)奶牛分子選育的準(zhǔn)確性提供重要信息。

1 材料與方法

1.1 材料

本研究試驗(yàn)個(gè)體來(lái)自于北京三元金銀島牛場(chǎng)，按照305 d產(chǎn)奶量、乳脂量和乳蛋白量3個(gè)性狀表型信息及育種值信息對(duì)所有個(gè)體進(jìn)行篩選，最終選取3頭高產(chǎn)荷斯坦奶牛和3頭低產(chǎn)荷斯坦奶牛，產(chǎn)奶性狀表型信息如表1所示。所有試驗(yàn)?zāi)膛＞鶠榻】祩€(gè)體，在同一飼養(yǎng)環(huán)境下飼養(yǎng)。試驗(yàn)牛進(jìn)行尾椎靜脈采血，血液收集到抗凝管后立即與RNA保護(hù)試劑混合，用于RNA的提取。

表1 試驗(yàn)群體

1.2 方法

1.2.1 RNA提取及測(cè)序 使用Trizol法提取所有樣品的總RNA。利用Nanodrop核酸分析儀和Qubit對(duì)RNA進(jìn)行純度和濃度的檢測(cè)。使用Agilent 2100對(duì)RNA的完整性進(jìn)行檢測(cè)。每個(gè)樣品的濃度均大于100 ng/μL，RIN值＞7.5。cDNA 文庫(kù)的構(gòu)建參照Illumina TruSeqTMRNA（Illumia，USA）樣品制備試劑盒操作說(shuō)明進(jìn)行。cDNA文庫(kù)質(zhì)檢合格后，使用Illumina Hiseq 2500測(cè)序平臺(tái)進(jìn)行的雙末端測(cè)序。

1.2.2 測(cè)序數(shù)據(jù)質(zhì)量控制 測(cè)序后的原始序列（Raw reads）需去除帶接頭、含有 poly-N以及低質(zhì)量reads形成Clean reads。將Clean reads比對(duì)到牛的參考基因組序列（UMD3.1）和相應(yīng)的基因注釋文件（UMD3.1）。

1.2.3 遺傳變異檢測(cè) 使用BWA［4］軟件將轉(zhuǎn)錄組雙端測(cè)序數(shù)據(jù)比對(duì)到參考基因組上，然后使用SAMtools［5］與BCFtools［6］軟件進(jìn)行遺傳變異的檢測(cè)。使用vcffilter軟件以DP＞5、QUAL＞30、MQ＞40為條件對(duì)檢測(cè)到的變異進(jìn)行過(guò)濾。使用BCFtools注釋已存在于dbSNP數(shù)據(jù)庫(kù)中的遺傳變異。使用snpEff［7］軟件對(duì)所有變異進(jìn)行功能性預(yù)測(cè)分析。

1.2.4 遺傳變異篩選及功能基因注釋 基于檢測(cè)到的遺傳變異，篩選出突變等位基因頻率在高產(chǎn)組中≥1/2，且低產(chǎn)組中≤1/3的遺傳變異；以及在低產(chǎn)組中突變等位基因頻率≥1/2，且高產(chǎn)組中≤1/3的遺傳變異。在這些遺傳變異中，進(jìn)一步篩選出被snpEff預(yù)測(cè)為影響程度高（HIGH）和中等（MODERATE）的變異，并統(tǒng)計(jì)這些變異所在的基因。

1.2.5 基因的差異表達(dá)分析 利用 BWA軟件進(jìn)行轉(zhuǎn)錄組測(cè)序數(shù)據(jù)的比對(duì)。使用featureCounts［8］統(tǒng)計(jì)counts數(shù)后，利用DESeq2［9］軟件進(jìn)行基因差異表達(dá)分析。以P＜0.05、差異表達(dá)倍數(shù)＞1.5為顯著性閾值篩選高低產(chǎn)組之間的差異表達(dá)基因。

1.2.6 基因的富集分析 利用DAVID（https：//david.ncifcrf.gov/）數(shù)據(jù)庫(kù)對(duì)所有差異表達(dá)基因進(jìn)行GO功能注釋和KEGG通路分析，以P＜0.05為顯著性閾值篩選顯著的GO條目和富集的通路。

1.2.7 候選功能基因的篩選 將篩選的重要遺傳變異涉及的基因與差異表達(dá)基因取交集，篩選潛在影響基因表達(dá)的遺傳變異。利用DESeq2軟件輸出的標(biāo)準(zhǔn)化表達(dá)矩陣，使用R語(yǔ)言對(duì)篩選的候選基因進(jìn)行繪圖，并對(duì)這些基因進(jìn)行生物信息學(xué)分析。同時(shí)將這些基因與Animal QTL數(shù)據(jù)庫(kù)（https：//www.animalgenome.org/cgi-bin/QTLdb/BT/index）中包含的產(chǎn)奶性狀QTLs進(jìn)行比對(duì)，篩選與產(chǎn)奶量、乳脂量、乳蛋白量、乳脂率、乳蛋白率QTLs重疊的候選功能基因，進(jìn)一步挖掘與奶牛產(chǎn)奶性狀相關(guān)的功能基因。

2 結(jié)果

2.1 測(cè)序數(shù)據(jù)質(zhì)量檢測(cè)結(jié)果

原始測(cè)序數(shù)據(jù)經(jīng)質(zhì)控后最終獲得38.1G的Clean reads。將這些Clean reads比對(duì)到牛的參考基因組上，每個(gè)個(gè)體的比對(duì)率在 83.6%-92.8%之間，結(jié)果如表2所示。這些reads比對(duì)到mRNA中的比例在84.5%-94.7%之間。其中大部分reads都比對(duì)到基因的編碼區(qū)，比對(duì)率在62.7%-74.9%之間。

表2 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估

2.2 遺傳變異檢測(cè)結(jié)果

測(cè)序數(shù)據(jù)經(jīng)分析后獲得6個(gè)個(gè)體轉(zhuǎn)錄本中的所有遺傳變異，統(tǒng)計(jì)結(jié)果如表3所示。所有個(gè)體轉(zhuǎn)錄本檢測(cè)到的遺傳變異平均數(shù)為145 750，其中SNP的數(shù)量占總遺傳變異數(shù)的91.23%。與dbSNP數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果顯示已被數(shù)據(jù)庫(kù)收錄的遺傳變異比例在80.95%到85.28%之間，其余為本研究新發(fā)現(xiàn)的遺傳變異。

表3 遺傳變異數(shù)目統(tǒng)計(jì)

使用snpEff軟件對(duì)所有遺傳變異進(jìn)行功能性預(yù)測(cè)分析，將其影響程度分為HIGH、MODERATE、LOW和MODIFIER四種情況，結(jié)果如表4所示。被預(yù)測(cè)為HIGH的遺傳變異對(duì)蛋白質(zhì)產(chǎn)物有高度的破壞性，可能導(dǎo)致蛋白質(zhì)功能喪失或截?cái)?；被預(yù)測(cè)為MODERATE的遺傳變異可能對(duì)蛋白質(zhì)的性質(zhì)和功能具有影響；被預(yù)測(cè)為L(zhǎng)OW的遺傳變異不會(huì)使蛋白質(zhì)發(fā)生改變；而被預(yù)測(cè)為MODIFIER的遺傳變異是非編碼的遺傳變異或影響非編碼基因的變異［3］。將這些遺傳變異所在的基因組功能區(qū)域進(jìn)行注釋?zhuān)ū?），結(jié)果顯示每個(gè)個(gè)體在基因的外顯子區(qū)檢測(cè)到的遺傳變異數(shù)平均為15 642個(gè)，在3' UTR和5' UTR檢測(cè)到的變異數(shù)量平均為7 883個(gè)和942個(gè)。一些基因內(nèi)含子上的遺傳變異被預(yù)測(cè)可能導(dǎo)致新可變剪切體的產(chǎn)生。

表4 遺傳變異影響程度的預(yù)測(cè)結(jié)果

表5 遺傳變異在基因組功能區(qū)的分布

2.3 遺傳變異篩選結(jié)果

根據(jù)突變位點(diǎn)在高低組中不同的等位基因頻率，對(duì)所有遺傳變異進(jìn)行篩選。結(jié)果顯示在高產(chǎn)組中突變位點(diǎn)等位基因頻率≥1/2，且在低產(chǎn)組中等位基因頻率≤1/3的遺傳變異共有27 143個(gè)，其中22 301個(gè)已存在于dbSNP數(shù)據(jù)庫(kù)中（圖1），其中被預(yù)測(cè)為對(duì)蛋白功能有重要（HIGH）影響的遺傳變異有330個(gè)，對(duì)蛋白功能有中等程度（MODERATE）影響的遺傳變異有1 174個(gè)（圖1），這些變異共涉及979個(gè)基因（圖2-A）。

同時(shí)，我們篩選到在低產(chǎn)組中突變等位基因頻率≥1/2，且在高產(chǎn)組中等位基因頻率≤1/3的遺傳變異共40 650個(gè)，其中33 687個(gè)已存在于dbSNP數(shù)據(jù)庫(kù)中（圖1），其中被預(yù)測(cè)對(duì)蛋白功能具有重要（HIGH）影響的遺傳變異有120個(gè)，對(duì)蛋白功能為中等程度（MODERATE）影響的遺傳變異有927個(gè)（圖1），這些變異共涉及643個(gè)基因（圖2-B）。

圖1 遺傳變異篩選結(jié)果分析

2.4 差異表達(dá)基因分析結(jié)果

利用DESeq2軟件進(jìn)行基因的差異表達(dá)分析，以P＜0.05、|Fold Change|≥ 1.5作為顯著性閾值，最終在高低組之間檢測(cè)到431個(gè)差異表達(dá)基因，其中190個(gè)基因在低產(chǎn)組的表達(dá)量高于在高產(chǎn)組的表達(dá)量，241個(gè)基因在低產(chǎn)組的表達(dá)量低于在高產(chǎn)組的表達(dá)量（圖2）。

2.5 產(chǎn)奶性狀重要候選基因的篩選

圖2 差異表達(dá)基因與含有重要遺傳變異基因的整合分析

將差異表達(dá)基因與上述篩選的含有重要遺傳變異的基因取交集，共獲得47個(gè)基因（圖2）。其中在高產(chǎn)組中突變等位基因頻率≥1/2，且低產(chǎn)組中等位基因頻率≤1/3，同時(shí)預(yù)測(cè)對(duì)蛋白功能具有HIGH或MODERATE影響的遺傳變異所涉及的基因與差異表達(dá)基因取交集，共篩選到28個(gè)基因，其中24個(gè)基因在高產(chǎn)組中的表達(dá)量顯著高于在低產(chǎn)組（圖2-A）。在低產(chǎn)組中等位基因頻率≥1/2，且在高產(chǎn)組中等位基因頻率≤1/3，并且預(yù)測(cè)對(duì)蛋白功能具有HIGH或MODERATE影響的遺傳變異所涉及的基因與差異表達(dá)基因取交集，共獲得19個(gè)基因，其中14個(gè)基因在低產(chǎn)組中的表達(dá)量顯著高于在高產(chǎn)組（圖2-B）。

將篩選到的47個(gè)基因進(jìn)一步進(jìn)行生物信息學(xué)分析，檢測(cè)到1個(gè)顯著的GO條目（GO：0004252～serine-type endopeptidase activity）和1條顯著的通路（bta01100：Metabolic pathways），共包含9個(gè) 基 因：ASS1、CKB、GGT1、UPP1、MGAM、SDSL、HP、LTF和MMP9（表6）。

將47個(gè)基因與Animal QTL數(shù)據(jù)庫(kù)中產(chǎn)奶性狀QTLs進(jìn)行比對(duì)，篩選到4個(gè)已報(bào)道與產(chǎn)奶性狀相關(guān)的基因，即DEFB4A、LTF、PGLYRP1、MS4A8。將上述12個(gè)重要候選基因內(nèi)的遺傳變異進(jìn)行匯總，共獲得14個(gè)遺傳變異位點(diǎn)（表7），其中8個(gè)突變會(huì)引起氨基酸的錯(cuò)義突變。使用DESeq2軟件得到的標(biāo)準(zhǔn)化表達(dá)量后，對(duì)這14個(gè)突變位點(diǎn)不同基因型個(gè)體進(jìn)行基因表達(dá)量的統(tǒng)計(jì)，結(jié)果顯示不同基因型個(gè)體之間基因的表達(dá)量存在明顯差別，如圖3所示。由于本研究試驗(yàn)群體有限，未來(lái)需在更大規(guī)模的奶牛群體中進(jìn)一步驗(yàn)證。

表6 篩選基因的GO和KEGG分析

3 討論

奶牛的產(chǎn)奶性狀一直是育種工作者最為關(guān)注的性狀之一，隨著人們生活水平的日益提高，人們對(duì)牛奶的需求量不斷地增加。準(zhǔn)確的遺傳評(píng)估無(wú)疑為加快奶牛產(chǎn)奶性狀遺傳進(jìn)展發(fā)揮了重要作用。當(dāng)前，基因組選擇是繼奶牛育種中的表型選擇、育種值選擇之后的最先進(jìn)的分子育種方法。如果可以獲得更多的與產(chǎn)奶性狀有關(guān)的分子標(biāo)記信息，并將其加入到基因組選擇當(dāng)中，就可以進(jìn)一步提高奶牛分子選育的準(zhǔn)確性。因此，許多科研工作者通過(guò)全基因組關(guān)聯(lián)分析、轉(zhuǎn)錄組、蛋白質(zhì)組等多種策略挖掘影響奶牛產(chǎn)奶性狀的功能基因和遺傳變異［10-13］。目前，利用高低產(chǎn)奶牛的乳腺和肝臟組織的轉(zhuǎn)錄組研究已有報(bào)道，這些研究著重利用測(cè)序數(shù)據(jù)進(jìn)行差異表達(dá)基因、差異表達(dá)非編碼RNA的篩選，并利用生物信息學(xué)手段對(duì)基因和非編碼RNA之間的調(diào)控關(guān)系進(jìn)行預(yù)測(cè)［12，14］，從而進(jìn)一步挖掘影響奶牛產(chǎn)奶性狀的候選基因。本研究利用二代轉(zhuǎn)錄組測(cè)序技術(shù)對(duì)高低產(chǎn)組奶牛血液組織中的差異表達(dá)基因進(jìn)行檢測(cè)。鑒于奶牛在泌乳過(guò)程中許多重要的營(yíng)養(yǎng)物質(zhì)是通過(guò)血液運(yùn)送到乳腺組織合成牛奶中的各種乳成分，對(duì)血液組織的研究具有重要的意義。此外，為了實(shí)現(xiàn)分子標(biāo)記應(yīng)用于奶牛育種實(shí)踐的目標(biāo)，本研究側(cè)重于高低產(chǎn)奶牛基因組中的遺傳變異的檢測(cè)和分析，同時(shí)結(jié)合基因在高低組中的表達(dá)以及在不同基因型個(gè)體中的表達(dá)進(jìn)行重要候選基因及遺傳變異的篩選。研究結(jié)果為揭示奶牛產(chǎn)奶性狀表型差異的遺傳基礎(chǔ)及分子選育提供了重要信息。

表7 12個(gè)重要候選基因的遺傳變異情況

本研究共篩選到12個(gè)影響產(chǎn)奶性狀的候選功能基因，其中ASS1、DEFB4A、GGT1、HP、LTF、MMP9、MGAM、UPP1、PGLYRP1基因內(nèi)含有在高產(chǎn)組中突變等位基因頻率較高，而在低產(chǎn)組中等位基因頻率較低的遺傳變異。這9個(gè)重要候選基因中包含11個(gè)重要突變位點(diǎn)，其中7個(gè)為錯(cuò)義突變。位于PGLYRP1基因的SNP為同義突變，Wang等［15］的研究顯示該突變位點(diǎn)與荷斯坦奶牛的305 d產(chǎn)奶量和體細(xì)胞數(shù)密切相關(guān)，并且突變型奶牛的305 d產(chǎn)奶量相比野生型個(gè)體的產(chǎn)奶量顯著上升，體細(xì)胞評(píng)分顯著下降，推斷該突變可能通過(guò)改變mRNA的穩(wěn)定性影響蛋白質(zhì)的表達(dá)和功能，同時(shí)也可能通過(guò)密碼子的偏好性改變蛋白質(zhì)的折疊結(jié)構(gòu)，從而影響蛋白功能的發(fā)揮。研究發(fā)現(xiàn)這些基因中包含多個(gè)重要的酶類(lèi)，與營(yíng)養(yǎng)物質(zhì)的吸收、合成及代謝過(guò)程緊密相關(guān)。例如，ASS1基因所編碼的精氨酸琥珀酸合成酶與氨基酸的合成、代謝和尿素循環(huán)等都有密切聯(lián)系［16］。GGT1基因所編碼的精-谷氨酰轉(zhuǎn)移酶是一種細(xì)胞膜結(jié)合酶，在γ-谷氨酰循環(huán)中發(fā)揮重要作用，與氨基酸的吸收密切相關(guān)［17］。MGAM基因編碼的麥芽糖-葡萄糖化酶是一種位于小腸絨毛刷狀緣的酶，負(fù)責(zé)將多糖分解為葡萄糖［18］，與碳水化合物的消化和吸收、半乳糖、蔗糖和淀粉的代謝有關(guān)。此外，位于這些基因上的遺傳變異不同基因型下基因的表達(dá)量存在明顯差異，都表現(xiàn)為突變型純合子個(gè)體基因的表達(dá)量最高，野生型純合子個(gè)體的表達(dá)量最低，雜合子個(gè)體表達(dá)量普遍居于兩者之間。由于這些基因包含的遺傳變異在高產(chǎn)組中的基因頻率明顯高于低產(chǎn)組，同時(shí)顯示突變型個(gè)體基因的表達(dá)量高于野生型個(gè)體，并且這些基因在高產(chǎn)組的表達(dá)量顯著高于低產(chǎn)組。因此，本研究推測(cè)這些遺傳變異對(duì)產(chǎn)奶性狀表型有正向作用，突變型為優(yōu)勢(shì)基因型，可進(jìn)一步驗(yàn)證后作為分子標(biāo)記應(yīng)用于奶牛的分子育種中。

此外，本研究篩選到重要的一個(gè)乳成分基因，即乳鐵蛋白基因（LTF）。LTF編碼的乳鐵蛋白是牛奶中重要的營(yíng)養(yǎng)成分，在先天免疫系統(tǒng)中起著重要的作用［19-20］，參與鐵代謝、抗腫瘤、抗細(xì)菌等多種免疫調(diào)節(jié)過(guò)程［21-23］。LTF基因在許多研究中被證明與奶牛的生產(chǎn)性能有關(guān)。2010年，O' Halloran等［24］發(fā)現(xiàn)LTF上一個(gè)位于轉(zhuǎn)錄起始位點(diǎn)上游-28 bp的SNP突變（A/C）對(duì)乳蛋白量有一定影響。2015年，Mao等［25］的研究表明LTF上位于轉(zhuǎn)錄起始位點(diǎn)上游-270 bp從T到C的SNP突變對(duì)產(chǎn)奶量、乳脂率和乳蛋白率均有正向的影響，而位于-190 bp從G到A的突變則對(duì)乳脂率和乳蛋白率有不利的影響。2017年，Viale等［26］發(fā)現(xiàn)LTF中一個(gè)SNP（rs43765462）與荷斯坦奶牛的乳脂率性狀正相關(guān)。2018年，Raschia等［27］發(fā)現(xiàn)LTF上另外一個(gè)SNP（rs43706485）與奶牛的305 d產(chǎn)奶量有潛在的相關(guān)性。本研究首次在LTF基因中發(fā)現(xiàn)一個(gè)可能引起可變剪接的SNP突變（T/A）。該突變位點(diǎn)不同基因型個(gè)體LTF基因的表達(dá)量表現(xiàn)為突變純合子基因型個(gè)體（AA）遠(yuǎn)大于突變雜合子基因型個(gè)體（AT），大于野生型純合子個(gè)體（TT），并且該突變等位基因頻率在高產(chǎn)組明顯高于低產(chǎn)組。因此，本研究認(rèn)為該突變?yōu)閷?duì)奶牛的產(chǎn)奶量、乳脂量和乳蛋白量具有正向影響。由于牛奶中的天然乳鐵蛋白含量極低，每毫升牛奶中僅含0.02 mg-0.35 mg［28］。隨著乳鐵蛋白在嬰幼兒免疫調(diào)節(jié)以及抗腫瘤等多種功能作用的不斷解析［23，29-30］，目前天然乳鐵蛋白做為食品添加劑和營(yíng)養(yǎng)強(qiáng)化劑的價(jià)格仍較為昂貴。因此，培育牛奶中具有較高含量乳鐵蛋白成分的奶牛具有很高的經(jīng)濟(jì)價(jià)值。本研究發(fā)現(xiàn)的該分子標(biāo)記有望將來(lái)用于選育牛奶中富含乳鐵蛋白的奶牛品系。

同時(shí)，本研究分別在DEFB4A基因和MS4A8基因中檢測(cè)到一個(gè)重要的遺傳變異。Bagnicka 等［31］的研究結(jié)果表明DEFB4A基因上位于轉(zhuǎn)錄位點(diǎn)起始位置上游的2 239 bp從C到T的突變對(duì)奶牛的乳脂率有不利影響。本研究發(fā)現(xiàn)DEFB4A基因上存在一個(gè)從C到A的錯(cuò)義突變（rs43108924），該突變位點(diǎn)在高產(chǎn)組中的等位基因頻率高于在低產(chǎn)組中，同時(shí)該突變位點(diǎn)不同基因型個(gè)體中DEFB4A的表達(dá)量呈現(xiàn)為突變純合子基因型個(gè)體表達(dá)量大于野生型個(gè)體，因此推測(cè)該突變是對(duì)奶牛的產(chǎn)奶量、乳脂量和乳蛋白量有正向作用的有利突變。Cochran等［32］的研究發(fā)現(xiàn)MS4A8中的錯(cuò)義突變（rs109761676）對(duì)奶牛的產(chǎn)奶量、乳脂率和乳蛋白率有負(fù)面影響。而該突變位點(diǎn)同樣在本研究中被篩選到。本研究結(jié)果顯示該突變位點(diǎn)（G/T）在奶牛低產(chǎn)組中的等位基因頻率明顯高于在高產(chǎn)組中的基因頻率，再次證明了該突變對(duì)產(chǎn)奶性狀是一個(gè)潛在的不利突變。對(duì)該遺傳變異位點(diǎn)進(jìn)行不同基因型個(gè)體表達(dá)量的分析，發(fā)現(xiàn)在該基因突變型純合子個(gè)體（GG）的表達(dá)量低于雜合子基因型（GT）的個(gè)體和野生型個(gè)體（TT），說(shuō)明該突變可能是引起MS4A8表達(dá)量降低的隱性突變，通過(guò)降低基因的表達(dá)水平最終導(dǎo)致奶牛產(chǎn)奶生產(chǎn)水平的下降。

值得關(guān)注的是，CKB基因編碼的肌酸激酶B與肌酸代謝、氨基酸代謝、尿素循環(huán)等均有一定的關(guān)系［33］。SDSL基因所編碼的絲氨酸脫水酶類(lèi)似物是一種與絲氨酸脫水酶類(lèi)似的，與氨基酸的合成與代謝有關(guān)的蛋白質(zhì)［34］。本研究發(fā)現(xiàn)這兩個(gè)基因上均存在一個(gè)遺傳變異在奶牛低產(chǎn)組中突變等位基因頻率明顯高于高產(chǎn)組中等位基因頻率，且突變型純合子個(gè)體基因的表達(dá)量高于或接近于突變型雜合子表達(dá)量高于野生型表達(dá)量。這說(shuō)明一些突變也可能是通過(guò)上調(diào)所在基因表達(dá)量影響奶牛的重要代謝過(guò)程，最終對(duì)產(chǎn)奶性能產(chǎn)生不利影響。

4 結(jié)論

本研究利用奶牛產(chǎn)奶性狀表型兩尾極端個(gè)體的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)進(jìn)行遺傳變異檢測(cè)和基因的差異表達(dá)分析，旨在篩選與奶牛產(chǎn)奶性狀相關(guān)的重要候選功能基因及遺傳變異。通過(guò)對(duì)突變位點(diǎn)在高低產(chǎn)組中等位基因頻率分布的比較、基因差異表達(dá)分析、生物信息學(xué)分析以及與產(chǎn)奶性狀QTL數(shù)據(jù)庫(kù)的比對(duì)，最終篩選到12個(gè)重要的產(chǎn)奶性狀候選基因和位于這些基因內(nèi)的14個(gè)重遺傳變異。經(jīng)分析，這些變異包括潛在的有利突變和不利突變，可能通過(guò)改變基因的表達(dá)調(diào)控奶牛產(chǎn)奶性狀的表型，這些信息對(duì)將來(lái)更加準(zhǔn)確的進(jìn)行奶牛分子育種具有重要意義。