Immunogenetic Basis of Chicken’s Heterophil/Lymphocyte Ratio Revealed by Genome-wide Indel Variants Analysis

doi:10.21203/rs.3.rs-1448551/v1

Download PDF

Research Article

Immunogenetic Basis of Chicken’s Heterophil/Lymphocyte Ratio Revealed by Genome-wide Indel Variants Analysis

https://doi.org/10.21203/rs.3.rs-1448551/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Enhancing host immunity is an effective way to reduce morbidity in chickens. Heterophil/lymphocyte (H/L) ratio is associated with host disease resistance in birds. Chickens with different H/L ratio levels show different disease resistances, caused by a difference in the functional capacity of heterophils. In this study, an H/L directional breeding chicken line was constructed and a combined statistical analysis was used to explore the candidate genes and pathways related to H/L ratio based on whole-genome sequencing of 368 individuals, divided into a non-selection line (NS) and selection lines including G5, G9 and F₂ of G9.

Results

Heterophils isolated from the H/L selection line had stronger phagocytosis and oxidative burst function than non-selection line. The 22.44 Mb genomic regions which annotated 300 protein-coding genes were selected in the genome of G9 compared to NS based on a genome-wide selective sweep. Several selective regions were identified containing genes like interferon induced with helicase C domain 1 (IFIH1) and moesin (MSN) associated with the intracellular receptor signaling pathway, C-C motif chemokine receptor 6(CCR6)、dipeptidyl peptidase 4(DPP4) and hemolytic complement (HC) associated with the negative regulation of leukocyte chemotaxis and tight junction protein 1(TJP1) associated with actin cytoskeleton organization. In addition, 45 genome-wide significant indels containing 29 protein-coding genes were also identified as associated with the H/L ratio based on genome-wide association study (GWAS). The expression of protein tyrosine phosphatase non-receptor type 5(PTPN5) and oxysterol binding protein like 5(OSBPL5) were positively correlated with H/L ratio and significantly different between the high and low H/L groups of the Beijing-You chicken. The A/A allelic frequency of indel 5_13108985 (P-value= 3.85E-06) within OSBPL5 gradually increased from the non-selection line to the G5 and G9 and the individuals with A/A exhibited lower H/L ratio.

Conclusions

The phagocytosis and oxidative burst functions were enhanced and 22.44 Mb genomic regions were selected with the directional selection of H/L ratio. PTPN5 and OSBPL5 genes were identified as candidate genes. These findings revealed the complex genetic mechanism of H/L ratio related to immunity and will allow selection for improving chicken immunity based on the H/L ratio.

chicken

H/L ratio

heterophil

indel

selective sweep

GWAS

Avian diseases represent an ongoing threat to poultry industries and human health worldwide and have become more significant with policies to reduce the use of antibiotics in animal production. Innate immunity can make up for the lack of vaccination and medication when chickens are challenged by diseases and the development of resistant chickens through selective breeding could prevent spread. This selection for immune function can be part of the breeding standard of commercially important traits in poultry industry, but is a challenging proposition for commercial breeding companies due to a deficiency of indicators accurately reflecting the immune strength [1].

Counts of relative leucocytes is the simplest but most robust measure to assess animal immune function using blood smears [2] and the heterophil and lymphocyte are the two most numerous immune cell types in avian leucocyte profiles [3, 4]. The heterophil/lymphocyte (H/L) ratio has been shown to be correlated with the strength of innate immune response in poultry [5–7]. The classification of resistant and sensitive lines according to H/L ratio was reported in earlier studies. Al-Murrani used the H/L ratio as a selection criterion for heat stress resistance in chickens and found that the resistant (R) group exceeded the egg production of the sensitive (S) group by 21%, had 4% more egg weight and showed lower progeny mortality at 18·6% compared with 28·4% [8]. Al-Murrani also showed that the H/L ratio was a criterion for selection of resistance to Salmonella typhimurium in chickens [9]. Thiam reported that the H/L ratio was associated with the intestinal barrier and immune response for Salmonella enteritidis clearance and the chickens with low H/L ratio had enhanced intestinal immunity [10]. The high heritability of the H/L ratio supported its use as an indicator for disease resistance in poultry breeding [11].

The baseline H/L ratio of individuals has been shown to be highly stable long-term [12]. But tissue injury due to infection can cause a short-term H/L ratio change, as lymphocytes migrate from the circulating blood to the skin, spleen and lymph nodes, while heterophils migrate from bone marrow into the blood [13, 14], so the rapid increase of H/L ratio in blood is mainly affected by the increase in heterophils. Heterophils play an indispensable role in the innate cellular immunity of the avian host and their biology resembles that of mammalian neutrophils [15]. They are the first line of defense against various invasive pathogens and irritants using toll-like receptors (TLR), Fc and complement receptors and other pathogen recognition receptors[16]. On detection of pathogens, heterophils respond and relay information to other immune cell types through a network of intracellular signaling pathways and the release and response to cytokines and chemokines [17]. Measures of heterophil function include phagocytosis, degranulation, oxidative burst, DNA and histone extracellular matrix (HETs) [18, 19]. Phagocytosis is the main process used by heterophils to kill invading pathogens, immediately followed by degranulation and production of an oxidative burst [20]. Heterophil function is under immune control and it is possible to genetically select for the enhancement of heterophil function [18, 21, 22].

The immune genetic mechanism of H/L ratio has not been fully analyzed, so its use is restricted due to the lack of theoretical information. To determine the genomic basis related to H/L ratio in immune function, this study performed whole-genome resequencing of the H/L ratio selection and non-selection lines constructed from Jingxing yellow chickens. The authors previous studies have shown that H/L ratio selection line had a higher survival rate after Salmonella infection than non-selection line (to be published). The H/L ratio selection line may be able to explain the genetic mechanism that contributes to the more resistant state, so identifying pathways or genes associated with the H/L ratio in poultry provides the opportunity for its wider application in genetic modification and selective breeding.

Phagocytosis and oxidative burst function of heterophils were different between the H/L selection line and non-selection line

Phagocytosis and oxidative burst are two important indicators to measure the microbial killing ability of heterophils. Phagocytosis and the oxidative burst function of heterophils isolated from selection and non-selection lines were compared at seven days post hatch. The phagocytosis was evaluated for the ability to phagocytize FITC-dextran by flow cytometry and showed a higher percentage of positive heterophils isolated from selection line contained FITC-dextran (35.50±2.56%) compared to heterophils isolated from the non-selection line (29.75±2.18%), though this difference was not significant (Figure1A).

Dydrorhodamine123 (DHR123) was used to determine if there was a difference in oxidative burst function of heterophils isolated from the selection and non-selection lines using flow cytometry (Figure1B). Heterophils isolated from the selection line generated significantly (P＜0.05) higher oxidative burst levels than heterophils isolated from non-selection line at 45.50±0.94% and 31.90±1.90%, respectively. The results showed that heterophils isolated from the H/L selection line had stronger phagocytic and oxidative burst function than non-selection lines and indicated the functions of the heterophil could be enhanced though H/L directional breeding.

The population structure differentiated with the directional selection of H/L

There were 1068153 insertion/deletion (indels) polymorphisms identified on 1 to 28 autosomes after filtering using PLINK (version 1.9) software[23] with a call rate more than 90％ and a minor allelic frequency (MAF) more than 0.05. There were 227 individuals that passed filters and quality control for genetic selection analysis. Principle component analysis (PCA) revealed that the G5 and G9 of selection and non-selection lines clustered in separate populations (Figure2A) and the G9 of the selection line clustered further with the non-selection line, indicating that they had greater genetic distance. The G5 of the selection line occupied an intermediate position along the axis and appeared to be more mixed with the G9 of selection line than with the non-selection line. To investigate the genetic relationships among different populations, a phylogenetic tree was constructed using the distance-based neighbor-joining method. The phylogenetic tree (Figure2B) was consistent with the PCA results and supported both the similar genetic makeup of the G5 and G9 of selection line, as well as the presence of two major clusters, consisting of the H/L ratio selection line and non-selection line. The best K=2 was next determined based on the cross-validation error for different numbers of ancestral genetic backgrounds, and the admixture analysis (K=2) was in line with the results of the PCA and phylogenetic tree analysis supporting a two-clade population structure and there was more gene flow between G5 and G9 of the selection line. There was no significant genetic improvement in G9 of the selection line compared with G5 of the selection line due to shorter breeding generations (Figure2C).

Identification of signature of selection between H/L selection line and non-selection line

By positive selection using the genomes of 184 G9 individuals made up of 92 from the selection line and 92 from the non-selection line were scanned for signatures that might underpin their distinctive features. Two complementary statistics were calculated using a sliding-window approach as 40-kb windows sliding in 20-kb steps. The first statistic was the pairwise genetic differentiation Fst index for identifying those genomic regions maximizing allelic frequency differences between the two populations. The second statistic measured variability in nucleotide diversity between the two groups, respectively. The θπ ratios θπ.NS/θπ.G9 were calculated based on the θπ values and then log2-transformed the data. The windows with the top 5% Fst where Fst exceeded 0.06 and log2 θπ ratio exceeded 0.37 simultaneously were considered as candidate regions under H/L directional selection in the selection line, which provided a total of 561 windows encompassing 22.44 Mb genomic regions and 300 protein-coding genes as seen in Figure3 and Additional file 2(Table S1 and S2). These selection candidates showed limited genetic diversity in G9 of the selection line, but not in the non-selection line, suggesting positive selection for H/L-related phenotypes. Overlapped putative genes were significantly functionally enriched in the intracellular receptor signaling pathway (GO:0030522), negative regulation of leukocyte chemotaxis (GO:0002689) and positive regulation of G2/M transition of mitotic cell cycle (GO:0010971), which may play important roles during the process of cellular immune function in chickens, seen in Figure4 and Additional file 2(Table S3).

Identification of genes associated with H/L ratio

For genome-wide association analysis of the H/L ratio, 249 samples and 983198 indels were kept after filtering using the same parameters as genetic selection analysis. Descriptive statistics of the phenotypic measurements of the H/L ratio traits used for gnome-wide association study (GWAS) are provided in Additional file 2 (Table S4). This analysis accounted for the population structure by adjustment of PCA covariates and genetic relatedness among all pairs of individuals as seen in Additional file 1FigureS1. Manhattan and quantile-quantile (Q-Q) plots are shown in Figure 5A and B. The genomic inflation factor (GIF) was 1.02, which suggested the population structure was well controlled. A total of 45 genome-wide significant indels associated with H/L were identified, of which three indels reached the 5% genome-wide significance (5_8074623, P-value=7.06E-08; 8_16086712, P-value=1.64E-07; 5_8015233, P-value=2.16E-07). All genome-wide significant indels were annotated to 29 protein-coding genes and the detailed information was shown in Additional file 2 (Table S5). The most significant indel 5_8074623 (P-value= 7.06E-08) was an intron1-2 of TEA domain transcription factor 1(TEAD1) and accounted for 11.08% genetic variance for H/L ratio.

To further identify candidate genes affecting H/L in other breeds of chicken, the Beijing-You chicken were selected as the verification population. The correlation between GWAS-related genes expression of the liver and H/L ratio in eight 28-day-old of Beijing-You chickens was compared. The Pearson correlation coefficient heatmap showed that the expression of PCDH17 was negatively correlated with H/L(r=-0.65), while the expression of PTPN5 (r=0.75) and OSBPL5 (r=0.89) were positively correlated with H/L as seen in Figure6A. Chickens were divided into two groups according to H/L value, where three chickens were deemed H/L-high (0.46±0.16) and three chickens as H/L-low (0.06±0.02) as seen in Figure6B. Except for lactate dehydrogenase A (LDHA)and transmembrane protein 86A (TMEM86A), the expression levels of the other 27 genes in the liver of Beijing-You chickens were low and the expression of sodium/potassium transporting ATPase interacting 2(NKAIN2), ras association domain family member 10 (RASSF10) and spalt like transcription factor 1(SALL1)genes were not detected. The differences in expression levels of PCDH17, PTPN5 and OSBPL5 in H/L-high and H/L-low groups were compared, and PTPN5 and OSBPL5 genes expressions were significantly different (P≤0.05)( Figure6C).

Four indels, 5_12975134(9.48kb downstream of PTPN5, P-value=6.48E-06), 5_13108985(P-value= 3.85E-06), 5_13151117(P-value=4.95E-06) and 5_13170170(P-value=7.23E-06) intron1-2 of OSBPL5, were identified as candidate variants. Allelic frequency analysis in the non-selection line, G5 and G9 showed the 5_13108985 variant allelic frequency gradually increased or decreased with the selection of H/L in different generations (Figure7A). Jingxing-yellow chickens with homozygous A/A for 5_13108985 variant showed significantly lower H/L phenotype than individuals with heterozygote A/ATCT and homozygous ATCT/ATCT (P<0.05), respectively (Figure7B).

In this study, cell biology assay proved directional selection of H/L ratio enhanced the phagocytosis and oxidative burst functions of heterophils in chickens. To reveal the immunogenetic basis of the H/L ratio, combined statistical analysis was used to search for putative genes by genomes of 360 chickens from different populations termed NS, G5, G9 and F₂ population of G9. Extensive genome data set allowed the identification of 3970941 indel variants that helped to unveil the genetic evolution pattern and phenotype related genes of the H/L ratio for artificial breeding in chickens.

The gene expression and function of the heterophil is under genetic control. Redmond and colleagues proved heterophils from commercially selected and non-selected genetic lines expressed cytokines differently in vitro exposure to Salmonella enteritidis, where functional phenotypes of phagocytosis and bacterial killing varied depending on the breed or cross of chicken and the SLC11A1 gene was associated with phagocytosis of Salmonella enteritidis [22]. Swaggerty demonstrated genes associated with fast and slow feathering located on the Z sex chromosome were associated with the biochemical killing mechanism of degranulation and oxidative burst of heterophil in avian species [24]. This study found the phagocytosis and oxidative burst of heterophils isolated from selection lines performed more powerfully compared to the non-selection line, indicating heterophil function related genes were selected under the H/L ratio selection pressure. To further explore the relevant genetic basis affecting the function of heterophil, combined statistical analysis was performed based on whole genome resequencing data.

Interested candidate genes are listed in Additional file 2(Table S2) based on genome-wide selective sweep. These candidate genes were widely involved in differential cell molecular processes, which implied complex genetic mechanisms of immune evolution in the process of H/L ratio directional breeding. The intracellular receptor signaling pathway GO:0030522 was the most significant pathway enriched by GO and KEGG pathway enrichment analysis and there were 13 selected genes regulating this pathway. Interferon induced with helicase C domain 1, also known as IFIH1 or MDA5 engaged in innate immune response and response to virus. In chickens, MDA5 recognized infectious bursal disease virus infection and this interaction resulted in the activation of chMDA5-related innate immune genes and upregulation of chicken MHC class I [25]. The MDA5 signaling pathway and innate immune cytokines were induced after avian reovirus infection, which would contribute to the avian reovirus-host interaction, especially at the early infection stage [26]. Moesin (MSN) is the most abundant ERM, the plasma membrane (PM) organizers in leukocytes and is localized to filopodia and other membranous protrusions that are important for cell-cell recognition and signaling and for cell movement, playing a significant role in leukocyte polarization, migration and intercellular adhesion [27].

Leukocyte chemotaxis is a crucial feature of the innate immune response and the prerequisite for the elimination of pathogens. Several genes were found to participate in the negative regulation of leukocyte chemotaxis (GO:0002689), CCR6 encoded chemokine (C-C motif) receptor 6 and enabled C-C chemokine binding activity and C-C chemokine receptor activity [28]. Dipeptidyl peptidase 4 (DPP4) engaged in the negative regulation of neutrophil chemotaxis [29]. HC encoded a component of the complement system, a part of the innate immune system that played a significant role in inflammation, host homeostasis, and host defense against pathogens and was reported to act upstream of or within neutrophil homeostasis [30]. The migration of heterophils to inflammatory sites is an important early stage of heterophil signal recognition of microorganisms or other antigens. Recognition of pathogen-associated molecular patterns (PAMPs) by pattern recognition receptors (PRRs) on the heterophil triggers the process of phagocytosis and secretion of cytokines and chemokines, which directs acquired immunity to elicit the appropriate host response to the microbe. Among known PRRs, the Toll-like receptors (TLRs) are the most important receptor of the heterophil in innate immune response to bacterial infection, of which TLRs 1/6/10, TLR 2 (types 1 and 2), TLR 3, TLR 4, TLR 5, and TLR 7 have been identified in the chicken heterophil [31], and TLR5 showed limited genetic diversity in G9, but not in non-selection line, suggesting positive selection for H/L-related phenotypes. The actin cytoskeleton organization signaling pathway played a key role in the cross endothelial migration and phagocytosis by heterophils and tight junction protein 1 (TJP1) was identified associated with the pathway, where it regulated the movement of ions and macromolecules between endothelial and epithelial cells [32]. Abdala-Valencia demonstrated that overexpression of TJP1 blocked leukocyte trans-endothelial migration [33].

Genome-wide association analysis was performed to identify the candidate genes of the H/L ratio. A previous study reported caspase recruitment domain family member 11 (CARD11), biogenesis of ribosomes BRX1 (BRIX1) and BTG3 associated nuclear protein (BANP) were associated with H/L ratio and played a role in its cellular functions [34]. In this study, a total of 45 genome-wide significant indels associated with H/L were identified and annotated to 29 protein-coding genes. The most significant indel 5_8074623 (P-value = 7.06E-08) was in intron1-2 of TEA domain family member 1 (TEAD1) involved in positive regulation of cell growth [35]. In immunity, TEAD1 directly bounded genomic regions adjacent to TLR1, TLR2, TLR3, TLR4, TLR5, TLR6, TLR7, with TLR9 and YAP/TEAD1 complex acting as a default repressor of cardiac toll-like receptor genes [36]. The activation of the TLR induced signaling pathway played a significant role in the antibacterial mechanism of the avian heterophil [31]. The expression level of 29 candidate genes in the liver of the Beijing-You chicken was verified according to grouping of high and low H/L ratio. Only the expression levels of PTPN5 and OSBPL5 were positively correlated with H/L ratio and were differentially expressed in high and low H/L groups. Protein tyrosine phosphatase, non-receptor type 5 (PTPN5) could enable protein tyrosine phosphatase activity, involved in negative regulation of MAP kinase activity [37] and p38 MAPK was reported as a substrate of PTPN5 [37, 38]. Heterophil extracellular traps (HETs) formation, a novel defense of chicken heterophils playing a significant role against pathogen infection, was related to ROS production dependent on the activation of p38 MAPK signaling pathways [39]. The Fc-receptor-mediated degranulation of chicken heterophils was regulated by the phosphorylation of p38 MAPK signaling pathways which was stimulated by the activation of the Syk tyrosine kinase [40]. The above evidence suggested that PTPN5 might play a significant role in the functional regulation of avian heterophil. Another gene, OSBPL5 (also known as ORP5) encoding oxysterol binding protein-like 5 was predicted to enable several functions, including cholesterol binding activity, phosphatidylserine transfer activity and phospholipid binding activity based on NCBI database. The ORP5/8 regulated Ca₂⁺ signaling at specific membrane contact site foci appear to play roles in processes such as mitochondrial respiration and cell proliferation [41], where Ca₂⁺ channels directly sense oxidative bursts and respond with Ca₂⁺ transients adjacent to active mitochondria [42]. The oxidative burst was catalyzed by a multi-protein complex NADPH oxidase existing in the inactive, dissociated state in the resting cells. The NADPH oxidase in the plasma membrane was activated with combinations of cytosolic proteins, Ca₂⁺, calmodulin and protein kinase following stimulation [43].The mechanism involved in activation of NADPH oxidase was conserved between chicken heterophils and mammalian neutrophils. The A/A homozygote individuals of indel 5_13108985 (P-value = 3.85E-06), within intron1-2 of OSBPL5 had a lower H/L ratio and the allelic frequency gradually increased with H/L directional selection (NS-G5-G9). Whether the insertion of three bases ATCT into the A allele was a deleterious mutation requires further verification.

The directional selection of H/L ratio enhanced the phagocytosis and oxidative burst functions of heterophils in chicken. There were several genes involved in intracellular receptor signaling pathways, negative regulation of leukocyte chemotaxis and actin cytoskeleton organization pathways, which were selected during the directional breeding of the H/L ratio. The PTPN5 and OSBPL5 genes were identified as candidate genes for the H/L ratio. which played a significant role in regulating the function of heterophils. These genes and pathways provided evidence for the genetic effect of artificial directional selection on the evolution of the H/L ratio in chicken. Although more research is needed to confirm the results, this study provides a new insight into the research of genetic improvement of immune traits such as the H/L ratio.

Experimental birds

Jingxing-yellow chickens, a domestic broiler breed, were obtained from Changping Base of Chinese Academy of Agricultural Sciences, from which a directional selection line was established based on the 56-day-old peripheral blood H/L ratio of Jingxing-yellow chickens. Thirty roosters and 90 female chickens with low H/L phenotype were selected as parents in each generation, based on the principles that one rooster was mated with three hens, roosters and hens in the same family were not full siblings and hens within a family were not half siblings.

Phenotype determination

The H/L ratio was evaluated at 56 days old in each generation. Fresh blood was extracted from the wing vein and stored in EDTA tubes. One drop of blood was smeared on to a glass slide for each bird. Once air-dried, smears were stained using a Wright’s Giemsa solution (G1020, Solarbio, Beijing) according to manufacturer's instructions for microscopic observation. One hundred leucocytes including lymphocyte, heterophil and monocyte cells were counted in different fields on each slide, using a Leica DM500 microscope with a magnification of 1000× and the H/L ratio was calculated by dividing the number of heterophils by the number of lymphocytes [44].

Detection of phagocytosis and oxidative burst function of heterophils

Heterophils were isolated from the peripheral blood of chicks at seven days post-hatch using the commercially-available chicken heterophils isolation kit (Solar Technology Co., China). Phagocytic ability was determined by flow cytometry by the detection of fluorescein isothiocyanate dextran (FITC-dextran), 40 kDa (Chondrex, USA) using a previously published method with modification [45]. Approximately 1×10⁶ isolated heterophils were suspended in one mL RPMI 1640 including 10% fetal bovine serum (FBS) and 1% chicken serum and incubated together with one mg/ml FITC-dextran and 100nM lipopolysaccharide (LPS) at 37 °C with 5% CO₂ for two hours. Incubation was stopped by mixing with ice-cold phosphate buffered solution (PBS), after which the cells were washed two times with cold PBS and analyzed using a FACS Calibur flow cytometer (FACS Calibur, BD Bioscience, USA). The endocytic activity was expressed as a positive rate of cells with FITC⁺.

Production of oxide during the oxidative burst was determined by oxidation of intracellular dihydrorhodamine 123 (DHR123) to green, fluorescent rhodamine 123 [46]. Approximately 1×10⁶ isolated heterophils were suspended in one mL RPMI 1640 including 10% FBS and 1% chicken serum. The heterophils were stimulated by one μg/ml phorbol myristate acetate (PMA), after which one μg/ml DHR123 was added to the medium and co-incubated at 37°C with 5% CO₂ for two hours. Heterophils were washed twice with cold PBS, resuspended in PBS, and analyzed with a BD FACS flow cytometer (BD Bioscience, USA). Oxidative burst activity was reported as a positive rate of cells with DHR123⁺. Flow cytometry data were analyzed with FlowJo software version 10 (TreeStar, Ashland, OR) and results per group were compared using Student’s t-test. All data were shown as mean ± standard error mean (SEM).

Whole genome sequencing and indel detection

Whole-genome resequencing of 368 individuals comprising G5 of selection line (G5), G9 of selection line (G9), non-selection line (NS) and F₂ of G9 constructed by G9 of the selection line and non-selection line, was performed. Forty-three chickens obtained from G5, 92 chickens obtained from G9 and 92 chickens obtained from NS, were used in this study. A total of 141 F₂ chickens were also used in the current study, obtained from the cross of G9 of selection line and non-selection line. In cross design, the F₁ population were established by 15 roosters (G9) mated with 30 hens (NS) and the F₂ population was established by using the F₁ hybrids based on the principle of random mating. The chickens’ genomic DNAs were extracted from blood samples from wing veins using the standard phenol-chloroform extraction procedure. The quality of the DNA was examined using the NanoPhotometer® (Implen, CA, USA), the concentration of the DNA was determined by Qubit® 3.0 Flurometer (Life Technologies, CA, USA) and the integrity of the DNA was evaluated using agarose gel electrophoresis. For genome sequencing, 0.5 μg of genomic DNA from each sample was used to construct paired-end sequencing libraries using an Annoroad® Universal DNA Fragmentase kit V2.0(AN200101-L) and an Annoroad® Universal DNA Library Prep Kit V2.0 (AN200101-L) following the manufacturer’s recommendations and index codes and sequenced on the Illumina HiSeq X Ten Sequencer (Illumina Inc) at a mean depth of ×10.

Raw data were processed with Perl scripts to ensure the quality of data used in further analysis. The adopted filtering criteria were the removal of the adaptor-polluted reads where reads containing more than five adapter-polluted bases were regarded as adaptor-polluted reads and would be filtered out, removal of the low-quality reads, that had a number of low-quality bases and a phred quality value less than 19 meaning that more than 50% of total bases were regarded as low-quality reads and removal of reads with the number of N bases accounting for more than 5 %. For paired-end sequencing data, both reads were filtered out if any read of the paired-end reads was adaptor-polluted. The obtained clean data after filtering was analyzed statistically on its quantity and quality, including Q30, data quantity and base content statistics.

After sequence quality filtering, the Burrows–Wheeler aligner (BWA 0.7.12) [47] was used to map the clean reads to the chicken reference genome GRCg6a, GCF_000002315.6. Samtools v1.2 [48] were used to sort reads and duplicate reads resulting from polymerase chain reaction (PCR) were removed using Picard tools v1.13 (http://broadinstitute.github.io/picard/). Indels were filtered before the following analysis with GATK VariantFiltration protocol [49].The filtering settings were as follows: QD<2.0, ReadPosRankSum<-20.0, FS>200.0, QUAL <30.0, DP<4.0. The annotation was performed using ANNOVAR [50] for all the qualified variants based on the GFF file.

Population structure and phylogenetic tree

Principal component analysis (PCA) was conducted using the PLINK software (Version: 1.90)[23]. For the first PCA plot, all 227 chickens were used, with the first three principal components cumulatively explaining 15.41% of the total variance. Thirty individuals were selected in each population for subsequent analysis. Population structure was evaluated using ADMIXTURE (Version 1.3.0) (https://dalexander.github.io/admixture/) and assuming two ancestral populations (K). After converting the PLINK indel matrix to distance matrix, the “meg.” format obtained by the “plink.distance.matrix.to.mega.pl” script, then MEGA7.0 [51] and iTOL software [52] were used to visualize the phylogenetic trees.

Detection of selective sweeps

All indels passing the quality control were used to detect signature of selection in G9 of selection line and non-selection line within the 184 genomes included 92 individuals in each population, respectively. The Fst fixation and nucleotide diversity ratio (θπ) were calculated within 40 kb windows sliding in 20 kb steps. The θπ ratios of the two population were calculated and then log2-transformed the data. Those windows comprised within the top 5% quantile of all two statistics were considered as candidate outliers and annotated using the ANNOVAR based on the GFF file. Functional enrichment for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses was assessed with Metascape [53].

Genome-wide association study

The univariate linear mixed model (LMM) implemented in GEMMA version 0.98.1 software (https://github.com/genetics-statistics/GEMMA/releases) was used for the association analysis [54].The first three principal component values PCA eigenvectors derived from the whole-genome indels were fitted as fixed effects in the mixed model to correct for population stratification [55]. The parameter ‘--indep-pairwise 25 5 0.2’ was used to infer effective independent tests. Considering the over-conservation of 5% Bonferroni correction method, the suggestive and whole-genome significance cutoff was defined as the Bonferroni threshold, which is 1.00/ independent indels (-log10 [P] =5.01) and 0.05/ independent indels (-log10 [P] =6.31). Manhattan and quantile–quantile (Q-Q) plots were constructed by the “qqman” package (http://cran.r-project.org/package=qqman) in RStudio-1.3.1093.

H/L: Heterophil/lymphocyte ratio; INDEL: Insertion-deletion; PCA: Principal component analysis; MAF: Minor allele frequency; Fst: The population differentiation index; θπ: Nucleotide diversity; Q-Q: Quantile–quantile; GWAS: Genome-wide association study

Data availability statement

The genome sequence data reported in this article is being uploaded to the genome sequence file of the BIG Data Center of the Beijing Institute of Genomics, Chinese Academy of Sciences and is publicly available from http://bigd.big.ac.cn (CRA003684).

Funding

This research was supported by grants from National Natural Science Foundation of China (No. 32072708), National Key R&D Program of China (2018YFE0128000) and Major scientific research projects (CAAS-ZDRW202005) of Chinese Academy of Agricultural Sciences.

Authors' information

Affiliations

¹Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing100193, China.

Authors' contributions

G.Z., J.Z. and J.W. conceived the project and designed the research. J.Z., J.W., J.D. and Z.W contributed to collecting samples and measuring phenotypic data. G.Z. and J.Z wrote the manuscript. J.W., M.T., Q.L., Q.W. revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Guiping Zhao

Ethics declarations

Ethics approval and consent to participate

All experimental procedures involving Jingxing-yellow chickens were performed according to the guidelines for experimental animals established by the Ministry of Science and Technology (Beijing, China). Ethical approval on animal survival was given by the animal welfare and ethics committee of the Institute of Animal Sciences (IAS), Chinese Academy of Agricultural Sciences (CAAS, Beijing, China) with the following reference number: IASCAAS-AE20140615. This study was carried out in compliance with the ARRIVE guidelines for reporting animal research.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

Not applicable

Fulton JE. Selection for avian immune response: a commercial breeding company challenge. Poult Sci2004, 83(4):658-661. doi:10.1093/ps/83.4.658
Ken N, Evans MR. Ecological immunology: life history trade-offs and immune defense in birds. Behavioral Ecology2000(1):19-26(18).
Lee KA. Linking immune defenses and life history at the levels of the individual and the species. Integrative and comparative biology2006, 46(6):1000-1015. doi:10.1093/icb/icl049
Palacios MG, Cunnick JE, Vleck D, Vleck CM. Ontogeny of innate and adaptive immune defense components in free-living tree swallows, Tachycineta bicolor. Dev Comp Immunol2009, 33(4):456-463. doi:10.1016/j.dci.2008.09.006
Davis AK, Maerz DLMC. The use of leukocyte profiles to measure stress in vertebrates: a review for ecologists. Functional Ecology2008, 22(5):760-772.
MacColl E, Vanesky K, Buck JA, Dudek BM, Eagles-Smith CA, Heath JA, Herring G, Vennum C, Downs CJ. Correlates of immune defenses in golden eagle nestlings. Journal of experimental zoology Part A, Ecological and integrative physiology2017, 327(5):243-253. doi:10.1002/jez.2081
Krams I, Vrublevska J, Cirule D, Kivleniece I, Krama T, Rantala MJ, Sild E, Hõrak P. Heterophil/lymphocyte ratios predict the magnitude of humoral immune response to a novel antigen in great tits (Parus major). Comparative biochemistry and physiology Part A, Molecular & integrative physiology2012, 161(4):422-428. doi:10.1016/j.cbpa.2011.12.018
al-Murrani WK, Kassab A, al-Sam HZ, al-Athari AM. Heterophil/lymphocyte ratio as a selection criterion for heat resistance in domestic fowls. Br Poult Sci1997, 38(2):159-163. doi:10.1080/00071669708417962
Al-Murrani WK, Al-Rawi IK, Raof NM. Genetic resistance to Salmonella typhimurium in two lines of chickens selected as resistant and sensitive on the basis of heterophil/lymphocyte ratio. Br Poult Sci2002, 43(4):501-507. doi:10.1080/0007166022000004408
Thiam M, Barreto Sánchez AL, Zhang J, Zheng M, Wen J, Zhao G, Wang Q. Association of Heterophil/Lymphocyte Ratio with Intestinal Barrier Function and Immune Response to Salmonella enteritidis Infection in Chicken. Animals : an open access journal from MDPI2021, 11(12). doi:10.3390/ani11123498
Campo JL, Davila SG. Estimation of heritability for heterophil:lymphocyte ratio in chickens by restricted maximum likelihood. Effects of age, sex, and crossing. Poult Sci2002, 81(10):1448-1453. doi:10.1093/ps/81.10.1448
Hõrak, Peeter, Saks, Lauri, Ots, Indrek, Kollist, Hannes. Repeatability of condition indices in captive greenfinches (Carduelis chloris). Canadian Journal of Zoology2002.
Dhabhar FS, McEwen BS. Acute stress enhances while chronic stress suppresses cell-mediated immunity in vivo: a potential role for leukocyte trafficking. Brain, behavior, and immunity1997, 11(4):286-306. doi:10.1006/brbi.1997.0508
Dhabhar FS, Miller AH, McEwen BS, Spencer RL. Stress-induced changes in blood leukocyte distribution. Role of adrenal steroid hormones. Journal of immunology (Baltimore, Md : 1950)1996, 157(4):1638-1644.
Harmon BG. Avian heterophils in inflammation and disease resistance. Poult Sci1998, 77(7):972-977. doi:10.1093/ps/77.7.972
Genovese KJ, He H, Swaggerty CL, Kogut MH. The avian heterophil. Dev Comp Immunol2013, 41(3):334-340. doi:10.1016/j.dci.2013.03.021
Swaggerty CL, Kogut MH, Ferro PJ, Rothwell L, Pevzner IY, Kaiser P. Differential cytokine mRNA expression in heterophils isolated from Salmonella-resistant and -susceptible chickens. Immunology2004, 113(1):139-148. doi:10.1111/j.1365-2567.2004.01939.x
Swaggerty CL, Pevzner IY, Lowry VK, Farnell MB, Kogut MH. Functional comparison of heterophils isolated from commercial broiler chickens. Avian pathology : journal of the WVPA2003, 32(1):95-102. doi:10.1080/0307945021000070769
Chuammitri P, Ostojić J, Andreasen CB, Redmond SB, Lamont SJ, Palić D. Chicken heterophil extracellular traps (HETs): novel defense mechanism of chicken heterophils. Veterinary immunology and immunopathology2009, 129(1-2):126-131. doi:10.1016/j.vetimm.2008.12.013
Maxwell MH, Robertson GW. The avian heterophil leucocyte: A review. World's Poultry Science Journal1998, 54(2):155-178.
Swaggerty CL, Pevzner IY, Ferro PJ, Crippen TL, Kogut MH. Association between in vitro heterophil function and the feathering gene in commercial broiler chickens. Avian pathology : journal of the WVPA2003, 32(5):483-488. doi:10.1080/0307945031000154071
Redmond SB, Chuammitri P, Andreasen CB, Palić D, Lamont SJ. Genetic control of chicken heterophil function in advanced intercross lines: associations with novel and with known Salmonella resistance loci and a likely mechanism for cell death in extracellular trap production. Immunogenetics2011, 63(7):449-458. doi:10.1007/s00251-011-0523-y
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJet al. PLINK: a tool set for whole-genome association and population-based linkage analyses. American journal of human genetics2007, 81(3):559-575. doi:10.1086/519795
Swaggerty CL, He H, Genovese KJ, Kaiser P, Pevzner IY, Kogut MH. The feathering gene is linked to degranulation and oxidative burst not cytokine/chemokine mRNA expression levels or Salmonella enteritidis organ invasion in broilers. Avian pathology : journal of the WVPA2006, 35(6):465-470. doi:10.1080/03079450601028829
Lee CC, Wu CC, Lin TL. Chicken melanoma differentiation-associated gene 5 (MDA5) recognizes infectious bursal disease virus infection and triggers MDA5-related innate immunity. Archives of virology2014, 159(7):1671-1686. doi:10.1007/s00705-014-1983-9
Xie L, Xie Z, Wang S, Huang J, Deng X, Xie Z, Luo S, Zeng T, Zhang Y, Zhang M. Altered gene expression profiles of the MDA5 signaling pathway in peripheral blood lymphocytes of chickens infected with avian reovirus. Archives of virology2019, 164(10):2451-2458. doi:10.1007/s00705-019-04340-8
García-Ortiz A, Serrador JM. ERM Proteins at the Crossroad of Leukocyte Polarization, Migration and Intercellular Adhesion. International journal of molecular sciences2020, 21(4). doi:10.3390/ijms21041502
Röhrl J, Yang D, Oppenheim JJ, Hehlgans T. Specific binding and chemotactic activity of mBD4 and its functional orthologue hBD2 to CCR6-expressing cells. The Journal of biological chemistry2010, 285(10):7028-7034. doi:10.1074/jbc.M109.091090
Herlihy SE, Pilling D, Maharjan AS, Gomer RH. Dipeptidyl peptidase IV is a human and murine neutrophil chemorepellent. Journal of immunology (Baltimore, Md : 1950)2013, 190(12):6468-6477. doi:10.4049/jimmunol.1202583
Pickering MC, Warren J, Rose KL, Carlucci F, Wang Y, Walport MJ, Cook HT, Botto M. Prevention of C5 activation ameliorates spontaneous and experimental glomerulonephritis in factor H-deficient mice. Proc Natl Acad Sci U S A2006, 103(25):9649-9654. doi:10.1073/pnas.0601094103
Kogut MH, Iqbal M, He H, Philbin V, Kaiser P, Smith A. Expression and function of Toll-like receptors in chicken heterophils. Dev Comp Immunol2005, 29(9):791-807. doi:10.1016/j.dci.2005.02.002
Rodgers LS, Beam MT, Anderson JM, Fanning AS. Epithelial barrier assembly requires coordinated activity of multiple domains of the tight junction protein ZO-1. Journal of cell science2013, 126(Pt 7):1565-1575. doi:10.1242/jcs.113399
Abdala-Valencia H, Kountz TS, Marchese ME, Cook-Mills JM. VCAM-1 induces signals that stimulate ZO-1 serine phosphorylation and reduces ZO-1 localization at lung endothelial cell junctions. J Leukoc Biol2018, 104(1):215-228. doi:10.1002/jlb.2ma1117-427rr
Zhu B, Li Q, Liu R, Zheng M, Wen J, Zhao G. Genome-Wide Association Study of H/L Traits in Chicken. Animals : an open access journal from MDPI2019, 9(5). doi:10.3390/ani9050260
Zhao B, Ye X, Yu J, Li L, Li W, Li S, Yu J, Lin JD, Wang CY, Chinnaiyan AMet al. TEAD mediates YAP-dependent gene induction and growth control. Genes & development2008, 22(14):1962-1971. doi:10.1101/gad.1664408
Gao Y, Sun Y, Ercan-Sencicek AG, King JS, Akerberg BN, Ma Q, Kontaridis MI, Pu WT, Lin Z. YAP/TEAD1 Complex Is a Default Repressor of Cardiac Toll-Like Receptor Genes. International journal of molecular sciences2021, 22(13). doi:10.3390/ijms22136649
Deb I, Manhas N, Poddar R, Rajagopal S, Allan AM, Lombroso PJ, Rosenberg GA, Candelario-Jalil E, Paul S. Neuroprotective role of a brain-enriched tyrosine phosphatase, STEP, in focal cerebral ischemia. The Journal of neuroscience : the official journal of the Society for Neuroscience2013, 33(45):17814-17826. doi:10.1523/jneurosci.2346-12.2013
Poddar R, Deb I, Mukherjee S, Paul S. NR2B-NMDA receptor mediated modulation of the tyrosine phosphatase STEP regulates glutamate induced neuronal cell death. Journal of neurochemistry2010, 115(6):1350-1362. doi:10.1111/j.1471-4159.2010.07035.x
Han Z, Zhang Y, Wang C, Liu X, Jiang A, Liu Z, Wang J, Yang Z, Wei Z. Ochratoxin A-Triggered Chicken Heterophil Extracellular Traps Release through Reactive Oxygen Species Production Dependent on Activation of NADPH Oxidase, ERK, and p38 MAPK Signaling Pathways. Journal of agricultural and food chemistry2019, 67(40):11230-11235. doi:10.1021/acs.jafc.9b03155
Kogut M, Lowry VK, Farnell M. Selective pharmacological inhibitors reveal the role of Syk tyrosine kinase, phospholipase C, phosphatidylinositol-3'-kinase, and p38 mitogen-activated protein kinase in Fc receptor-mediated signaling of chicken heterophil degranulation. International immunopharmacology2002, 2(7):963-973. doi:10.1016/s1567-5769(02)00050-4
Pulli I, Lassila T, Pan G, Yan D, Olkkonen VM, Törnquist K. Oxysterol-binding protein related-proteins (ORPs) 5 and 8 regulate calcium signaling at specific cell compartments. Cell calcium2018, 72:62-69. doi:10.1016/j.ceca.2018.03.001
Booth DM, Várnai P, Joseph SK, Hajnóczky G. Oxidative bursts of single mitochondria mediate retrograde signaling toward the ER. Molecular cell2021, 81(18):3866-3876.e3862. doi:10.1016/j.molcel.2021.07.014
Winterbourn CC, Kettle AJ. Redox reactions and microbial killing in the neutrophil phagosome. Antioxid Redox Signal2013, 18(6):642-660. doi:10.1089/ars.2012.4827
Campbell TW. Evaluation of the Blood Film. Clinics in laboratory medicine2015, 35(3):703-721. doi:10.1016/j.cll.2015.05.016
Yi X, Zeng C, Liu H, Chen X, Zhang P, Yun BS, Jin G, Zhou A. Lack of RNase L attenuates macrophage functions. PLoS One2013, 8(12):e81269. doi:10.1371/journal.pone.0081269
Avendaño A, Sales-Pardo I, Marin L, Marin P, Petriz J. Oxidative burst assessment and neutrophil-platelet complexes in unlysed whole blood. Journal of immunological methods2008, 339(2):124-131. doi:10.1016/j.jim.2008.09.003
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics (Oxford, England)2009, 25(14):1754-1760. doi:10.1093/bioinformatics/btp324
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The Sequence Alignment/Map format and SAMtools. Bioinformatics (Oxford, England)2009, 25(16):2078-2079. doi:10.1093/bioinformatics/btp352
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly Met al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res2010, 20(9):1297-1303. doi:10.1101/gr.107524.110
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic acids research2010, 38(16):e164. doi:10.1093/nar/gkq603
Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Mol Biol Evol2016, 33(7):1870-1874. doi:10.1093/molbev/msw054
Letunic I, Bork P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic acids research2021, 49(W1):W293-w296. doi:10.1093/nar/gkab301
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun2019, 10(1):1523. doi:10.1038/s41467-019-09234-6
Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet2012, 44(7):821-824. doi:10.1038/ng.2310
Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet2006, 38(8):904-909. doi:10.1038/ng1847

No competing interests reported.

Additionalfile1.FigureS1.pdf
Additional file 1. Figure S1 PCA pots of the first two components of 249 chickens in GWAS. G5: G5 of selection line; G9: G9 of selection line; F₂: F₂ of G9.
Additionalfile2.TableS16.xlsx
Additional file 2. Table S1 Top5％ windows of genome-wide selective sweep; Table S2 Genes annotated in the top 5％ window; Table S3 Top 20 clusters of GO biological processes and KEGG Pathway across putative genes; Table S4 Descriptive statistics for H/L traits in three population; Table S5 Significant indels associated with H/L ratio.

Download PDF

Version 1

posted

You are reading this latest preprint version

Immunogenetic Basis of Chicken’s Heterophil/Lymphocyte Ratio Revealed by Genome-wide Indel Variants Analysis

Status:

Version 1

Abstract

Figures

Background

Results

Phagocytosis and oxidative burst function of heterophils were different between the H/L selection line and non-selection line

The population structure differentiated with the directional selection of H/L

Identification of signature of selection between H/L selection line and non-selection line

Identification of genes associated with H/L ratio

Discussion

Conclusion

Materials And Methods

Experimental birds

Phenotype determination

Detection of phagocytosis and oxidative burst function of heterophils

Whole genome sequencing and indel detection

Population structure and phylogenetic tree

Detection of selective sweeps

Genome-wide association study

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1