整合多组学数据以识别影响后期产蛋母鸡产蛋性能的关键基因

《Poultry Science》：Integrated multi-omics data to identify the key genes affecting egg-laying performance in late-phase hens

　　杨洪瑞|王磊|李荣|南九红|蔡 Jinping|赵云霞|朱桂瑜|李世军
中国华中农业大学动物科学技术学院与兽医学院，教育部农业动物遗传、育种与繁殖重点实验室，武汉，430070

**摘要**
在整个产蛋周期中，后期产蛋阶段是最长的时期，它直接决定了总产蛋量并影响淘汰时间的选择。为了研究这一阶段产蛋性能下降的机制，本研究选择了文昌鸡作为研究对象，这种鸡具有较低的产蛋率和较高的表型变异。我们研究了872只文昌鸡在后期产蛋阶段的产蛋模式，并对3只高产蛋率（产蛋率>76.0%）的母鸡和3只低产蛋率（产蛋率<50.0%）的母鸡进行了机制研究。通过对这些母鸡的卵巢间质组织进行RNA-seq和ATAC-seq分析，我们发现了800个差异表达基因（p<0.05且|log?FC|≥1），其中包括317个上调基因和483个下调基因，这些基因富集在免疫反应、内肽酶活性负调控、紧密连接以及神经活性配体-受体相互作用等生物过程中。ATAC-seq分析在两组中分别发现了1065个和428个潜在的差异可及区域，其中9.8%和5.8%的区域位于启动子区域。综合分析确定了三个可能影响后期产蛋性能的关键基因，包括MICA、MHCIA6和MHCIY。这些基因都属于非经典MHC-Ⅰ家族，与免疫调节有关。此外，还预测有四种转录因子（HOXB6、GATA3、TFAP2A和WT1）可能调控这些关键基因。总体而言，我们的研究结果表明，后期产蛋阶段卵巢免疫稳态的失衡可能是产蛋量下降的重要因素。这项研究不仅提供了维持高产蛋性能的机制见解，还为文昌鸡的遗传改良提供了潜在的分子靶点，并推进了对家禽生殖调控的理解。

**引言**
母鸡的产蛋周期包括三个不同的阶段：前期、高峰期和后期产蛋阶段。值得注意的是，后期产蛋阶段占总周期时间的50%以上。这一阶段的产蛋率限制了总产蛋量，是提高整体产蛋量的关键因素。然而，在后期产蛋阶段，随着年龄的增长，产蛋率会逐渐下降，尤其是在肉鸡品种中（Zhang等人，2021年）。先前的研究表明，这一阶段的产蛋率具有很强的遗传特性和群体间的差异。在一條肉鸡杂交系中，52周龄时产蛋率的遗传力和变异系数分别为0.426和33.3%（Chomchuen等人，2022年）。全基因组关联研究（GWAS）也鉴定了多个影响鸡产蛋性状的SNP和候选基因，如SHROOM2、SYNE3和CNNM2（Tan等人，2025年）。因此，针对这一阶段的遗传研究对于提高文昌鸡的产蛋性能至关重要。

卵巢是鸟类生殖系统的关键组成部分，其生理状态在调节生殖性状中起着重要作用。鸡的卵巢由卵泡和卵巢间质组成。后者对维持正常的卵泡发育至关重要，分为皮质和髓质两部分。皮质位于外层，是卵泡生长和成熟的地方；髓质位于内部，包含血管、神经、淋巴管、免疫细胞、细胞外基质及其他卵巢特异性成分。这些元素共同构成了一个复杂的微环境，在卵泡发育、排卵调节和维持卵巢免疫稳态中发挥关键作用（Kinnear等人，2020年）。在后期产蛋阶段，卵巢老化加剧了氧化应激和局部炎症，导致卵泡闭锁、颗粒细胞凋亡，最终影响产蛋量（Hao等人，2020年；Zeng等人，2024年）。多项研究表明，LECT2、BMP4和COL1A1等基因可以通过调节卵巢间质内的炎症反应、缓解内质网应激和促进细胞周期进展来促进颗粒细胞增殖，从而保证正常的卵泡发育（Yao等人，2020年；Chen等人，2024b；Li等人，2025年）。尽管有证据表明卵巢功能下降是导致产蛋量减少的关键因素，但这一过程中产蛋量下降的分子机制仍有待进一步阐明。

目前，RNA-seq在禽类产蛋性状研究中得到了广泛应用，已用于鉴定许多关键基因和途径，包括GRIA1、PRLR、THBS2和IGF1以及神经活性配体-受体相互作用、ECM-受体相互作用、免疫功能、MAPK信号传导和钙信号传导等途径（Zhang等人，2019年；Mu等人，2021年；Bhavana等人，2022年；He等人，2022年）。然而，单一组学方法无法完全阐明后期产蛋阶段卵巢功能下降的调控机制。基因表达的差异通常源于上游调控元件活性的改变，而染色质的可及性是其功能能力的关键决定因素。ATAC-seq能够绘制全基因组染色质的可及性图谱，识别出活性调控元件，如启动子和增强子（Buenrostro等人，2013年）。将ATAC-seq与RNA-seq结合使用，可以通过将染色质可及性与基因转录水平相关联来阐明基因调控机制。这种多组学方法越来越多地用于揭示复杂的性状（如生殖和器官发育）的遗传调控机制。Zhang等人（2023年）发现了与性反转鸡生殖腺发育相关的FGFR3基因。Li等人（2022年）表明，鸡颗粒细胞的转录活性受染色质结构动态变化的调控。Dinh等人（2023年）证明排卵刺激显著影响小鼠的细胞染色质可及性。因此，这种多组学整合方法为阐明后期产蛋母鸡卵巢内的调控网络提供了有效途径。

文昌鸡是中国一种著名的黄羽品种，以其优良的肉质而受到重视。然而，该品种的产蛋性能较差，这对育种企业的经济效益产生了不利影响。因此，研究如何提高文昌鸡的产蛋量对家禽产业具有重要意义。此外，其相对较短的驯化历史保留了丰富的遗传多样性，表现为产蛋量的显著表型变异，表明具有较大的选择育种潜力（Xing等人，2020年；Ren等人，2024年）。因此，文昌鸡是研究卵巢差异如何调节产蛋率的理想模型。

**材料与方法**
**伦理声明**
本研究中的所有动物实验均获得了华中农业大学动物福利委员会的批准（批准编号：HZAUCH-2024-0044）。

**动物及卵巢组织的采集**
共收集了872只文昌鸡，来自中国海南（Tanniu）文昌鸡有限公司，用于后续实验。这些鸡在12周龄后被转移到同一鸡舍的单个笼子里。它们按照中国国家标准“黄羽肉鸡的营养需求”（NY/T 3645 2020）提供玉米-大豆粉饲料，并自由饮水。产蛋开始后，每周光照时间增加1小时，直至达到16小时，并保持不变。从44周至49周龄期间，每天手动记录每只鸡的产蛋量。根据这些数据，随机选择了3只产蛋率>76.0%的母鸡组成高产蛋组（HEL），以及3只产蛋率<50.0%的母鸡组成低产蛋组（LEL）。HEL组包括H-1、H-2和H-3，LEL组包括L-1、L-2和L-3。对这些鸡实施了安乐死并采集了血液样本。解剖腹部后，取出了卵巢。去除卵泡中的卵黄，用冷PBS清洗剩余组织。取两部分卵巢间质组织用于后续文库构建和数据分析。一部分组织切碎后加入Tritol试剂（15-596-018，Invitrogen，美国卡尔斯巴德），冷冻在液氮中并储存在-80°C以提取总RNA；另一部分组织也冷冻在液氮中并储存在-80°C用于ATAC-seq文库构建。剩余组织放入4%甲醛（PFA）固定液中（BL539A，Biosharp，中国安徽）进行形态学检查。

**染色和石蜡包埋**
HEL组和LEL组的卵巢组织样本在4% PFA固定液中固定24小时，然后用流动水冲洗过夜。随后进行一系列步骤，包括分级乙醇脱水、透明处理、石蜡渗透和包埋。使用切片机从石蜡块中切割出4微米厚的切片，并将其安装在玻璃载玻片上。切片用苏木精和伊红（H&E）染色，然后用中性香脂封片。制备好的卵巢组织切片在立式显微镜下观察。

**RNA-seq文库构建**
使用酚氯仿法从卵巢组织中提取总RNA。利用NanoDrop 2000分光光度计（Thermo Fisher Scientific，美国马萨诸塞州沃尔瑟姆）评估RNA浓度和质量。然后用Agilent 2100 Bioanalyzer仪器（Agilent Technologies，美国加利福尼亚州帕洛阿尔托）检测RNA完整性。RIN值大于8.0的样本用于mRNA文库构建（根据Illumina公司的TruSeq mRNA Sample Preparation Kit），并在Illumina NovaSeq 6000测序系统（Illumina，美国加利福尼亚州圣地亚哥）上进行双端测序。样本H-1、H-2和H-3来自HEL组，L-1、L-2和L-3来自LEL组。原始序列数据已存入中国国家生物信息中心（CNCB）的基因组序列档案（GSA）（Zhang等人，2025年），BioProject ID为PRJC054005。

**RNA-seq数据分析**
使用FastQC（v0.12.1）（Wingett和Andrews，2018年）对原始测序数据进行质量控制。保留Q30>80%且GC含量在40%至60%之间的文库用于后续分析。清洗后的读段与鸡参考基因组（galGal6，GCF_000002315.5）对齐，然后使用RSEM（v1.3.1）（Li和Dewey，2011年）和 splice-aware STAR对齐器（v2.7.10a）（Dobin等人，2012年）进行转录本水平定量。定量参数设置为：–paired-end –star –calc-ci –strandedness reverse，基因表达水平以每百万碱基对的片段数（FPKM）表示。使用SAMtools（v1.17）（Danecek等人，2021年）对结果对齐文件进行排序并转换为BAM格式，以便后续可视化。使用R中的prcomp函数（v4.5.1）对样本进行主成分分析（PCA）。

**差异表达基因的鉴定和功能富集分析**
使用edgeR包（v4.0.16）（Robinson等人，2010年）进行差异表达分析。基于倍数变化（FC）评估基因表达。差异表达基因（DEGs）的阈值设定为p<0.05且|log?(FC)|≥1。使用KOBAS（Bu等人，2021）对DEGs进行GO和KEGG通路富集分析，以鸡（Gallus gallus）为参考物种。从STRING数据库（Szklarczyk等人，2022年）中筛选出相互作用得分≥0.40（中等置信度）的DEGs的蛋白质-蛋白质相互作用（PPI）网络。将PPI网络导入Cytoscape软件（v3.9.1）（Shannon等人，2003年）进行可视化分析。

**RT-qPCR**
为了验证差异表达基因，从鉴定出的DEGs中随机选择了9个基因，以GAPDH作为内源性参考基因进行RT-qPCR。引物的特异性通过NCBI Primer-BLAST分析和常规PCR后的琼脂糖凝胶电泳进行了验证（图S2，表S10）。总RNA使用HiScript II Q RT SuperMix for qPCR（+gDNA wiper）试剂盒（R223-01，Vazyme，江苏，中国）逆转为cDNA。然后使用ChamQ通用SYBR qPCR混合液（Q711-02，Vazyme，江苏，中国）在ABI实时PCR系统（Applied Biosystems，加利福尼亚州卡尔斯巴德，美国）上进行RT-qPCR。统计分析基于ΔCt值（ΔCt = Ct目标基因 - Ct GAPDH）。两组之间的统计显著性通过非配对双尾 Student's t检验确定。由于两组之间只进行了一次比较，因此未应用多重比较校正。

ATAC-seq文库的构建
冷冻组织在液氮中粉碎，并处理大约0.05克的粉末。除非另有说明，所有后续步骤均在4°C下进行。样品重新悬浮在1毫升冰冷的1×DPBS溶液中，然后在2500 rcf下离心5分钟以去除杂质。与1毫升冰冷的裂解缓冲液孵育3分钟后释放细胞核。使用血细胞计数器和台盼蓝染色在显微镜下计数细胞核。每个样品共50,000个细胞核，用Tn5转座酶（A0215，YINGZI GENE，湖北，中国）在37°C下处理1小时。产物使用DNA清洁试剂盒（TD413，GENSTONE BIOTECH，北京，中国）进行纯化，并通过Qubit荧光法（Thermo Fisher Scientific，马萨诸塞州沃尔瑟姆，美国）进行定量。PCR扩增后，产物分别使用0.5×、0.8×和1.1×DNA清洁珠子纯化并回收，以构建最终的ATAC-seq文库。配对末端测序在Illumina NovaSeq 6000上进行。样品H-1和H-2来自HEL组，L-1和L-2来自LEL组。所有样品也进行了RNA-seq。原始序列数据已存放在中国国家生物信息中心（CNCB）的基因组序列存档（GSA）（Zhang等人，2025）网站上，BioProject ID：PRJC054005。

ATAC-seq的数据分析
初步过滤后，使用FastQC（v0.12.1）（Wingett和Andrews，2018）评估数据质量。使用Trimmomatic（v0.39）（Bolger等人，2014）去除接头和低质量碱基，参数设置为：ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36。修剪后，进行第二轮FastQC以确认保留的读取序列的Q30 > 80%，GC含量为40%至60%，且无残留的接头污染。清洁后的配对末端ATAC-seq读取序列使用BWA-MEM（Li，2013）与鸡的参考基因组（galGal6，GCF_000002315.5）对齐。SAM文件使用SAMtools转换为坐标排序的、索引的BAM文件。然后通过MAPQ ≥ 30的适当配对对齐的读取序列进行过滤，并在下游分析前去除线粒体（chrM）读取序列。这些BAM文件随后使用deepTools（v3.5.1）（Ramírez等人，2014）的bamCoverage转换为BigWig格式，并通过每百万映射读取数的RPKM进行归一化，以便在IGV（v2.10.0）（Thorvaldsdóttir等人，2013）中可视化。使用multiBamSummary总结全基因组读取分布，并使用plotCorrelation计算和可视化样本间的相关性。为了生成代表性的信号轨迹以进行可视化和聚合信号比较，使用SAMtools（Danecek等人，2021）合并和索引每个条件下的生物重复样本。使用MACS2（v2.2.8）（Zhang等人，2008）分两阶段进行峰值调用：第一阶段使用参数：–nomodel –shift -37 –extsize 73 -g 1065370000，第二阶段使用估计的平均片段大小，参数：-q 0.05 -g mm -B –nomodel –extsize 260。使用MACS2 bdgdiff和参数-g 60 -l 260获得潜在的差异可访问性区域（候选DARs）。我们注意到，这种聚合信号方法没有明确模拟重复样本间的生物学变异性，因此在这里用于探索性识别候选区域。最后，使用ChIPseeker（v1.38.0，参数：tssRegion = c(-2000, 2000), addFlankGeneInfo=TRUE, flankDistance=5000）（Yu等人，2015）对候选DARs进行注释。

整合分析以识别关键基因
使用UCSC基因组浏览器（Perez等人，2024）批量检索鉴定出的DEGs的基因组坐标，并生成包含其启动子区域的BED文件。然后使用BEDTools（v2.31.0）（Quinlan和Hall，2010）识别与这些启动子区域重叠的候选差异染色质可访问性峰值。其启动子显示出差异可访问性的DEGs被视为候选基因。随后通过整合功能富集和PPI网络分析的结果来选择关键基因。

潜在转录因子的预测
使用KOBAS（Bu等人，2021）对DEGs（p < 0.05，|log?(FC)| ≥ 0）进行GO富集分析。在“DNA结合转录因子活性，RNA聚合酶II特异性”（GO:0000981，padj = 0.0424；表S11）术语中富集的DEGs被定义为差异表达的转录因子（TFs）（表S8）。使用JASPAR数据库（Rauluseviciute等人，2023）分析先前识别出的峰值区域内的TF结合基序（MICA，chr16:2,068,675-2,068,985；MHCIA6，chr16:2,044,792-2,045,249；MHCIY，chr16:1,915,878-1,916,593）的富集情况。然后将得到的候选TFs集与它们的鸡同源物通过UniProt（The UniProt Consortium，2025）对齐，选择序列同源性大于79%的TFs作为关键基因的潜在调节因子。枢纽基因的功能富集分析。
**类别描述：**
上调基因
下调基因

**GO-BP：**
免疫反应：MHCIA6, HLA-F10
AL3, MICA, MHCIY

**内肽酶活性的负调控：**
OVAL, OVALY, KNG1, C4A

**对皮质类固醇的反应：**
OVAL, BPIFB2

**通过抗菌肽介导的体液免疫反应：**
CATH1, TF

**对革兰氏阴性细菌的防御反应：**
CATH1, TFG

**细胞外空间：**
OVALY, KNG1, MHCIA6, C4A, OVAL, HLA-F10AL3, MICA, MHCIY

**DNA修复：**
DBH, GAL, CATH1, TF, CG-16

**质膜的重要组分：**
OVALY, MHCIA6, MHCIY, CYP3A5, HLA-F10AL3, MICA

**跨膜转运：**
DBH, SLC17A8, EGF, CD247, AGTR1

**质膜相关：**
EGF, CD247, AGTR1

**运动活动：**
MYH1B, MYH1D, MYH1F

**丝氨酸型内肽酶活性：**
OVAL, OVALY

**紧密连接：**
MYH1B, MYH1D, MYH1F

**酪氨酸代谢：**
DDC, DBH, TH

**神经活性配体-受体相互作用：**
KNG1, AGTR1

**血管平滑肌收缩：**
ACTG2, AGTR1

**黏附相关：**
VTN, EGF

**候选差异可及染色质区域的鉴定与注释：**
为了从表观基因组角度探究产蛋量变化的原因，我们使用实验组中产蛋量最高和最低的两只鸡的卵巢基质组织构建了ATAC-seq文库。数据质量统计显示，每个样本产生了超过3580万个清晰的可读序列，映射到基因组的比例超过98%（表S6）。
所有样本之间的Spearman相关性分析显示，两个HEL样本之间的相关系数为0.80，两个LEL样本之间的相关系数为0.83。此外，通过聚类分析可以明显区分这两个组（图4A）。

**下载：**
- 下载高分辨率图像（775KB）
- 下载全尺寸图像

**图4. **染色质可及性变异区域的鉴定与注释。**
(A) 不同样本间染色质可及性的Spearman相关性。每个单元格的颜色强度和数字表示相关系数。
(B) HEL组和LEL组之间候选DARs的文氏图。数字表示组内独有的或共有的候选DARs。
(C) 候选DARs的全基因组分布。HEL组和LEL组的特定候选DARs分别用红色和蓝色表示。每一行代表一条染色体；x轴表示基因组位置，y轴表示峰值强度。
(D) 候选DARs的基因组注释。x轴表示每个注释类别中的峰值百分比。彩色条形代表与特定基因组特征相关的候选DARs的比例。
(E) 候选DARs相对于TSS的分布。x轴表示候选DARs与TSS之间的距离（上游和 downstream）。

**进一步探讨卵巢表观基因组差异对产蛋量的影响：**
我们鉴定了HEL组和LEL组之间染色质的候选DARs。HEL组检测到17,014个，LEL组检测到16,377个染色质可及区域，分别包含1,065个和428个候选DARs（图4B）。候选DARs在染色体上分布均匀，无明显偏差（图4C）。这些候选DARs的基因组注释（图4D）显示，HEL组在启动子区域（距离转录起始位点（TSS）2 kb以内）的比例（9.86%）高于LEL组（5.84%）。内含子区的候选DARs在HEL组也更为丰富（10.51%），而在LEL组中较少（3.04%）。此外，在HEL组中，小部分候选DARs位于外显子（0.38%）、UTR区（0.38%）以及TSS下游300 bp处（0.09%）。两组中其余的候选DARs位于基因间区域。分析转录因子结合位点（TFBL）相对于TSS的分布（图4E）显示，HEL组呈平衡分布，而LEL组则偏下游。

**关键基因的鉴定与转录因子的筛选：**
为了探究差异基因表达的潜在机制，我们鉴定了其启动子区域内染色质可及性发生变化的DEGs（DEGs）。这揭示了18个含有候选DARs的DEGs（图4A；表S7）。对这些DEGs的功能富集分析显示，RFT1、MICA、MHCIA6和MHCIY主要与“免疫反应”、“细胞外空间”和“膜的重要组成部分”等术语相关，并且它们在HEL组中均上调。进一步的PPI网络分析将这些基因中的MICA、MHCIA6和MHCIY确定为前20%的枢纽基因，因此将它们定义为关键候选基因。它们的启动子候选DARs显示在图4B和表2中。

**表2. 关键基因列表。**

**基因olog2FC（RNA-seq）**
**p值（RNA-seq）**
**染色体**
**峰值起始**
**峰值结束**
MICA：4.67, 1.05, E-04, 162,068, 675, 2,068, 985
MHCIA6：5.39, 1.21, E-04, 162,044, 792, 2,045, 249
MHCIY：3.27, 2.13, E-03, 161, 915, 878, 1,916, 593

**阐明差异染色质可及性相关关键基因表达的机制：**
我们对p < 0.05且|log?FC| ≥ 0.0的基因进行了GO富集分析。这一结果揭示了在“DNA结合转录因子活性，RNA聚合酶II特异性”（GO:0000981，padj = 0.040）术语中富集的基因，这些基因被认为是差异表达的转录因子（表S8）。这三种关键基因的启动子内的候选DARs被筛查潜在的转录因子结合位点，预测出七个具有高结合亲和力的转录因子（相对分数>90%）：HOXB6、E2F1、GATA3、TFAP2A、WT1和KLF10（表S9）。鉴于JASPAR中的转录因子结合模型主要基于人或小鼠数据，我们评估了这些转录因子与其鸡类同源物的蛋白序列同源性。HOXB6（79.2%）、GATA3（92.0%）、TFAP2A（90.8%）和WT1（86.4%）显示出高度保守性，而其余转录因子的同源性较低（<70.0%）（图S3）。因此，我们认为HOXB6、GATA3、TFAP2A和WT1可能是上调关键基因表达的潜在调节因子（图4C）。这些转录因子可能通过结合其启动子区域内的特定基序来增强关键基因的转录起始。

**图5. **关键基因的定义与潜在转录因子的预测。**
(A) 候选基因筛选的文氏图。DEG启动子（RNA-seq）与启动子注释的候选DARs（ATAC-seq）的交集。
(B) 关键基因启动子区域的染色质可及性剖面。灰色阴影区域突出显示候选DARs。
(C) 预测与关键基因启动子结合的转录因子。

**讨论：**
产蛋后期是影响鸡整体产蛋量的关键时期。在本研究中，我们整合了RNA-seq和ATAC-seq，研究了高产蛋量和低产蛋量文昌鸡在这一阶段的卵巢中的基因表达和染色质可及性差异（图6）。我们鉴定了三个关键基因，包括MICA、MHCIA6和MHCIY。这些基因都属于非经典MHC-Ⅰ（MHC-Ⅰb）家族，位于鸡的16号染色体的MHC-Y区域。这类基因具有高度多态性和异种免疫原性，广泛表达在细胞表面，可能介导特定的免疫反应（Miller和Taylor，2016；Zhang等人，2020）。特别是MHC-Ⅰb分子是T细胞发育、分化和整体T细胞介导的免疫反应的关键调节因子（Won等人，2023；Kim等人，2024）。

**与以往研究的共识一致：**
先前的研究已经表明MHC基因与鸡的产蛋量有关（Ouyang等人，2000）。最近的研究进一步表明，MHC I基因的上调表达可能提高55周龄母鸡的产蛋量和免疫功能（Geng等人，2026）。总体而言，我们的结果支持卵巢免疫稳态在生殖表现中起关键作用的观点。卵巢包含多种免疫细胞，包括单核细胞、巨噬细胞、T细胞和自然杀伤细胞（Shen等人，2023）。尽管这些细胞仅占卵巢细胞总数的很小一部分，但它们通过与卵泡池的细胞因子介导的相互作用建立了必要的免疫微环境（Cacciottola等人，2025）。随着年龄的增长，卵巢中的免疫浸润增加，并与MHC-I通路的上调相关。这一过程可能增强T细胞与邻近细胞之间的相互作用，从而重塑免疫微环境（Wei等人，2026）。在产蛋后期，持续的排卵导致卵巢内逐渐积累炎症反应（Zhang等人，2022）。这种慢性炎症与卵巢T细胞数量减少、间质纤维化增加、卵泡储备耗竭和卵泡闭锁率升高有关，共同影响卵巢的免疫功能（Barua和Yoshimura，1999；Wang等人，2020；Zeng等人，2024）。我们的发现表明，MHC-Ib分子可能有助于缓解这种慢性炎症并维持卵巢免疫稳态，从而支持持续的产蛋能力。我们的ATAC-seq数据进一步支持了这一观点，这些数据揭示了高产蛋量和低产蛋量母鸡在这些关键MHC-Ib基因启动子区域内的染色质可及性差异，表明它们的转录调控受到局部染色质结构变化的影响。

**其他相关基因：**
除这些免疫相关基因外，我们的综合分析还确定了四个可能的调节因子，包括HOXB6、GATA3、TFAP2A和WT1。HOXB6是一种ANTP类同源框蛋白，通过结合增强子参与细胞分化（Yang等人，2024），并与女性生殖系统病理有关，包括卵巢肿瘤（Stone等人，2003）、卵巢癌（Kar等人，2015）和子宫内膜异位症（Geng等人，2022）。GATA3是一种锌指转录因子，对T细胞发育和炎症反应至关重要（Popp等人，2017；Bacha等人，2025），还调节卵巢中的颗粒细胞凋亡和自噬（Luo等人，2023）。此外，GATA3可以与TFAP2A合作共同调节染色质可及性和靶基因表达（Liu等人，2023）。TFAP2A是一种TFAP2家族转录因子，在卵母细胞成熟调节中起重要作用（Lin等人，2022）。WT1是另一种锌指转录因子，参与细胞增殖、凋亡、黏附和细胞间信号传导等多种过程（Zhu等人，2024），对颗粒细胞发育至关重要（Cen等人，2020）。在鸡中，WT1的表达在卵泡发育阶段特异，其在未成熟卵泡中表达较高，在排卵前显著下调（Chun等人，1999；Tang等人，2022）。因此，我们在HEL组观察到的WT1表达上调可能表明有更多的未成熟卵泡，反映了更大的产蛋潜力。重要的是，我们的ATAC-seq分析在这些转录因子结合位点附近发现了差异可及区域，表明染色质重塑可能影响控制这三个关键免疫相关基因表达的调控网络。

**多组学方法的局限性：**
尽管我们的DEG鉴定基于名义p值阈值（p < 0.05）且未进行假发现率（FDR）校正，但使用了edgeR包得到确切的p值（Rapaport等人，2013）。尽管我们的鉴定方法得到了RT-qPCR验证的支持，并在相关研究中得到了广泛应用（Chen等人，2024a；Sun等人，2024；Zhou等人，2024；Huang等人，2025；Liu等人，2025；Duan等人，2026；Li等人，2026），但我们意识到这可能会增加假阳性发现的风险。不同产蛋组之间的染色质可及性差异仅基于探索性分析，未考虑重复样本之间的生物学变异。此外，本研究中确定的关键基因和潜在转录因子仅基于相关性证据，缺乏功能验证。未来的研究应优先应用更严格的标准（例如FDR）进行DEG鉴定，以获得更保守的估计，并采用考虑重复性的统计框架（例如DESeq2或edgeR）来验证这些候选DARs并估计真实的生物学变异。对于关键基因和潜在转录因子，可以采用染色质免疫沉淀和电泳迁移率测定来确认它们的相互作用，并进行增益或失功能实验以阐明关键基因在下游调控途径中的作用。此外，整合其他表观基因组和三维基因组技术将有助于更全面地理解产蛋后期的卵巢调控机制，从而为文昌鸡的精准育种提供可靠的目标。

**结论：**
总之，本研究表明，卵巢免疫稳态在维持鸡产蛋后期期间起着关键作用。通过整合RNA-seq和ATAC-seq，我们鉴定了与染色质可及性变化相关的关键免疫相关基因，并揭示了可能调控它们的转录因子。这种多组学方法使得单一组学分析无法解决的区域性调控网络可以得到系统性的描述。我们的研究结果为文昌鸡产蛋性能的遗传改良提供了潜在的分子靶点，并有助于更深入地理解影响禽类卵巢功能和生殖衰老的调控机制。
作者贡献声明：
杨宏瑞：撰写－审稿与编辑、撰写－初稿、数据可视化、数据验证、数据管理。
王雷：资料准备、实验设计。
荣莉：数据可视化、统计分析。
南久红：统计分析。
蔡金萍：实验设计。
赵云霞：项目指导、资金争取。
朱桂玉：实验方法设计、概念构思。
李世军：撰写－审稿与编辑、项目监督、实验方法设计、资金争取、概念构思。