利用参数函数从荷斯坦牛每天的产奶量推导出泌乳曲线参数的基因组背景

《Journal of Dairy Science》：Genomic background of lactation curve parameters derived from daily milk yield in Holstein cattle using parametric functions

　　帕特里克·R·福尔索-肯莫涅（Patrick R. Fotso-Kenmogne）、保罗·L·S·卡内罗（Paulo L.S. Carneiro）、德尔万·A·席尔瓦（Delvan A. Silva）、西尔琳·F·拉扎罗（Sirlene F. Lázaro）、佩德罗·F·C·阿丰特（Pedro F.C. Aponte）、海梅·A·科布奇（Jaime A. Cobuci）、希纳亚·R·奥利维拉（Hinayah R. Oliveira）和路易斯·F·布里托（Luiz F. Brito）
巴西亚西南部州立大学动物科学研究生项目，伊塔佩廷加（Itapetinga），45700-000，巴伊亚州（BA），巴西

**摘要**
由于泌乳曲线参数（LCP）与泌乳曲线的形状及其生物学解释密切相关，因此在优化奶牛育种计划中至关重要。本研究的主要目标是基于三种非线性函数（Wood [WD]、Wilmink [WL] 和 Ali-Schaeffer [AS]）对美国荷斯坦奶牛的多种 LCP 进行全基因组关联研究和功能富集分析。我们使用了2012年至2019年间出生的11,139头初次产犊牛和6,735头经产牛的2,754,840条和1,642,653条每日产奶量记录。共有14,464头动物进行了基因分型，检测了60,277个单核苷酸多态性（SNP）。通过GBLUP方法估计了SNP效应、它们解释的总加性遗传方差的比例及其近似P值。使用修正的Bonferroni多重检验校正方法识别出显著SNP，以考虑独立染色体片段的数量。随后检查了位于这些显著SNP上游或下游100 kb范围内的基因和数量性状位点，并对每个LCP的候选基因进行了功能富集分析。对于初次产犊牛，WD参数鉴定出81个、WL参数鉴定出128个、AS参数鉴定出196个显著SNP；而对于经产牛，WD参数鉴定出120个、WL参数鉴定出125个、AS参数鉴定出128个显著SNP。这些显著SNP分布在18条常染色体上。一个基因组区域（BTA14: 1,801,116 bp）在所有参数函数（WD、WL和AS）中均存在。位于BTA15和BTA19上的基因组标记仅与WL参数相关，而位于BTA3和BTA20上的标记仅与AS参数相关。在BTA14上发现的显著SNP解释了超过1%的总加性遗传方差。有20个候选基因（ARC、ADGRB1、C8orf90、CYP11B1、DENND3、GML、GPR20、JRK、LY6D、LY6E、LY6L、LYNX1、LYPD2、PSCA、PTK2、PTP4A3、SLURP1、SLC45A4、THEM6和TSNARE1）在所有参数函数中均存在。对于初次产犊牛的WD参数c（泌乳峰值后的斜率下降参数）以及经产牛的AS参数d和f（与斜率上升相关的参数），未发现显著的基因本体术语。先前的研究报告已经鉴定出这些基因组区域内的候选基因，这些基因具有与细胞凋亡、基因表达调控、牛奶生产、临床乳腺炎、牛奶乳糖含量、体细胞评分、脂肪产量和乳房形态相关的生物学功能。本研究有助于识别与LCP相关的多个候选基因，加深了我们对美国荷斯坦奶牛LCP基因组结构的理解。

**解释性总结**
本研究旨在使用三种参数函数（Wood、Wilmink和Ali-Schaeffer）的随机回归模型来识别影响泌乳曲线参数的基因组区域。研究使用了初次和经产牛的每日产奶量记录，并对所有动物进行了基因分型。基于GBLUP方法，估计了基因组标记效应、它们解释的总加性遗传方差的比例及其近似P值。所鉴定的候选基因有助于我们理解LCP背后的基因组结构。

**引言**
由于牛奶生产在该行业具有重要的经济价值，因此它曾一直是奶牛育种计划的主要选择目标（Brito等人，2021；Pedrosa等人，2021）。泌乳曲线的数学模型可以描述牛奶产量的时间变化。一般来说，奶牛的产奶量在分娩后迅速增加达到峰值，随后在一定时期内保持相对稳定（受泌乳持续时间影响），然后逐渐减少直至干奶期。因此，了解泌乳曲线（而不仅仅是总产奶量）非常重要。泌乳曲线的建模有助于表征诸如泌乳峰值和持续时间等重要特征，以及整体抗逆性指标（例如Chen等人，2023；Lazaro等人，2024；Oliveira等人，2024）。已经开发了多种参数函数来建模泌乳曲线，包括Wood（WD，1967）、Wilmink（WL，1987）和Ali & Schaeffer（AS，1987）提出的函数。这些函数具有生物学意义的可解释参数，如a参数（WD和WL）表示分娩后的初始产奶量近似值；b参数（WD和WL）表示直至泌乳峰值产奶量的上升斜率；c参数（WD和WL）表示泌乳峰值后的斜率下降值；a参数（AS）与泌乳峰值产奶量相关；b和c参数（AS）与斜率下降相关；d和f参数（AS）与斜率上升相关（Strucken等人，2012；Lazaro等人，2024）。除了描述泌乳曲线外，揭示其参数的基因组背景也非常重要，因为这些参数具有明确的生物学意义，可以纳入选择指标以改变泌乳曲线的形状。此外，泌乳曲线参数（LCP）是可遗传的（Oliveira等人，2024），这意味着通过选择育种可以改变一个群体的平均泌乳曲线。根据整个泌乳曲线选择动物可以同时考虑泌乳生物学和不同曲线形状的生理后果，有助于减少早期和泌乳高峰期的能量不足，同时提高动物健康、福利和牛奶生产效率（Oliveira等人，2019a）。与此选择策略一致，多项研究报告了牛奶生产的遗传背景。例如，Oliveira等人（2019b）鉴定了加拿大荷斯坦奶牛中与产奶量相关的407个SNP，而Strucken等人（2012）研究了德国奶牛中牛奶生产性状的时变遗传效应，发现多个基因组区域影响牛奶生产。识别与LCP相关的基因组区域和候选基因对于阐明泌乳期间牛奶生产变异的生物学机制至关重要。全基因组关联研究（GWAS）已广泛用于这一目的，但大多数研究集中在牛奶生产性状上，而非LCP。例如，Lazaro等人（2024）使用参数函数对穆拉水牛的泌乳长度、泌乳持续时间、马苏里拉奶酪产量和牛奶生产性状进行了GWAS，并鉴定出多个候选基因。然而，关于奶牛LCP的类似研究仍然有限。因此，本研究的主要目标是：（1）对LCP性状（WD、WL和AS）进行GWAS，以识别影响这些参数的关键基因组区域和候选基因；（2）进行功能基因组分析，以加深对这些性状的生物学途径的理解。全基因组贝叶斯值（GEBV）是使用GLUP方法获得的，该方法遵循了Aguilar等人（2010年）描述的方法。对于每头动物，在整个泌乳期间，GEBV向量是通过以下方式获得的：
GEBVi=Φδ^i，
其中
δ^i 是每头动物i的预测基因组系数向量；而
Φ 是与WD、WL或AS功能相关的每个DIM的独立协变量矩阵。在基于RRM计算出GEBV之后，使用postGSf90软件（Aguilar等人，2014年）来反向求解由SNP效应的加性基因组随机回归系数预测的GEBV。为回归系数计算的SNP效应用于估计随时间变化的SNP效应（从泌乳期的第5天到第305天），详细信息见Oliveira等人（2019b）和Lazaro等人（2024年）的研究。SNP效应可以如下描述（Wang等人，2012年）：
u^c=Z'ZZ'-1a^g，
其中
u^c 是第c个随机回归系数的估计SNP解向量；Z 是包含中心化基因型的矩阵；而
a^g 是通过GBLUP分析估计的第s个随机回归系数的GEBV向量。该向量包含了所有已基因分型动物的第s个随机回归系数。随后，将同一SNP k的所有随机回归系数（对于WD和WL功能，c = 1, 2, 3；对于AS功能，c = 1, 2, …, 5）的SNP解组合成一个向量
u^k=u^k1u^k2u^k3u^k4u^k5'，并使用它来估计所有DIM的SNP效应（从第5天到第305天），如Oliveira等人（2019c）所提出的：
SNPk^=Tu^K，
其中
SNPk^ 是包含第k个SNP每个DIM的SNP效应的向量；T 是每个DIM的协变量矩阵，与WD、WL和AS功能相关；而
u^K 是与第k个SNP相关的所有随机回归系数的SNP解向量。SNPj解释的加性遗传方差百分比计算如下：
%σA,j2=2pj1-pjα^j2σa2×100，
其中
α^ 是SNP效应，
σa2 是总加性遗传方差，pj 是假设与当前种群中等位基因频率相同的参考等位基因频率（VanRaden，2008）。为了识别显著的SNP，SNP效应的近似P值按以下方式计算（Aguilar等人，2019）：
P-valuesij=21-Φa^ijsda^ij，
其中
P-valuesij 是与参数i中的第j个SNP相关的P值，
a^ij 是SNP效应，
而
Φ 是累积标准正态函数（Gualdrón Duarte等人，2014）。基因组膨胀因子（λ）定义为观察到的卡方检验统计量的中位数与卡方分布预期中位数的比率 [λ = median(χ2)/0.456]，用于评估种群分层。SNP效应根据其与每个参数的关联显著性进行排名，使用Bonferroni多重检验校正来定义显著关联的阈值，该阈值为6.08 [?log10 (0.05/60,277)]，其中60,277代表经过QC后的SNP总数。然而，这种更严格的Bonferroni校正被认为是过于保守的。它没有考虑由于SNP之间的强连锁不平衡而导致的测试间的相关性，并且倾向于增加II型错误（假阴性）的风险，可能会遗漏生物学上有意义的关联。因此，为了降低这种假阴性的可能性，还应用了一种修改后的Bonferroni校正。这种方法使用独立染色体片段的数量（Me）而不是SNP的总数。Me是在染色体水平上估计的（Li等人，2015年），基于Goddard等人（2011年）提出的公式：
Me = 2NeL/log(NeL)，
其中Ne是有效种群大小，L是摩尔根氏单位中的染色体长度。采用了Ne=66（Makanjuola等人，2020）。当所有随机回归系数（c = 1, 2, 3对应WD和WL功能；或c = 1, 2, …, 5对应AS功能）的SNP解的近似P值低于修改后的全基因组Bonferroni阈值（0.05/Me）时，这些SNP解被宣布为统计学上显著。在这项研究中，Me估计为2,044.7，得到的?log10显著性阈值为4.61。这种方法在先前的研究中也有应用（Pedrosa等人，2021；Chen等人，2024；Sousa Junior等人，2024）。使用Ensembl BioMart工具根据Bos taurus ARS-UCD1.3参考基因组（Durinck等人，2009）识别可能与选定SNP相关的候选基因。保留位于显著SNP上游或下游100 kb以内的基因和QTL以进行后续分析。QTL信息来自Animal QTLdb（版本56；Hu等人，2022），并使用GALLO R包（Fonseca等人，2020）评估QTL的功能富集。随后，使用DAVID平台（Huang等人，2009；Sherman等人，2022）对每个LCP识别的候选基因进行功能富集分析。当假发现率（FDR；Benjamini和Hochberg，1995）等于0.05时，基因和基因本体（GO）术语被认为是显著富集的。

表1显示了第一胎和第二胎母牛的WD模型中的参数a、b和c；WL模型中的参数a、b和c；以及AS模型中的参数a、b、c、d和f的平均值（±SE）、最小值和最大值。第一胎和第二胎母牛的遗传率估计如图1所示。对于第一胎母牛，当拟合WD功能时，遗传率估计范围为0.14至0.21，WL为0.14至0.20，AS为0.10至0.21；对于第二胎母牛，WD为0.10至0.16，WL为0.11至0.16，AS为0.11至0.16。

描述性统计：每个参数和遗传率
表1显示了美国荷斯坦牛每日产奶量（DMY）的Wood、Wilmink和Ali-Schaeffer模型中参数a、b和c的平均值（±SE）、最小值和最大值。
模型 PMean ± SEM 最小值 最大值
第一胎母牛 Wood a 25.3 10 ± 0.01 18.5 0 34.6
b 2.7 0.001 2.0 3.5
c 4.789 × 10?4 ± 2.37 × 10?7 0 0.3
d 4.789 × 10?4 ± 2.37 × 10?7 0 0.6
e 4.789 × 10?4 ± 2.37 × 10?7 0 0.3
Wilmink a 12.3 9.7 16.3 0 30.3
b 23.1 0.010 17.4 0 30.3
c 0.239 × 10?3 ± 1.16 × 10?7 0 0.3
d 0.239 × 10?3 ± 1.16 × 10?7 0 0.3
e 0.239 × 10?3 ± 1.16 × 10?7 0 0.3
Ali-Schaeffer a 26.7 0.011 20.0 36.1 0 27.3
b 17.7 0.013 11.1 0 27.3
c 33.6 0.026 20.2 51.3 0 15.9
d 11.2 0.006 7.6 15.9 0 15.9
f 0.286 0.000 20.1 0 0.4

图1. 使用Wood（WD）、Wilmink（WL）和Ali-Schaeffer（AS）函数在美国荷斯坦牛中估计的遗传率。(a)代表第一胎母牛，(b)代表第二胎母牛。

关联分析：λ值的范围为0.91至0.99，所有LCP的分位数-分位数（Q-Q）图显示在补充文件S1中（补充文件S1：图S1 - S6；https://doi.org/10.6084/m9.figshare.31045159）。根据Q-Q图和λ值，没有未校正的种群分层的迹象。补充文件S1还提供了所有LCP的曼哈顿图（补充文件S1：图S7 - S12；https://doi.org/10.6084/m9.figshare.31045159）。在显著SNP上游和下游100 kb范围内发现的QTL分别提供在补充文件S2（第一胎母牛）和S3（第二胎母牛）中。补充文件S4和S5分别提供了第一胎和第二胎母牛的解释加性遗传方差百分比的显著SNP。

表2总结了主要的显著SNP和候选基因，完整结果分别在表3和表4中展示。对于第一胎母牛，分别识别出81个、128个和196个WD、WL和AS参数的显著SNP（表3）。对于第二胎母牛，分别识别出120个、125个和128个WD、WL和AS参数的显著SNP（补充表S5）。此外，BTA14上有最多的显著SNP，第一胎母牛有358个，第二胎母牛有350个（表3，表4）。此外，当基于SNP总数考虑更严格的Bonferroni校正阈值时，显著SNP大多位于BTA14上，且一个基因组区域（BTA14: 1,801,116 bp；PSCA, JRK, ARC, ADGRB1）是所有参数函数共有的。总体而言，与基因组区域重叠的大部分QTL之前与产奶性状相关，但也有一些与其他性状（生长、繁殖和健康）相关（补充文件S2和S3；https://doi.org/10.6084/m9.figshare.31045159）。对于第一胎母牛，位于BTA14上的49个、89个和139个显著SNP分别解释了WD、WL和AS参数总加性遗传方差的1%以上；而对于第二胎母牛，分别识别出93个、101个和115个显著SNP（补充文件S4和S5；https://doi.org/10.6084/m9.figshare.31045159）。功能注释分析中最显著术语的摘要显示在表5和表6中，完整结果在补充文件S6中（https://doi.org/10.6084/m9.figshare.31045159）；发现了共同的生物学功能。

表2. 与美国荷斯坦牛每个泌乳曲线参数相关的最显著SNP的总结
模型 PChr 位置（bp） P值 基因符号 方差（%）
第一胎母牛 Wood BTA14 1,801,116 1.84 × 10?8 PSCA, JRK, ARC, ADGRB1
b BTA14 1,801,116 3.98 × 10?2 PSCA, JRK, ARC, ADGRB1
c BTA14 35,7 2,773 1.37 × 10?6 TRPA1
Wilmink BTA14 1,801,116 5.28 × 10?2 PSCA, JRK, ARC, ADGRB1
d BTA14 1,801,116 1.31 × 10?2 PSCA, JRK, ARC, ADGRB1
e BTA14 1,801,116 4.86 PSCA, JRK, ARC, ADGRB1
f BTA14 1,736,599 7.41 × 10?11 LY6D, LYNX1, LYPD2, SLURP1, THEM
PSCA, JRK, ARC, ADGRB1
Ali-Schaeffer BTA14 1,801,116 3.95 × 10?3 PSCA, JRK, ARC, ADGRB1
b BTA14 1,801,116 2.46 × 10?2 PSCA, JRK, ARC, ADGRB1
c BTA14 1,801,116 2.82 × 10?2 PSCA, JRK, ARC, ADGRB1
d BTA14 1,801,116 1.55 PSCA, JRK, ARC, ADGRB1
e BTA14 1,801,116 3.50 × 10?13 PSCA, JRK, ARC, ADGRB1
f BTA14 1,801,116 4.23 × 10?7 PSCA, JRK, ARC, ADGRB1

表3. 在美国荷斯坦牛中，位于泌乳曲线参数显著窗口±100 kb内的候选基因
模型 PChr 位置（bp） P值 基因符号 方差（%）
第一胎母牛 Wood BTA14 1,801,116 1.84 × 10?8 PSCA, JRK, ARC, ADGRB1
b BTA14 1,801,116 3.98 × 10?2 PSCA, JRK, ARC, ADGRB1
c BTA14 35,7 2,773 1.37 × 10?6 TRPA1
Wilmink BTA14 1,801,116 5.28 × 10?2 PSCA, JRK, ARC, ADGRB1
d BTA14 1,801,116 1.31 × 10?2 PSCA, JRK, ARC, ADGRB1
e BTA14 1,801,116 4.86 PSCA, JRK, ARC, ADGRB1
f BTA14 1,736,599 7.41 × 10?11 LY6D, LYNX1, LYPD2, SLURP1, THEM
PSCA, JRK, ARC, ADGRB1
Ali-Schaeffer BTA14 1,801,116 3.95 × 10?3 PSCA, JRK, ARC, ADGRB1
b BTA14 1,801,116 2.46 × 10?2 PSCA, JRK, ARC, ADGRB1
c BTA14 1,801,116 2.82 × 10?2 PSCA, JRK, ARC, ADGRB1
d BTA14 1,801,116 5.85 PSCA, JRK, ARC, ADGRB1
e BTA14 1,801,116 3.50 × 10?13 PSCA, JRK, ARC, ADGRB1
f BTA14 1,801,116 4.23 × 10?7 PSCA, JRK, ARC, ADGRB1SNPs
起始位置
结束位置
基因符号

Wooda
BTA14
201,636,599
3,040,147
LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, GPR20, SLC45A4, DENND3, PTK2, MIR151
A
BTA17
0171,181,463
71,381,463
ZNF70, C17H22
orf15, CHCHD10, MMP11, SMARCB1, DERL3, SLC2A11, MIF, GSTT1, CABIN1, bta-mir-2893
BTA29
0148,353,418
48,553,418
MRGPRG, OSBPL5, U6
b
BTA50
195,076,510
95,276,510
ART4, C5H12
orf60, WBP11, H2AJ, H4C16, GUCY2C, PLBD1

BTA144
81,389,496
35,832,773
LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, GPR20, SLC45A4, DENND3, PTK2, MIR151

表4. 在美国荷斯坦奶牛（第二次产犊的奶牛）中，与泌乳曲线参数显著相关的候选基因，这些基因位于±100 kb的窗口内

模型
Chr
SNPs数量
起始位置
结束位置
基因符号

Wooda
BTA8
017,877,499
8,077,499
FAM16
7A, TDH, MTMR9
BTA9
0289,312,796
89,947,787
FBXO5, MTRF1L, RGS17
b
BTA14
251,488,879
2,189,613
GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1
BTA10
151,733,017
151,933,017
COL6A5
BTA13
0155,500,716
55,700,716
CDH4
BTA146
21,208,359
3,040,147
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A
BTA6
0188,819,352
89,019,352
CXCL5, CXCL2
BTA142
61,208,359
2,247,133
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TRPA1

BTA18
013,646,145
3,846,145
CNTNAP4
c

BTA20
112,092,595
126,292,595
SLC9A1, TENT5B, TRNP1, KDF1, NUDC
BTA70
3105,312,109
105,565,008
—
BTA90
194,246,873
94,446,873
ZDHHC14, SNX9
BTA140
235,632,773
68,649,959
TRPA1

BTA18
0237,228,355
61,538,125
—
Wilminka

BTA6
183,751,252
89,019,352
LCORL, TMPRSS11
D, TMPRSS11
A, ANKRD17, ALB, AFP, AFM, RASSF6, CXCL5, CXCL2
BTA143
41,208,359
2,858,369
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TRPA1

BTA180
13,646,145
3,846,145
CNTNAP4
c

BTA20
126,092,595
127,914,244
SLC9A1, TENT5B, TRNP1, KDF1, NUDC, SYF2
BTA70
3105,312,109
105,565,008
—
BTA90
194,246,873
94,446,873
ZDHHC14, SNX9

BTA140
235,632,773
68,649,959
TRPA1

BTA180
237,228,355
61,538,125
—
Wilminka

BTA6
183,751,252
89,019,352
LCORL, TMPRSS11
D, TMPRSS11
A, ANKRD17, ALB, AFP, AFM, RASSF6, CXCL5, CXCL2

BTA143
41,208,359
2,858,369
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TRPA1

BTA180
13,646,145
3,846,145
CNTNAP4
c

BTA20
112,092,595
126,292,595
SLC9A1, TENT5B, TRNP1, KDF1, NUDC, SYF2
BTA70
3105,312,109
105,565,008
—
BTA90
194,246,873
94,446,873
ZDHHC14, SNX9

BTA140
235,632,773
68,649,959
TRPA1

BTA140
237,228,355
61,538,125
—
Wilminka

BTA50
195,076,510
95,276,510
ART4, C5H12
orf60, WBP11, H2AJ, H4C16, GUCY2C, PLBD1

BTA144
81,389,496
66,871,994
LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TSPYL5

BTA180
13,646,145
3,846,145
CNTNAP4
c

BTA20
126,092,595
127,914,244
SLC9A1, TENT5B, TRNP1, KDF1, NUDC, SYF2
BTA70
3105,312,109
105,565,008
—
BTA142
11,208,359
69,609,386
TP53
INP1, CCNE2, INTS8, DPY19L4, ESRP1, GML, CYP11B1, LY6D, LYNX1, SLURP1, THEM6, PSCA, ZC3H3, JRK, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, ARC, ADGRB1, LY6L, TSNARE1, MAFA, LYPD2

Ali-Schaeffer

BTA50
195,076,510
95,276,510
ART4, C5H12
orf60, WBP11, H2AJ, H4C16, GUCY2C, PLBD1

BTA144
71,389,496
66,871,994
LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TSPYL5

BTA180
13,646,145
3,846,145
CNTNAP4
c

BTA20
212,092,595
127,914,244
SLC9A1, TENT5B, TRNP1, KDF1, NUDC, SYF2
BTA70
3105,312,109
105,565,008
—
BTA142
11,208,359
69,609,386
TP53
INP1, CCNE2, INTS8, DPY19L4, ESRP1, GML, CYP11B1, LY6D, LYNX1, SLURP1, THEM6, PSCA, ZC3H3, JRK, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, ARC, ADGRB1, LY6L, TSNARE1, MAFA, LYPD2

表4. 在美国荷斯坦奶牛（第二次产犊的奶牛）中，与泌乳曲线参数显著相关的候选基因，这些基因位于±100 kb的窗口内

模型
Chr
SNPs数量
起始位置
结束位置
基因符号

Wooda
BTA80
17,877,499
8,077,499
FAM16
7A, TDH, MTMR9
BTA9
0289,312,796
89,947,787
FBXO5, MTRF1L, RGS17
b
BTA14
251,488,879
2,189,613
GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1
BTA10
151,733,017
151,933,017
COL6A5
BTA13
0155,500,716
55,700,716
CDH4
BTA146
21,208,359
3,040,147
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TRPA1

BTA180
237,228,355
61,538,125
—
Wilminka

BTA6
183,751,252
89,019,352
LCORL, TMPRSS11
D, TMPRSS11
A, ANKRD17, ALB, AFP, AFM, RASSF6, CXCL5, CXCL2

BTA143
41,208,359
2,858,369
ZC3H3, MAFA, RHPN1, ZNF696, GLI4, GPIHBP1, LY6H, LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, C8
orf90, PTP4A3, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TRPA1

BTA144
81,389,496
35,832,773
LY6L, LY6E, GML, CYP11B1, LY6D, LYNX1, LYPD2, SLURP1, THEM6, PSCA, JRK, ARC, ADGRB1, TSNARE1, GPR20, SLC45A4, DENND3, PTK2, MIR151A, TSPYL5

BTA180
13,646,145
3,846,145
CNTNAP4
c

BTA20
126,092,595
127,914,244
SLC9A1, TENT5B, TRNP1, KDF1, NUDC, SYF2
BTA70
3105,312,109
105,565,008
—
BTA142
11尽管基因组的遗传率较低到中等，但它们更为准确。第一胎和第二胎母牛的遗传率估计值在分娩次数（DIM）和胎次之间存在差异（图1）。这些差异可能反映了哺乳期间的生理变化，而不同胎次之间的差异可能是由于选择决策改变了MY的遗传变异。使用相同的数据集和四阶勒让德多项式，Fotso-Kenmogne等人（2025年）报告的日奶产量的遗传率估计值在0.03到0.22之间。这种轻微的差异可以归因于用于估计（协）方差分量的系数。许多位于BTA14上的重要SNP（补充文件S4和S5）解释了总加性遗传变异的相当大比例，表明BTA14是LCP遗传结构的重要组成部分。一些候选基因在Wood、Wilmink和Ali-Schaeffer功能中都有重叠（图2）。这可以解释为这些功能（WD、WL和AS）描述的是相同的性状。在这里，我们主要讨论位于最显著区域的候选基因（表2）。

**Wood函数的参数**
对于第一胎和第二胎母牛，共发现了148个候选基因（表3、表4）。a参数和b参数的候选基因数量多于c参数。总体而言，所有参数（a、b和c）都鉴定出了12个相似的候选基因（ADGRB1、ARC、CYP11B1、GML、JRK、LY6D、LYNX1、LYYPD2、PSCA、SLURP1、THEM6和TSNARE1）。有8个候选基因（SLC9A1、TENT5B、TRNP1、KDF1、NDUC、ZDHHC14、SNX9和TRPA1）仅适用于第一胎母牛的c参数，这表明它们可能在该参数所基于的生物过程中起作用。在最显著的区域中，BTA14包含PSCA（前列腺干细胞抗原）、JRK（Jrk螺旋-转角-螺旋蛋白）、ARC（活性调控的细胞骨架相关蛋白）、ADGRB1（粘附G蛋白偶联受体B1）和TRPA1（瞬时受体电位阳离子通道亚家族A成员）基因。TRPA1基因是仅属于Wood函数的6个候选基因之一（图2）。PSCA基因属于LY-6表面蛋白家族，参与主要组织相容性复合物的调节，并且据报道与首次产犊的美国荷斯坦奶牛的临床乳腺炎有关（Tiezzi等人，2015年）。此外，ARC基因也与临床乳腺炎相关特征和体细胞评分有关（Narayana等人，2023年）。鉴于ARC与乳腺炎相关，它可能与哺乳高峰产量后的下降斜率参数有关。c参数代表哺乳高峰产量后的下降斜率，可以被认为是哺乳持续性的一个衡量指标。JRK基因与多胎母牛的305天日奶产量和高峰产量有关（Atashi等人，2020年）。此外，ADGRB1基因与牛奶中的乳糖含量有关（Costa等人，2019年），而TRPA1与生长调节有关（Aponte等人，2024年）。生长调节可能促进乳腺细胞增殖，从而增加日奶产量和哺乳持续性。与此一致的是，Wu等人（2019年）证明TRPA1在大家畜品种（即嘉兴红牛、秦川牛和泸西牛）中调节身体和器官组织的生长与发育。在BTA2上，KDF1（角质形成细胞分化因子1）和NUDC（核分布C，动力蛋白复合体调节因子）被确定为表皮发育和乳房形态的关键调节因子（Miles等人，2021年）。先前的研究表明，各种QTL（补充文件S2和S3）位于该显著SNP周围的基因组区域，包括milking speed。目前尚不清楚是否存在基因型与挤奶方式之间的交互作用。高产奶牛的生产性能、寿命和福利高度依赖于它们所处的环境与其遗传背景之间的交互作用。因此，有必要探索潜在的基因型与挤奶方式之间的交互作用。

**Wilmink函数的参数**
对于第一胎和第二胎母牛，分别发现了114个和84个候选基因（表3、表4）。位于BTA14上的13个候选基因（ADGRB1、ARC、CYP11B1、GML、JRK、LY6D、LYNX1、LYYPD2、LY6L、PSCA、SLURP1、THEM6和TSNARE1）在第一胎母牛的所有参数（a、b和c）中都存在。LY6D（淋巴细胞抗原6家族成员D）基因与脂肪产量和奶产量有关（Suchocki等人，2016年）以及乳腺炎（Tiezzi等人，2015年）。根据Tiezzi等人（2015年）的研究，LY6D基因编码的是淋巴细胞抗原6复合物的组成部分，而中性粒细胞调节功能的控制也对抵抗乳房感染起着重要作用。LYNX1（Ly6/神经毒素1）基因被报告为与奶产量相关的候选基因（Hosseinzadeh等人，2024年）。LYYPD2（LY6/PLAUR结构域2）之前已被报道与瓦隆荷斯坦奶牛的奶酪生产特征和日奶产量相关（Atashi等人，2023年）。SLURP1（分泌型LY6/PLAUR结构域1）是淋巴细胞抗原-6复合体（LY6）的一部分，以其中性粒细胞调节功能而闻名，并与临床乳腺炎有关（Tiezzi等人，2015年）。THEM6（硫酯酶超家族成员6）基因与乳糖产量有关（Costa等人，2019年），并编码一种硫酯酶。TSNARE1（T-SNARE结构域1）基因在突触囊泡外泌和细胞内蛋白质转运中起作用，并且之前与荷斯坦奶牛在热应激期间的直肠温度有关（Luo等人，2021年）。这种效应可能与哺乳曲线高峰产量后的斜率有关，因为热应激下的奶牛往往产奶量减少。CXCL5（C-X-C基序趋化因子配体5）和CXCL2（C-X-C基序趋化因子配体2）基因与荷斯坦奶牛的乳腺炎有关（Islam等人，2020年；Lázaro等人，2024年）。Atashi等人（2020年）报告称，THEM6、SLURP1、LYNX1、LYYPD2和LY6D与荷斯坦多胎母牛的305天日奶产量和高峰产量有关。第一胎和第二胎母牛中最显著的术语与LY6 UPA受体类似。小鼠蛋白质组中大约有60种蛋白质被归类为Ly6/uPAR家族。这些蛋白质是膜锚定的，参与多种功能，包括受体结合和细胞粘附。Ly6/uPAR家族成员在细胞信号传导中起调节作用，尤其是在大脑和先天免疫系统的细胞中（Tirosh等人，2013年）。这些功能表明它们可能在细胞失衡中起作用，从而导致日奶产量减少（哺乳高峰产量后的下降斜率参数）。

**Ali和Schaeffer函数的参数**
共鉴定出286个候选基因，其中149个来自第一胎母牛，137个来自第二胎母牛（表3、表4）。十个候选基因（ARC、ADGRB1、JRK、LY6D、LYNX1、LYYPD2、PSCA、SLURP1、THEM6和TSNARE1）在第一胎母牛的所有参数（a、b、c、d和f）中都存在。在d和f参数的最显著区域中，BTA11包含C1D（核受体辅阻遏物）、DNAAF10（动力蛋白轴丝组装因子10）、PNO1（NOB1同源物的伴侣）、PPP3R1（蛋白质磷酸酶3调节亚基B，α）和CNRIP1（大麻素受体相互作用蛋白1）基因。C1D基因与凋亡过程和基因表达调控有关（Rothbarth等人，2001年）。凋亡是一种程序性细胞死亡形式，当细胞接收到内部或外部信号时启动，并通过一系列生化事件进行。DNAAF10基因已被证明编码参与轴丝动力蛋白组装的蛋白质（Braschi等人，2022年）。DNAAF10也可能与凋亡过程有关，这可能与我们观察到的高峰产量后的下降斜率有关，因为更多的凋亡可能导致日奶产量进一步减少。PNO1基因调节肝细胞的自噬和凋亡（Han等人，2021年）。PNO1是核糖体生物合成的关键成分，也被认为是与脂肪沉积相关的候选基因（Moscarelli等人，2020年）。PPP3R1基因编码一种磷酸酶蛋白，并与Nellore奶牛的肌肉生长、发育和降解有关（Santos Silva等人，2020年）。这些过程可能影响日奶产量，特别是在第一胎母牛的早期和高峰哺乳期，这一阶段的特征是负能量平衡。CNRIP1基因似乎参与调节CNR1（大麻素受体1）的部分功能，从而调节受体活性。CNRIP1也在绒毛膜上被鉴定出来，可能增强脂质代谢并调节氧化应激（dos Santos Silva等人，2024年）。因此，由于其在氧化应激相关过程中的作用，该基因可能是与b和c参数相关的潜在候选基因。

本研究鉴定了许多与美洲荷斯坦奶牛不同日奶产量（LCP）相关的SNP、QTL、基因组区域和候选基因。为了完整性，应用了基于SNP总数的严格Bonferroni多重检验校正以及考虑独立染色体段数量的修改后的Bonferroni校正。类似的校正策略已在之前的基因组研究中使用（Pedrosa等人，2021年；Chen等人，2024年；Sousa Junior等人，2024年）。传统的Bonferroni方法假设测试之间的独立性，但在存在强 linkage disequilibrium（LD）的情况下这一假设不成立。此外，这种方法通常被认为过于保守，可能会显著增加未能检测到生物学相关基因组区域的可能性，从而导致假阴性结果。因此，我们主要讨论了最显著的SNP，同时在补充文件S4和S5中报告了其余的关联。这项研究揭示了与LCP相关的候选基因，有助于我们理解LCP的遗传结构及其对乳品生产的影响。下一步应使用全基因组序列数据进行额外的GWAS分析，因为这种方法可能有助于发现更多的遗传关联；并且还应对其他牛品种进行类似的分析。此外，基于基因敲除或基因编辑验证这些关联将是另一个重要的步骤。

**结论**
从美洲荷斯坦奶牛的日奶产量中得出的哺乳曲线参数的基因组背景揭示了许多与Wood、Wilmink和Ali-Schaeffer参数显著相关的SNP，这些SNP分布在18条染色体上。有20个候选基因（LY6D、LYNX1、LYYPD2、SLURP1、THEM6、PSCA、JRK、ARC、ADGRB1、TSNARE1、GPR20、SLC45A4、DENND3、PTK2、LY6L、LY6E、GML、CYP11B1、C8orf90和PTP4A3）在所有参数函数（Wood、Wilmink和Ali-Schaeffer）中都存在。这些候选基因参与了多种生物途径，如调节功能、信号受体活性和凋亡过程。鉴定的基因组区域与之前研究中报道的几个QTL有重叠。这些发现有助于更好地理解哺乳曲线参数的基因组基础，这将有助于开发针对荷斯坦奶牛的定向育种策略，从而提高生产力和乳品行业的整体可持续性。