A total of 200 chickens were obtained from the Poultry Institute of the Chinese Academy of Agricultural Sciences (Yangzhou, China), including 100 Dagu chickens and 100 Wenchang chickens. The two chicken breeds were both from gene pool preservation populations. After hatching, all chicks were housed in an environmentally controlled room until 7 days of age and then transferred to sterilized isolation ventilated cages (IVCs) (IPQ-type 3 negative pressure isolator). The chicks were randomly divided into two groups: noninfected (25 individuals) and Salmonella typhimurium (ST)-infected (75 individuals) in each breed. These birds of the ST-infected group were randomly selected for Salmonella infection at 14 and 28 days of age. After Salmonella infection, the birds of each group were assigned to one IVC. The chicks received ad libitum access to feed and water throughout the experiment.

Salmonella typhimurium 21484 (China Industrial Microbial Culture Preservation Center, Beijing, China) was as infectious strain in this experiment. Before the infection assay, 100 µl cryopreserved bacterial solution was thawed at room temperature and then inoculated into 1000 ml Luria Bertani broth at 37°C with agitation (150 rpm) and cultured for 12 h. Then, 100 µl cultured bacterial liquid was inoculated into 1000 ml Luria Bertani broth and cultivated under the same conditions until the bacterial solution was recovered three times to achieve the highest activity. After centrifugation, the bacterial concentrate was resuspended using sterile phosphate-buffered saline (PBS). The colony-forming units (CFU) of Salmonella was determined by plating serial dilutions in triplicate on agar (37°C, overnight). In infection assay, the chicks from the ST-infected group were infected with 1.5×10¹³ CFU of ST/ml of PBS orally at 14 days old and 28 days old. The chicks from the noninfected group were given same volume of sterile PBS.

The peripheral blood white blood cell count method was performed as described previously (18). Briefly, fresh blood was collected from the wing vein and 10 µl of which smeared on microscopic glass slides for each bird. Blood smears were air-dried and then stained with Wright’s Giemsa solution (G1020, Solarbio, Beijing) according to the manufacturer’s instructions. According to schematic diagram, two hundred leukocytes (heterophils, lymphocytes, and monocytes) were counted using a Leica DM500 microscope with a magnification of 1000× immersion oil. The percentage of monocytes was equal to the number of macrophages to the total number of leukocytes (heterophils, lymphocytes, and monocytes). H/L was the ratio of heterophils to lymphocytes.

The concentrations of IL-17, IL-23 and G-CSF in the serum were measured using enzyme- linked immunosorbent assay (ELISA) kits according to the manufacturer’s instructions (Bioswamp). The assay was performed in duplicate and included ten to eleven non-Salmonella infected and Salmonella infected chickens from Dagu chickens and Wenchang chickens at 24 hours after infection of 28 days old. In brief, optical density (OD) was measured by a microplate reader, standard curve was generated by dilute standard, the concentrations of IL-17, IL-23, and G-CSF in the serum were calculated according to the standard curves.

The primary monocyte-derived macrophages of Wenchang chickens and Dagu chickens were extracted and cultured according to previous studies (19–21). Briefly, peripheral blood mononuclear cells (PBMCs) were isolated from fresh blood using chicken peripheral blood lymphocyte isolation kit (P8740, Beijing Solarbio Science & Technology Co., Ltd.) . After 6~8 h of static incubation in petri dish, adherent cells were washed twice with PBS to remove thrombocytes, nonadherent lymphocytes and other semi adherent cells. These adherent cells were primarily chicken monocytes. Subsequently, fresh Roswell Park Memorial Institute (RPMI)-1640 medium with 10% fetal bovine serum (FBS, Gibco), 5% chicken serum, 1% penicillin and streptomycin were added to the remaining monocytes. Chicken monocytes were then cultured in 12-well cell culture plates to generate macrophage differentiation for subsequent experiments. Treatment of all individuals during the experiment was carried out under consistent conditions.

Assessment of phagocytic activity in macrophages was performed using a flow cytometry-based method. Briefly, 0.5 ml RPMI-1640 (containing 10% fetal bovine serum, 5% chicken serum, 1% penicillin and streptomycin) was added to each well covering macrophages in cell culture plates. Macrophages were stimulated with LPS (500 ng/ml) for 6 h, and then 4 µl of latex beads (Sigma-Aldrich, cat. L3030-1ML, 2.0 µm mean particle size; St. Louis, MO, USA) were mixed, and the samples were incubated for 2 h at 37°C. After incubation, the supernatant was removed (including nonadherent cells and excess fluorescent microspheres), and the cells were washed twice with cold PBS to stop phagocytosis. Next, adherent cells were trypsinized and transferred from the 12-well plates to test tubes containing 0.5 ml PBS. The samples were filtered through a 100-μm filter and analyzed by flow cytometry.

The cell populations were distinguished based on forward scatter (FSC), side scatter characteristics (SSC) and GFP-positive cells. First, we set up the FSC-SSC 2D scatter plot, the macrophage gate was used to define the analyzed macrophage population by adjusting the voltage values of FSC and SSC, and the interference of other nucleated cells and cell debris was excluded to the maximum extent ( Supplementary Figure S1 ). The fluorescence intensity of macrophages was detected by the GFP pathway, the fluorescence pathway of light emitted by fluorescent microspheres, obtaining 10000 macrophages per sample. Data are displayed in two-dimensional scatter plots and histograms. In the two-dimensional scatter diagram, the macrophage population without phagocytosis of fluorescent microspheres and the macrophage population phagocytosis of fluorescent microspheres could be delineated by setting gates. In the histogram, the macrophage population that did not phagocytose fluorescent microspheres and the macrophage population that phagocytosed fluorescent microspheres could be calibrated by using a scale, and the proportion of macrophages that did not phagocytose fluorescent microspheres and macrophages that phagocytosed different numbers of fluorescent microspheres was analyzed. The phagocytosis percentage and phagocytosis index were calculated according to the following formula:

To characterize the expression differences of related genes during phagocytosis, a total of 16 macrophage samples were collected from Wenchang chickens and Dagu chickens (n=8 each breed). Macrophages of non-infection group were subjected to RNA-seq during phagocytosis testing at 28 days old. Total RNA was extracted using TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA). Using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA) to measure RNA concentration and purity. Using RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) to assessed RNA integrity.

1 μg RNA per sample was used as input material. Using poly-T oligo-attached magnetic beads to purify mRNA from total RNA. Sequencing libraries were constructed using the NEBNext UltraTM RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). Finally, polymerase chain reaction (PCR) products were purified (AMPure XP system). Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA) was used to assessed the library quality. Library preparations were sequenced on an Illumina platform, and paired-end reads were generated. Approximately 6.35 Gigabases of clean data were obtained for each sample.

In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N sequences and low-quality reads from the raw data. The Q20, Q30, GC content and sequence duplication level of the clean data were calculated. All downstream analyses were based on high quality clean data. The clean data were further analysis with a bioinformatic pipeline tool, BMKCloud (www.biocloud.net).

Quantification of gene expression was estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). EBSeq performed differential expression analysis between two groups. False discovery rate (FDR) < 0.05 and fold change ≥1.5 were set as the thresholds for significantly differentially expressed genes (DEGs).

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed by KOBAS 3.0 (22) to explore the function of differentially expressed genes (DEGs) and hub genes. P < 0.05 was the threshold for significant enrichment.

All 16 transcriptomes data of macrophages were used to perform weighted gene co-expression network analysis. FPKM was utilized to construct a gene expression matrix (FPKM ≥ 1, variation of FPKM: cv ≥ 0.5). Weighted correlation network analysis (WGCNA) in the R package with the ‘unsigned’ weighted network type was carried out (23). Soft thresholding power option was set to 15, topology fit index R² was first greater than 0.8, minimum module size was limited to 30. A minimum height of 0.22 was chosen to merge highly co-expressed modules. P < 0.05 was defined as the significance threshold. Module eigengenes (MEs) were referred to as ‘‘hub genes”, dependent on high intramodular connectivity values (|kME|≥0.8).

The HD11 cells were maintained in RPMI-1640 medium supplemented with 10% FBS, 5% chicken serum, 1% penicillin, 0.1% 2-Mercaptoethanol, 1% nonessential amino acid and 0.24% HEPES. The homogenized cell suspensions were seeded into 24-well plates at 5 × 10⁵ cells/well for a 24-well plate and maintained in humidified incubator at 37°C with 5% CO₂. Cells were stimulated with lipopolysaccharide (LPS) at 100ng/ml for 6 hours and 12 hours.

The total RNA of cells was extracted by Trizol reagent. RNA (1500ng) was reverse transcribed by cDNA synthesis kit (TIANGEN) for quantitative real time PCR. Primers were designed according to chicken coding region sequences and synthesized by The Beijing Genomics Institute (BGI), which are shown in Table 1 . ACTB was selected as house-keeping gene. Quantitative real-time PCR was performed in triplicate using the Invitrogen PowerUp™ SYBR^® Green Master Mix (ABI). with the following cycle profile: 95 °C for 3 min, 40 cycles of 95 °C for 3 s, annealing temperature for 34 s in the QuantStuio 7 Flex Real-Time PCR System (Waltham, MA, USA). The results were analyzed by 2^−ΔΔCt method (24).

DNA was extracted from the blood samples of 30 Dagu chickens and 30 Wenchang chickens by the conventional phenol-chloroform method. Sequenced on the Illumina HiSeq2500 platform after quality checking and quantification. The Genome Analysis Toolkit (GATK) (25) was employed for preprocessing and single nucleotide polymorphism (SNP) calling. Phylogenetic tree was built based on the weighted neighbor-joining method (26) and visualized using MEGA7.0 software (27) and iTOL software (28). SNPs were pruned based on linkage disequilibrium utilizing PLINK software (Version: 1.90) (29) with the options ‘–indep-pairwise 25 10 0.2 -mind 0.1 –maf 0.01’. Principal component analysis (PCA) was performed using PLINK software. Based on allele frequencies at variable sites, the average fixation index (FST) and nucleotide diversity(pi) were calculated to identify selection regions (–windows–40 kb –step 10 kb). For each window, the average FST values (30)and pi value (31) were assessed between the Dagu chickens and the Wenchang chickens. pi value was log₂-transformed. The windows with the top 5% FST values and log₂(pi-ratio) values were simultaneously shown as candidate outliers.

Student’s t-test analyzed the significance of differences among groups. P < 0.05 (*), P < 0.01 (**) or P < 0.001 (***) implied a statistically significant difference respectively. Data were visualized using GraphPad Prism version 8 (GraphPad Software, San Diego, CA, USA) and R version 4.0.5 and expressed as mean ± standard deviation (SD) unless otherwise indicated.

To assess the level of immune response to infection, 28-day-old Dagu chickens and Wenchang chickens were challenged with Salmonella. Dagu chickens were characterized by a significantly increased percentage of macrophages and heterophils and lymphocytes (H/L) ratio in peripheral blood compared to Wenchang chickens. Additionally, the percentage of macrophages remained largely unchanged, while the H/L ratio increased significantly after Salmonella infection, which was consistent with previous reports (32) ( Figures 1A, B ). Furthermore, we compared interleukin 17 (IL17), interleukin 23 (IL23) and granulocyte colony-stimulating factor (G-CSF) concentrations in serum between Dagu chickens and Wenchang chickens. Unlike the above two indicators, the concentrations of IL-17, IL-23, and G-CSF were significantly higher in Wenchang chickens than in Dagu chickens before and after Salmonella infection ( Figures 1C–E ).

To compare the phagocytic ability of peripheral blood macrophages between Dagu chickens and Wenchang chickens, 100 chicks of each breed were randomly divided into a control group and an infection group. At 14 days of age, the phagocytosis index of macrophages was significantly increased after Salmonella challenge, especially at 6 hours after infection (6 hpi) ( Figure 2A ). A similar trend existed for macrophage phagocytosis rate ( Figure 2B ). Remarkably, whether it was the phagocytosis index or phagocytosis rate of macrophages, Dagu chickens were always higher than Wenchang chickens, which was more significant after Salmonella challenge. Furthermore, at 28 days of age, under non-infection conditions, the phagocytosis index and phagocytosis rate of macrophages in Dagu chickens were significantly increased than Wenchang chickens, while the difference decreased after Salmonella challenge ( Figures 2C, D ). The results showed that the difference in the phagocytic ability of macrophages between Wenchang chickens and Dagu chickens at the early stage depended mainly on the stimulation of the external pathogenic environment. However, with age, differences in phagocytic ability were directly manifested and related to species specificity. Furthermore, the body weight of Dagu chickens was significantly increased compared to the body weight of Wenchang chickens at 14 days old, while the difference disappeared at 28 days old ( Supplementary Figure S2 ). Overall, the advantage of the macrophage phagocytic ability of Dagu chickens was more prominent.

Phagocytosis is an important cellular process that induces antimicrobial responses and modulates adaptive immunity (33). This process is coordinated by multiple gene signaling pathways to activate the internalization machinery (34). Analysis of the transcriptional diversity associated with phagocytosis to reveal the genetic mechanism of the difference in phagocytosis between Dagu chickens and Wenchang chickens should be possible. Based on this premise, we sought to integrated mRNA sequencing data to explore the functional factors involved in phagocytosis from peripheral blood macrophages of 8 Dagu chickens and 8 Wenchang chickens during phagocytosis. PCA showed that the two groups of individuals were distinctly clustered together ( Figure 3A ). 1136 differentially expressed genes (DEGs) were identified between the Dagu chickens and Wenchang chickens, of which Dagu chickens performed 436 upregulated genes and 700 downregulated genes compared to Wenchang chickens ( Figure 3B ; Supplementary Table S1 ). Hierarchical clustering was consistent with PCA based on all DEGs and showed that there were more downregulated than upregulated genes ( Figure 3C ). Functional enrichment analysis of the DEGs showed that many biological processes were related to phagocytosis, including the cell cycle, cytokine-cytokine receptor interaction, ECM-receptor interaction, apoptosis, and phagosome ( Figure 3D ; Supplementary Table S2 ).

Some interesting DEGs (expression fold changes greater than 2 between Dagu chickens and Wenchang chickens) were associated with immune responses, including CX3CR1, IL1R2, CENPF, CCR5, DIAPH3, and IQGAP3 ( Figure 3E ). Analysis of above genes expression patterns revealed that multiple genes were likely involved in chicken macrophage biological functions in immunity.

Weighted gene co-expression network analysis (WGCNA) was constructed to further elucidate biological pathways involved in phagocytosis. First, the best soft threshold of 15 was determined from the scale-free topological model and mean connectivity ( Figure 4A ). Then, using hierarchical clustering of genes based on the topological overlap matrix (TOM), a cluster dendrogram of co-expression network modules was generated ( Figure 4B ). A total of 22 co-expression modules were captured, three of which (MEdarkorange2, MEivory, and MElightyellow) were significantly correlated with the phagocytosis index, and four of which (MEcyan, MEdarkorange2, MEivory, and MElightyellow) were significantly correlated with the phagocytosis rate (Pearson correlation, P < 0.05). Notably, the phagocytosis index was strongly positively correlated with the phagocytosis rate (R² = 0.77) ( Supplementary Figure S3 ). The MEdarkorange2, MEivory, and MElightyellow modules were jointly identified modules, that were significantly negatively correlated with the phagocytosis index and phagocytosis rate ( Figure 4C ).

Functional enrichment of the genes in these significant modules revealed close linkages with phagocytosis. Three negatively correlated modules, MEdarkorange2, MEivory, and MElightyellow, together showed enrichment for metabolic pathways (most genes involved in this pathway) and oxidative phosphorylation ( Figures 5A–C ; Supplementary Table S3 ). Different modules were implicated in several different pathways: cellular senescence and phagosome for MElightyellow; MAPK signaling pathway and focal adhesion for MEivory; ubiquitin mediated proteolysis and cell cycle for MEdarkorange2.

A total of 584 protein-coding genes were identified as hub genes in three significant modules ( Supplementary Table S4 ). Furthermore, 2 MElightyellow hub genes, 1 MEivory hub gene and 19 MEdarkorange2 hub genes were differentially expressed. These genes were all significantly downregulated in Dagu chickens compared to Wenchang chickens ( Figure 6 ). Notably, the MElightyellow hub gene S100A12 is a member of the S100 family of proteins containing 2 EF-hand calcium-binding motifs involved in the innate immune response, including monocyte chemotaxis, neutrophil chemotaxis and positive regulation of NF-kappaB transcription factor activity (35). The MEdarkorange2 hub gene CUEDC2 acts the negative regulation of macrophage cytokine production in the inflammatory response (36). These findings suggested that the mutual regulation of multiple genes and pathways affected the immune response of chicken macrophages.

To verify whether the DEGs and hub genes screened in the above results were associated with inflammation, we detected the expression of these genes in macrophages after LPS stimulation in HD11 cells ( Figure 7 ). The rapid increase in the expression of IL1β indicated that the construction of the inflammatory model was successful. The expression of IL1R2 gene showed the similar trend as that of IL1β gene, which first increased and then decreased within 12 hours after LPS stimulation. Furthermore, CX3CR1, CENPF, AKR1E2 and H2AFZ genes expressions increased significantly, while CCR5 and DIAPH3 genes expressions decreased significantly at 12 hours after LPS stimulation. Notably, the expressions of CCR5, S100A12, IQGAP3 and ADH4 genes differed significantly within 6 hours after LPS stimulation, indicating that these genes correspondingly antigen stimulated performed a faster immune response.

Comparative population genomics evaluated the population differentiation based on genomes of 30 Dagu chickens and 30 Wenchang chickens. A phylogenetic tree and PCA clearly showed population differences between the Dagu chickens and Wenchang chickens ( Supplementary Figure S4 ). The top 5% windows based on the value of FST and log₂PI_(WC/DG) were defined as strong repeatable signals of selection and visualized as manhattan plots ( Figure 8A ). In particular, 26 DEGs and 13 hub genes intersected with selected regions ( Figure 8B ). However, only the H2AFZ gene was identified among the DEGs and hub genes ( Figure 8C ). H2AFZ is a highly conserved gene encoding H2A.Z.1, mainly regulates cell cycle signaling and DNA replication (37).