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

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

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

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

1 材料與方法

1.1 材料

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

表1 試驗群體

1.2 方法

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

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

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

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

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

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

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

2 結果

2.1 測序數(shù)據(jù)質量檢測結果

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

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

2.2 遺傳變異檢測結果

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

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

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

表4 遺傳變異影響程度的預測結果

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

2.3 遺傳變異篩選結果

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

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

圖1 遺傳變異篩選結果分析

2.4 差異表達基因分析結果

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

2.5 產奶性狀重要候選基因的篩選

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

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

將篩選到的47個基因進一步進行生物信息學分析，檢測到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ù)庫中產奶性狀QTLs進行比對，篩選到4個已報道與產奶性狀相關的基因，即DEFB4A、LTF、PGLYRP1、MS4A8。將上述12個重要候選基因內的遺傳變異進行匯總，共獲得14個遺傳變異位點（表7），其中8個突變會引起氨基酸的錯義突變。使用DESeq2軟件得到的標準化表達量后，對這14個突變位點不同基因型個體進行基因表達量的統(tǒng)計，結果顯示不同基因型個體之間基因的表達量存在明顯差別，如圖3所示。由于本研究試驗群體有限，未來需在更大規(guī)模的奶牛群體中進一步驗證。

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

3 討論

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

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

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

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

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

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

4 結論

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