lmmunogenetic basis of chicken’s heterophil to lymphocyte ratio revealed by genome-wide indel variants analysis

ZHANG Jin,WANG Jie,WANG Qiao,CUl Huan-xian,DlNG Ji-qiang,WANG Zi-xuan,Mamadou THlAM,Ll Qing-he,ZHAO Gui-ping#

1 Institute of Animal Sciences,Chinese Academy of Agricultural Sciences,Beijing 100193,P.R.China

2 Poultry Institute,Shandong Academy of Agricultural Sciences,Jinan 250023,P.R.China

Abstract

Enhancing host immunity is an effective way to reduce morbidity in chickens.Heterophil to lymphocyte ratio (H/L) is associated with host disease resistance in birds.Chickens with different H/L levels show different disease resistances.However,the utility of the H/L as an indicator of immune function needs to be further analyzed.In this study,a H/L directional breeding chicken line (Jingxing yellow chicken) was constructed,which has been bred for 12 generations.We compared the function of heterophils,and combined statistical analysis to explore the candidate genes and pathways related to H/L.The oxidative burst function of the heterophils isolated from the H/L selection line (G12) was increased (P=0.044) compared to the non-selection line (NS).The 22.44 Mb genomic regions which annotated 300 protein-coding genes were selected in the genome of G9 (n=92) compared to NS (n=92) 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 genomewide significant indels containing 29 protein-coding genes were also identified as associated with the H/L based on genome-wide association study (GWAS).The expression of protein tyrosine phosphatase non-receptor type 5 (PTPN5) (r=0.75,P=0.033) and oxysterol binding protein like 5 (OSBPL5) (r=0.89,P=0.0027) were positively correlated with H/L.Compared to the high H/L group,the expressions of PTPN5 and OSBPL5 were decreased (P＜0.05) in the low H/L group of Beijing you chicken.The A/A allelic frequency of indel 5_13108985 (P=3.85E–06) within OSBPL5 gradually increased from the NS to G5 and G9,and the individuals with A/A exhibited lower H/L than individuals with heterozygote A/ATCT (P=4.28E–04) and homozygous ATCT/ATCT (P=3.40E–05).Above results indicated oxidative burst function of heterophils were enhanced,and 22.44 Mb genomic regions were selected with the directional selection of H/L.In addition,PTPN5 and OSBPL5 genes were identified as H/L-related candidate genes.These findings revealed the complex genetic mechanism of H/L related to immunity and will allow selection for improving chicken immunity based on the H/L.

Keywords: H/L,heterophil,indel,selective sweep,GWAS

1.lntroduction

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 (Fulton 2004).

Count of relative leucocytes is the simplest but most robust measure to assess animal immune function using blood smears (Ken and Evans 2000),the heterophils and lymphocytes are the two most numerous immune cell types in avian leucocyte profiles (Lee 2006; Palaciosetal.2009).The heterophil to lymphocyte ratio (H/L) has been shown to be correlated with the strength of innate immune response in poultry (Davisetal.2008; Kramsetal.2012; MacColletal.2017).Al-Murranietal.(1997,2002) used the H/L as a selection criterion for heat stress resistance,and also showed that H/L was a criterion for selection of resistance toSalmonellatyphimuriumin chickens.Thiametal.(2021) reported that the H/L was associated with the intestinal barrier and immune response forSalmonellaenteritidisclearance and the chickens with low H/L had enhanced intestinal immunity.The high heritability of the H/L supported its use as an indicator for disease resistance in poultry breeding (Campo and Davila 2002).

The baseline H/L of individuals has been shown to be highly stable long-term (Hõraketal.2002).The rapid increase of H/L in blood is mainly affected by the increase in heterophils,for example,tissue injury due to infection can cause a short-term change of H/L,as lymphocytes migrate from the circulating blood to the skin,spleen and lymph nodes,while heterophils migrate from bone marrow into the blood (Dhabharetal.1996; Dhabhar and McEwen 1997).Heterophils play an indispensable role in the innate cellular immunity of the avian host and their biology resembles that of mammalian neutrophils (Harmon 1998).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 (Genoveseetal.2013).On detection of pathogens,heterophils respond and relay information to other immune cells through a network of intracellular signaling pathways and the release and response to cytokines and chemokines (Swaggertyetal.2004).Phagocytosis is the main process used by heterophils to kill invading pathogens,immediately followed by degranulation and production of an oxidative burst (Maxwell and Robertson 1998).Heterophil function is under immune control and it is possible to genetically select for the enhancement of heterophil function (Swaggertyetal.2003a,b; Redmondetal.2011).

However,the immunogenetic basis of H/L has not been fully elucidated to date.Therefore,we hypothesized that the function of heterophil might be enhanced with the directional selection with H/L,and selection signatures and key genes associated with H/L might be accurately identified through constructed H/L directional selection population.To determine the hypothesis,this study used a line of Jingxing yellow chickens with selection on H/L for 12 generations.We checked the function difference of heterophil between selection line and non-selection line.We also performed whole-genome resequencing of the H/L selection line and non-selection line to reveal selection signatures and H/L-associated genes.The H/L selection line might be able to explain the genetic mechanism that contributes to the more resistant state,so identifying pathways or genes associated with the H/L in poultry provides the opportunity for its wider application in genetic modification and selective breeding.

2.Materials and methods

2.1.Animals,experimental design and treatments

Jingxing yellow chickens,raised in Changping Base of Chinese Academy of Agricultural Sciences,were used for current study.The chickens were selected on H/L to form a H/L directional selection line,with the base population consisting of 200 males and 500 females.All individuals were measured H/L phenotype of peripheral blood at 56 days old.Individuals with low H/L were selected to build 30 families in each generation,with a male:female ratio of 1:3 (30 roosters and 90 hens were used to breed the next generation).In order to avoid inbreeding,roosters and hens in the same family were not full siblings and hens within a family were not half siblings in each generation.At the same time,a control population (nonselection line) was breed using pooled semen in each next generation.In cross design,the F1population was established by 15 roosters (G9) mated with 30 hens (NS) and the F2population was established by using the F1hybrids based on the principle of random mating.All these chickens were reared in semi-open chicken house and fed with traditional corn–soybean meal diet according to the Feeding Standard of Chicken established by the Ministry of Agriculture,Beijing,China (GB/T 7714-2015 2004).The feed formula was shown in Appendix A.The average body weight was about 1 200 g at 56 days old based on statistics for three consecutive years.Blood samples of 368 Jingxing yellow chickens (56 days old) were collected from wing vein and stored in 5 mL EDTA sodium anticoagulant tube at –20°C.

2.2.H/L determination

The H/L 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 1 000× and the H/L was calculated by dividing the number of heterophils by the number of lymphocytes (Campbell 2015).

2.3.Detection of phagocytosis and oxidative burst function of heterophils

Heterophils were isolated from the peripheral blood of chicks at 7 days post-hatch using the commerciallyavailable Chicken Heterophils Isolation Kit (Solar Technology Co.,China) but with improved.Briefly,disodium ethylenediaminetetraacetic acid (EDTA)-anticoagulated blood was diluted 1:1 with RPMI-1640 media containing 1% methylcellulose (25 centipopes; Sigma Chemical Co.,St.Louis,MO,USA) and centrifuged at 50×g for 30 min.The buffy coat layers were retained and suspended in Ca2+and Mg2+-free Hank’s balanced salt solution (HBSS,1:1),layered over a discontinuous gradient reagent A and reagent C in the kit,and centrifuged at 800×g for 30 min.After centrifugation,the heterophil layers were collected and washed twice in RPMI-1640 (1:1).The cells were then re-suspended in fresh RPMI-1640,counted on a haemocytometer.

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 (Yietal.2013).Approximately 1×106isolated heterophils were suspended in 1 mL RPMI 1640 including 10% fetal bovine serum (FBS) and 1% chicken serum and incubated together with 1 mg mL–1FITC-dextran and 100 nmol L–1lipopolysaccharide (LPS) at 37°C with 5% CO2for 2 h.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 (Avendañoetal.2008).Approximately 1×106isolated heterophils were suspended in 1 mL RPMI 1640 including 10% FBS and 1% chicken serum.The heterophils were stimulated by 1 μg mL–1phorbol myristate acetate (PMA),after which 1 μg mL–1DHR123 was added to the medium and co-incubated at 37°C with 5% CO2for 2 h.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’st-test.All data were shown as mean±standard error.

2.4.Whole-genome sequencing and indel detection

Whole-genome resequencing of 368 individuals comprising G5 of selection line (G5,n=43),G9 of selection line (G9,n=92),non-selection line (NS,n=92) and F2(n=141) constructed by G9 and NS.The chickens’ genomic DNAs were extracted 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 Annoroad®Universal DNA Fragmentase Kit V2.0 (AN200101-L) and 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.,USA) at a mean depth of 10×.The genome sequence data were uploaded to the genome sequence file of the BIG Data Center of the Beijing Institute of Genomics,Chinese Academy of Sciences and are publicly available from http://bigd.big.ac.cn (CRA003684).

Raw data were processed with Perl scripts to ensure the quality of data used in further analysis.The filtering criteria were as follows: (1) remove reads containing adapter sequence (more than five bases); (2) remove low-quality reads (more than 50% of bases having Phred Quality value less than 19); and (3) remove reads in which more than 5% of bases are N.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) (Li and Durbin 2009) was used to map the clean reads to the chicken reference genome (ftp://ftp.ensembl.org/pub/release-96/fasta/gallus_gallus/dna/) Samtools v1.2 (Lietal.2009) were used to sort reads and duplicate reads resulting from PCR were removed using Picard tools v1.13 (http://broadinstitute.github.io/picard/).Indels were filtered before the following analysis with GATK Variant Filtration protocol (McKennaetal.2010).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 (Wangetal.2010) for all the qualified variants based on the GFF file.

2.5.Population structure and phylogenetic tree

Principal component analysis (PCA) was conducted using the PLINK Software (Version 1.90) (Purcelletal.2007).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 (Kumaretal.2016) and iTOL Software (Letunic and Bork 2021) were used to visualize the phylogenetic trees.

2.6.Detection of selective sweeps

All indels passing the quality control were used to detect signature of selection in G9 of selection line and nonselection 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.

2.7.Enrichment analysis of genes under potential selection in H/L selection line and non-selection line

To investigate the function of selected genes between selection line and non-selection line by selective sweep analysis,the Metascape (Zhouetal.2019) was used to perform functional enrichment for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses.Thresholds were set toP＜0.01,a minimum count of three and an enrichment factor exceeding 1.5.

2.8.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 (Zhou and Stephens 2012).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 (Priceetal.2006).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.

2.9.Transcriptome sequencing and analysis

Total RNA was isolated from livers of Beijing you chicken (n=8).The purity and construction of RNA-seq library refer to the published method (Barreto Sánchezetal.2022).In brief,genomic DNA was removed by using the TIANGEN DNase Kit (Tiangen Biochemical Technology (Beijing) Co.,Ltd.,China).

RNA was purified using the kaiaoK5500®Spectrophotometer (Kaiao,Beijing,China),and the integrity and concentration of RNA were determined using the Bioanalyzer 2100 RNA Nano 6000 Assay Kit (Agilent Technologies,CA,USA).mRNA was purified using poly-T oligo-attached magnetic beads,and target products were identified.The purified double-stranded cDNA was repaired at the end,addition of base A and sequencing junction,and recovered by agarose gel electrophoresis for target size fragment and PCR amplification,and the library was completed.Finally,the quality control of the sequencing data was performed using FastQC (version 0.11.5).The mRNA sequencing data were deposited in the Genome Sequence Archive repository (https://ngdc.cncb.ac.cn/gsa/),accession number CRA007168.

2.10.Statistical analysis

Student’st-test was used for analysis of significant differences between groups.P＜0.05 implied a statistically significant difference.Each replication was an experimental unit.Data were visualized using GraphPad Prism version 8 (GraphPad Software,San Diego,CA,USA) and R version 4.0.5,and expressed as the mean±standard deviation (SD) unless otherwise indicated.

3.Results

3.1.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 line (G12) and nonselection line were compared at 7 days post-hatch.The phagocytosis was evaluated for the ability to phagocytize FITC-dextran by flow cytometry.Heterophils isolated from selection line showed a higher percentage of contained FITC-dextran (35.50±2.56)% compared to non-selection line (29.75±2.18)%,but this difference was not significant (P=0.521) (Fig.1-A).

Fig.1 Functional comparison of heterophils isolated from selection line (S) and non-selection line (NS).A,phagocytosis of heterophils isolated from S (n=6) and NS (n=4) lines when stimulated with the inflammatory agonist lipopolysaccharide (LPS).B,heterophils isolated from S (n=6) and NS (n=6) lines generate oxidative burst response when stimulated with the inflammatory agonist PMA.Data were presented as the mean±standard error.Scatter plot showed the data for one sample in each group.

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 (Fig.1-B).Heterophils isolated from the selection line (45.50±0.94)% generated significantly higher oxidative burst level than non-selection line (31.90±1.90)% (P=0.044).The result indicated the function of the heterophils could be enhanced though H/L directional breeding.

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

There were 1 068 153 insertion/deletion (indels) polymorphisms identified on 1 to 28 autosomes after filtering using PLINK (version 1.9) Software (Purcelletal.2007) 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 nonselection line clustered in separate populations (Fig.2-A) 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 (Fig.2-B) 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 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 (Fig.2-C).

Fig.2 Population genetics analyses of directed selection line and non-selection line.A,principal component analysis (PCA) pots of the first two components of 227 chickens.The fraction of the total variance explained was reported on each individual axis between parentheses.B,neighbor-joining tree constructed using P-distances between individuals.C,genetic structure of different population.The length of each colored segment represents the proportion of the individual’s genome from ancestral populations (K=2).G5,5th generation (G5) of selection line (n=43); G9,9th generation (G9) of selection line (n=92); NS,non-selection line (n=92).

3.3.ldentification of signature of selection between H/L selection line and non-selection line

Positive selection using the genomes of G9 (n=92) of selection line and non-selection line (n=92) was 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 index (Fst) 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% 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 Fig.3.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 Fig.4 and Appendix B.

Fig.3 Signatures of selection in selection line and nonselection line.The nucleotide diversity (θπ) ratio (x-axis) and the population genetic differentiation Fst values (y-axis) were calculated within 40-kb sliding windows (step size=20 kb).The significance threshold of selection signature was arbitrarily set to top 5% percentile outliers for each test.The top 5% quantiles of both statistics were shown in the top-right corner defined by grey dashed lines.

Fig.4 Top 20 clusters of GO biological processes and KEGG pathway across putative genes.P＜0.01,a minimum count of three and an enrichment factor exceeding 1.5.

3.4.ldentification of loci and genes associated with H/L based on GWAS

For genome-wide association study (GWAS) of the H/L,249 samples and 983 198 indels were kept after filtering using the same parameters as genetic selection analysis.Descriptive statistics of the phenotypic measurements of the H/L traits used for GWAS are provided in Appendix C.This analysis accounted for the population structure by adjustment of PCA covariates and genetic relatedness among all pairs of individuals as seen in Appendix D.Manhattan and quantile–quantile (Q–Q) plots are shown in Fig.5-A 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=7.06E–08; 8_16086712,P=1.64E–07; 5_8015233,P=2.16E–07).All genome-wide significant indels were annotated to 29 protein-coding genes and the detailed information was shown in Appendix E.The most significant indel 5_8074623 (P=7.06E–08) was an intron 1–2 of TEA domain transcription factor 1 (TEAD1) and accounted for 11.08% genetic variance for H/L.

Fig.5 Genome-wide association analysis of heterophil to ly(H/L) of 249 (G5,n=35; G9,n=73; F2,n=141) chickens.A,m for genome-wide association study (GWAS) for H/L.The and blue lines indicate the 5% genome-wide significance and suggestive thresholds (P=9.71E–06),respectively.quantile (Q–Q) plots of the GWAS for H/L.mphocyte ratio anhattan plot horizontal red (P=4.86E–07) B,quantile–

3.5.Liver transcriptome of Beijing you chicken revealed PTPN5 and OSBPL5 genes associated with H/L

To further identify candidate genes affecting H/L in other breeds of chicken,Beijing you chicken were selected as the verification population.The correlation between GWAS-related genes expression of the liver and H/L in eight 28-day-old of Beijing you chickens was compared.The Pearson correlation coefficient heatmap showed that the expression ofPCDH17was negatively correlated with H/L (r=–0.65,P=0.082),while the expression ofPTPN5(r=0.75,P=0.033) andOSBPL5(r=0.89,P=0.0027) were positively correlated with H/L as seen in Fig.6-A.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) (P=0.014) as seen in Fig.6-B.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 ofPCDH17,PTPN5andOSBPL5in high H/L and low H/L groups were compared,and the expression ofPTPN5andOSBPL5were decreased (P＜0.05) in the low H/L group than high H/L group in Beijing you chicken (Fig.6-C).

Fig.6 Screening of heterophil to lymphocyte ratio (H/L)-related candidate genes based on transcriptome.A,heat map of Pearson correlation coefficient between gene expression and H/L based on FPKM value (n=8).B,comparison of H/L phenotypic difference between H/L-high (n=3) and H/L-low groups (n=3).C,heat map of FPKM value of candidate genes in H/L-high (n=3) and H/L-low groups (n=3).D,comparison of gene expression difference between H/L-high (n=3) and H/L-low groups (n=3).

Interestingly,there were four indels near these two genes.5_12975134 (P=6.48E–06) was located 9.48 kb downstream ofPTPN5,and other three indels were located in intron1–2 ofOSBPL5,including 5_13108985 (P=3.85E–06),5_13151117 (P=4.95E–06) and 5_13170170 (P=7.23E–06).Allelic frequency analysis in the NS,G5 and G9 showed the 5_13108985 variant allelic frequency gradually increased or decreased with the selection of H/L in different generations (Fig.7-A).Jingxing yellow chickens with homozygous A/A for 5_13108985 variant showed significantly lower H/L phenotype than individuals with heterozygote A/ATCT (P=4.28E–04) and homozygous ATCT/ATCT (P=3.40E–05),respectively (Fig.7-B).

Fig.7 Allelic frequency analysis of candidate variants.A,changes in allele frequency of candidate variant in non-selection line (NS),G5 of selection line (G5) and G9 of selection line (G9).B,plot for H/L and genotype at 5_13108985 variant.

4.Discussion

The results supported our hypothesis that the function of heterophil might be enhanced with the directional selection with H/L,and selection signatures and key genes associated with H/L might be accurately identified through constructed H/L directional selection population.Cell biology assay proved directional selection of H/L enhanced the oxidative burst (P=0.044) function of heterophils in chicken.Numerous candidate genes,widely involved in differential cell molecular processes,were selected in the H/L selection line.In addition,H/Lrelated genes were identified and proven to be effective in Beijing you chicken.The findings contributed to our understanding of the genetic basis of selecting on H/L and its application in the process of disease resistance breeding in chicken.

We found the oxidative burst of heterophils isolated from selection line (G12) performed more powerfully compared to the non-selection line (P=0.044).It was inferred that this was due to genetic alterations.Previous studies reported the gene expression and function of the heterophil were under genetic control.Redmond and colleagues proved heterophils from commercially selected and non-selected genetic lines expressed cytokines differentlyinvitroexposure toSalmonellaenteritidis,where functional phenotypes of phagocytosis and bacterial killing varied depending on the breed or cross of chicken and theSLC11A1gene was associated with phagocytosis ofSalmonellaenteritidis(Redmondetal.2011).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 (Swaggertyetal.2006).Therefore,we deduced that heterophil function related genes were selected under the H/L 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 were widely involved in differential cell molecular processes,which implied complex genetic mechanisms of immune evolution in the process of H/L 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 asIFIH1orMDA5engaged in innate immune response and response to virus.In chickens,MDA5recognized 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 (Leeetal.2014).TheMDA5signaling 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 (Xieetal.2019).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 (García-Ortiz and Serrador 2020).

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),CCR6encoded chemokine (C–C motif) receptor 6 and enabled C–C chemokine binding activity and C–C chemokine receptor activity (Röhrletal.2010).Dipeptidyl peptidase 4 (DPP4) engaged in the negative regulation of neutrophil chemotaxis (Herlihyetal.2013).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 (Pickeringetal.2006).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 whichTLRs1/6/10,TLR2(types 1 and 2),TLR3,TLR4,TLR5,andTLR7have been identified in the chicken heterophil (Kogutetal.2005),andTLR5showed limited genetic diversity in G9,but not in non-selection line,suggesting positive selection for H/Lrelated 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 (Rodgersetal.2013).Abdala-Valenciaetal.(2018) demonstrated that overexpression ofTJP1blocked leukocytetrans-endothelial migration.

GWAS was performed to identify the candidate genes of the H/L.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 and played a role in its cellular functions (Zhuetal.2019).In this study,a total of 45 genome-wide significant indels associated with H/L were identified and annotated to 29 proteincoding genes.The most significant indel 5_8074623 (P=7.06E–08) was in intron1–2 of TEA domain family member 1 (TEAD1) involved in positive regulation of cell growth (Zhaoetal.2008).In immunity,TEAD1directly bounded genomic regions adjacent toTLR1,TLR2,TLR3,TLR4,TLR5,TLR6,TLR7,withTLR9andYAP/TEAD1complex acting as a default repressor of cardiac toll-like receptor genes (Gaoetal.2021).The activation of the TLR induced signaling pathway played a significant role in the antibacterial mechanism of the avian heterophil (Kogutetal.2005).The expression level of 29 candidate genes in the liver of Beijing you chicken was verified according to grouping of high and low H/L,only the expression levels ofPTPN5(r=0.75,P=0.033) andOSBPL5(r=0.89,P=0.0027) were positively correlated with H/L.Moreover,the expression ofPTPN5andOSBPL5were decreased (P＜0.05) in the low H/L group of Beijing you chicken.PTPN5could enable protein tyrosine phosphatase activity,involved in negative regulation of MAP kinase activity (Debetal.2013) and p38 MAPK was reported as a substrate ofPTPN5(Poddaretal.2010; Debetal.2013).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 (Hanetal.2019).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 (Kogutetal.2002).The above evidence suggested thatPTPN5might play a significant role in the functional regulation of avian heterophil.Another gene,OSBPL5(also known asORP5) was predicted to enable several functions,including cholesterol binding activity,phosphatidylserine transfer activity and phospholipid binding activity based on NCBI database.TheORP5/8regulated Ca2+signaling at specific membrane contact site foci appear to play roles in processes such as mitochondrial respiration and cell proliferation (Pullietal.2018),where Ca2+channels directly sense oxidative bursts and respond with Ca2+transients adjacent to active mitochondria (Boothetal.2021).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,Ca2+,calmodulin and protein kinase following stimulation (Winterbourn and Kettle 2013).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=3.85E–06),within intron1–2 ofOSBPL5had a lower H/L than individuals with heterozygote A/ATCT (P=4.28E–04) and homozygous ATCT/ATCT (P=3.40E–05),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.

5.Conclusion

The directional selection of H/L enhanced the oxidative burst function 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.ThePTPN5andOSBPL5genes were identified as candidate genes for the H/L,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 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.

Acknowledgements

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

Declaration of competing interest

The authors declare that they have no conflict of interest.

Ethical approval

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.

Appendicesassociated with this paper are available on https://doi.org/10.1016/j.jia.2022.12.012