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

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

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

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

1 材料與方法

1.1 材料

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

表1 試驗群體

1.2 方法

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

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

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

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

1.2.5 基因的差異表達(dá)分析 利用 BWA軟件進(jìn)行轉(zhuǎn)錄組測序數(shù)據(jù)的比對。使用featureCounts［8］統(tǒng)計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ù)庫對所有差異表達(dá)基因進(jìn)行GO功能注釋和KEGG通路分析，以P＜0.05為顯著性閾值篩選顯著的GO條目和富集的通路。

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

2 結(jié)果

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

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

表2 測序數(shù)據(jù)質(zhì)量評估

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 討論

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

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

本研究共篩選到12個影響產(chǎn)奶性狀的候選功能基因，其中ASS1、DEFB4A、GGT1、HP、LTF、MMP9、MGAM、UPP1、PGLYRP1基因內(nèi)含有在高產(chǎn)組中突變等位基因頻率較高，而在低產(chǎn)組中等位基因頻率較低的遺傳變異。這9個重要候選基因中包含11個重要突變位點，其中7個為錯義突變。位于PGLYRP1基因的SNP為同義突變，Wang等［15］的研究顯示該突變位點與荷斯坦奶牛的305 d產(chǎn)奶量和體細(xì)胞數(shù)密切相關(guān)，并且突變型奶牛的305 d產(chǎn)奶量相比野生型個體的產(chǎn)奶量顯著上升，體細(xì)胞評分顯著下降，推斷該突變可能通過改變mRNA的穩(wěn)定性影響蛋白質(zhì)的表達(dá)和功能，同時也可能通過密碼子的偏好性改變蛋白質(zhì)的折疊結(jié)構(gòu)，從而影響蛋白功能的發(fā)揮。研究發(fā)現(xiàn)這些基因中包含多個重要的酶類，與營養(yǎng)物質(zhì)的吸收、合成及代謝過程緊密相關(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)為突變型純合子個體基因的表達(dá)量最高，野生型純合子個體的表達(dá)量最低，雜合子個體表達(dá)量普遍居于兩者之間。由于這些基因包含的遺傳變異在高產(chǎn)組中的基因頻率明顯高于低產(chǎn)組，同時顯示突變型個體基因的表達(dá)量高于野生型個體，并且這些基因在高產(chǎn)組的表達(dá)量顯著高于低產(chǎn)組。因此，本研究推測這些遺傳變異對產(chǎn)奶性狀表型有正向作用，突變型為優(yōu)勢基因型，可進(jìn)一步驗證后作為分子標(biāo)記應(yīng)用于奶牛的分子育種中。

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

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

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

4 結(jié)論

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