RNA‐Seq‐based discovery of genetic variants and allele‐specific expression of two layer lines and broiler chicken

Abstract Recent advances in the selective breeding of broilers and layers have made poultry production one of the fastest‐growing industries. In this study, a transcriptome variant calling approach from RNA‐seq data was used to determine population diversity between broilers and layers. In total, 200 individuals were analyzed from three different chicken populations (Lohmann Brown (LB), n = 90), Lohmann Selected Leghorn (LSL, n = 89), and Broiler (BR, n = 21). The raw RNA‐sequencing reads were pre‐processed, quality control checked, mapped to the reference genome, and made compatible with Genome Analysis ToolKit for variant detection. Subsequently, pairwise fixation index (F ST) analysis was performed between broilers and layers. Numerous candidate genes were identified, that were associated with growth, development, metabolism, immunity, and other economically significant traits. Finally, allele‐specific expression (ASE) analysis was performed in the gut mucosa of LB and LSL strains at 10, 16, 24, 30, and 60 weeks of age. At different ages, the two‐layer strains showed significantly different allele‐specific expressions in the gut mucosa, and changes in allelic imbalance were observed across the entire lifespan. Most ASE genes are involved in energy metabolism, including sirtuin signaling pathways, oxidative phosphorylation, and mitochondrial dysfunction. A high number of ASE genes were found during the peak of laying, which were particularly enriched in cholesterol biosynthesis. These findings indicate that genetic architecture as well as biological processes driving particular demands relate to metabolic and nutritional requirements during the laying period shape allelic heterogeneity. These processes are considerably affected by breeding and management, whereby elucidating allele‐specific gene regulation is an essential step towards deciphering the genotype to phenotype map or functional diversity between the chicken populations. Additionally, we observed that several genes showing significant allelic imbalance also colocalized with the top 1% of genes identified by the FST approach, suggesting a fixation of genes in cis‐regulatory elements.


| INTRODUC TI ON
Poultry farming is one of the fastest-growing industries, and it is expected to continue growing as human population growth drives demand for meat and eggs (Mottet & Tempio, 2017;Talebi et al., 2020).
Traditionally, the domestic chicken populations were bred for two main reasons egg and meat production (Rubin et al., 2010). Over 120 million tons of meat and 1.2 trillion eggs are consumed globally each year from commercial broiler and layer suppliers. For this reason, the poultry industry keeps more than 70 billion chickens yearly to fulfill consumer demand (Qanbari et al., 2019). The selection progress in layer chicken was mainly driven by commercial egg production and efficient feed conversion, whereas in broilers, the emphasis was on rapid growth, that is, achieving a 50-60 times increase in body weight (BW) from hatch to market. This selective breeding has resulted in considerable advancements specifically for individual traits (Havenstein et al., 2003).
In comparison with other commercial broiler strains around the globe, the Cobb 500 broiler has a high growth performance and comparable breast meat production at different processing ages (Coneglian et al., 2010;Taschetto et al., 2012). The parents and grandparents of the Cobb 500 broiler belong to Cobb-Vantress Inc., which was established in 1916 and is one of the oldest broiler breeding companies in the world (www.cobbvantr ess.com). Over the past 50 years, LSL and LB layers have established themselves as the world's leading commercial laying lines. These lines have been specifically selected for their performance in egg production and carefully selected for their combination of capabilities (http://www.ltz.de/ en/layer s/). Interestingly, both lines show similar egg production performance, but genetic differences between the two lines are responsible for the physiological differences, which in turn has led to variations in certain traits such as body weight, immunity, bone metabolism, phytate degradation, and transcription abundance of genes that could be attributed to immune system processes and phosphorus metabolism (Habig et al., 2012;Iqbal et al., 2022;Reyer et al., 2021;Schmucker et al., 2021;Sommerfeld, Omotoso, et al., 2020).
Many studies suggested that the difference in growth performance between broilers and layers was mainly due to the higher feed intake in broilers (Hocking et al., 1997;Mahagna & Nir, 1996;Masic et al., 1974). Nevertheless, several recent studies compared layers to broilers and highlighted numerous breed-specific features, including growth-rate performance, physiological traits, metabolism, and immunological functions (Koenen et al., 2002;Nihashi et al., 2019;Qanbari et al., 2019;Talebi et al., 2020). One study revealed that broilers at 6 weeks of age have a 5-fold greater BW than layers, with most of the difference in BW attributed to increased growth in broilers during the first 2 weeks after hatch (Zhao et al., 2004). Scheuermann et al. (2004) found that the higher BW of broilers was primarily due to the size of skeletal muscles, as broiler muscles have more myofibers with larger diameters than layers. Comparison of liver function in broilers and layer chickens at various stages of embryonic development and hatch indicates that on day 18 of embryogenesis, broiler liver triacylglycerol level is 1.3 times higher in broilers than in layer, and 2.2 fold higher at chick hatch (Buzała et al., 2015;Sato et al., 2006).
Comparison of inflammatory responses to lipopolysaccharide (LPS) revealed that both breeds showed differences in cytokine expression and immune responses with the broiler chicken suppressing the inflammatory responses, indicating a positive correlation with growth rate (Leshchinsky & Klasing, 2001). The egg production performance of the two-layer strains LB and LSL is almost comparable, but they differ significantly in phenotypic traits, immunological function, gene expression, and metabolic activities. A multi-omics study discovered that the LB and LSL layer lines use divergent intrinsic mechanisms that shape their immune and metabolic function (Iqbal et al., 2022).
Another pan-omics longitudinal integration study at different production periods from pre-laying to the onset of laying indicated that they differed in their ability to activate immune and metabolic mechanisms, as well as unique gut-microbiota interactions (Ponsuksili et al., 2022).
Recently several studies demonstrated that computational approaches to identify variants and accurate mapping of RNA-seq reads is an effective and cost-efficient source for the detection of genomic variations (Adetunji et al., 2019;Cornwell et al., 2018;Jehl et al., 2021;Piskol et al., 2013;Tang et al., 2014). Previous studies have shown that nonsynonymous, synonymous, and noncoding SNPs derived from RNA-seq data can represent markers of genetic differentiation, particularly when they are targets for selection or involved in the regulation of gene expression (Page & Chapman, 2021;Sun et al., 2019). Allele-specific expression (ASE) is a phenomenon involving unbalanced expression results in a selection process in which one allele is preferentially expressed over another, with the potential functional consequence on phenotype (Pierre et al., 2022).
In this study, we used transcriptomic data to identify SNPs and short indels in expressed genes to investigate genetic variation in three different chicken populations. F ST screening was used to identify regions of genetic differentiation resulting from the domestication of chickens and subsequent specialization into broiler and layer lines. Focusing on SNPs detected in expressed regions, characterized variants affecting protein functions to study cis-regulated genes by analyzing allele-specific expression in the jejunum mucosa of LB and LSL strains at different time points (10,16,24,30, and 60 weeks).

| RNA-Seq data collection and preprocessing
This study used previously published RNA sequencing data that were deposited by us in the ArrayExpress database at EMBL-EBI under accession numbers: E-MTAB-9137 and E-MTAB-9109 for Lohmann Brown (LB) and Lohmann Selected Leghorn (LSL) laying hens. In

K E Y W O R D S
allele-specific expression, broiler, fixation index, laying hen, RNA-seq, SNPs addition, RNA sequencing data from the broiler (BR, Cobb500; ArrayExpress database at EMBL-EBI under accession numbers: E-MTAB-6169 and E-MTAB-12147), were used. Information on sampling and experimental design for the two-layer strains is documented in Omotoso et al. (2021) and Reyer et al. (2021).
Briefly, the broiler chicken population comprised mixed-sex Cobb 500 birds from two experiments (trial 1, n = 13; trial 2, n = 8) described previously in Metzler-Zebeli et al. (2018) and Reyer et al. (2018). RNA extraction and RNA-seq were performed from the breast, duodenum, ileum, and jejunum in trial 1 and from muscle and liver samples in trial 2. RNA-seq reads from tissues were combined and analyzed for variant discovery. A total of 200 chickens were selected from two-layer lines (Lohmann Brown [LB], n = 90; Lohmann Selected Leghorn [LSL], n = 89; Broiler [BR], n = 21). In the initial preprocessing step, all raw sequencing reads from LB, LSL, and BR were examined for quality control and trimming approach using FastQC (version 0.11.7). The information on the average number of reads/library or strains that were included in the analysis is shown in Table S1.

| Mapping and variant detection
Sequenced reads were processed in accordance with GATK (version 4.2.0.0) best practices for reliable SNP detection and genotype calling in RNAseq (Jehl et al., 2021). Initially, we used STAR (version 2.7.8a; Dobin et al., 2013) in 2-pass mode for mapping to the Ensembl Gallus gallus reference genome assembly (GRCg6a, with annotation version 104). In resulting BAM files, multiple alignments are rejected using samtools (version 1.12; Danecek et al., 2021). Further data cleanup steps included removing duplicate reads, splitting reads by cigar string to detect splicing events, and base recalibration to known variants from Ensembl v94's dbSNP (Hunt et al., 2018).
Subsequently, discovered SNP and indel variants were counted for each individual using GATK HaplotypeCaller using options for minimum confidence threshold = 20 and to avoid soft-clipped bases.
Finally, all individual gVCF files were combined and genotyped with a minimum phred-scaled confidence threshold = 20 (McKenna et al., 2010;Van der Auwera et al., 2013). The created VCF file was subjected to SNP variants selection and filtration with filters FisherStrand >30 and QualByDepth <2 to obtain SNPs for each individual. The numeric recording of genotypes was performed using the "Recode Genotypes" method within the JMP Genomics v.10.1 (SAS Institute).
Two individuals from the LB strain were excluded from the analysis after GATK SNP identification. One animal had a low SNP density, while the other animal did not cluster with other individuals from LB strain as confirmed with Bioconductor R package arrayQualityMetrics (version 3.52.0).

| Fixation index (F ST ) estimation and candidate gene selection
Initially, the fixation index (F ST ) for each SNP was calculated between the three different chicken breed groups (LB vs. LSL, LB vs. BR, and LSL vs. BR) by using Weir and Cockerham-based approach within the SNPRelate R package (version 1.24.0; Weir & Cockerham, 1984;Zheng et al., 2012). During the F ST calculation, SNPs on nonautosomes or monomorphic were excluded. Only autosomal SNPs with minor allele frequency (MAF ≥ 0.05) were shortlisted. Subsequently, the "aggregate" function within the R program-

| KEGG and IPA pathway analyses
Only the genes having the top 1% of Z (F ST ) values from all the three comparisons: (LB vs. LSL, LB vs. BR, and LSL vs. BR) were subjected to ClueGO (version 2.5.1) and Cluepedia (version 1.5.7) plugin in Cytoscape (version.3.6.1) environment (Bindea et al., 2009(Bindea et al., , 2013Shannon et al., 2003). In ClueGO, each comparison was considered a cluster, and right-sided hypergeometric tests were performed to calculate the p-value. Benjamin-Hochberg was used for multiple testing corrections and the genome of Gallus gallus was selected as a reference. KEGG pathways with p < 0.05 were considered significant. For ASE genes for each group of laying strain and age were subjected to pathway analysis using the Ingenuity Pathway Analysis (IPA) software (Qiagen). IPA categorizes genes based on annotated gene functions and statistically tests for the representation of functional terms within the gene list and then calculates adjusted p-values using the Benjamini-Hochberg approach.

| Variant effect predictor (VEP) analysis
The Ensemble variant effect prediction (VEP) tool (version 106) was used to annotate 282,587 SNPs for the species Gallus gallus. Siftb was used to determine both terms (tolerated and deleterious) and also to calculate the SIFT score (0 to 1) based on sequence homology and amino acid properties (McLaren et al., 2016). Qualitative predictions were derived from the SIFT score, such that substitutions with a score of 0.05 were termed "deleterious" and all others were considered "tolerated".

| Allele-Specific expression (ASE) analysis
To obtain allele-specific expression (ASE) for SNP variants, the GATK best practice pipeline was run a second time. Discovered SNPs from the first round (explained in Section 2.2) were Nmasked in the reference genome using bedtools (version v2.27.1) and were used for mapping to obtain unbiased STAR alignments.
Next, performing residual GATK best practice steps generated unbiased ASE counts. Only results of SNP variants that were initially masked in the reference were retained. A minor allele count of 3 or more reads, total SNP count of 30 or more, and minor allele ratio of a minimum of 1% were required per accepted SNP and sample.
Allele-specific expression (ASE)-based RNA sequencing detection was performed within LB and LSL laying strains at each age of laying hens (10, 16, 24, 30, and 60 weeks) using ASEP (Allele-Specific Expression Analysis in a Population, version 0.1.0) in the R programming environment (Fan et al., 2020). The concept of the "ASE_detection" function from the ASEP package was applied to detect significant (p-value <0.05) gene-level ASE effects within the population of each age and each strain combination. Initially, the parameters set for the ASE-detection were unphased, adaptive, and resample rate of 1e6.
Finally, ASE results for all 10 individuals per line and time point were analyzed, to identify strain-specific genes with ASE, as well as genes with ASE that were shared between weeks within strains and genes with ASE common between strains across time points by using UpSetR (version 1.4.0) (Conway et al., 2017). Statistical significance was considered for genes showing allele-specific expression with a p-value of <0.05.

| RE SULTS
An RNA sequencing-based approach was used to discover variants in three chicken breeds, LB, LSL, and BR. Initially, 200 individuals [LB (n = 90), LSL (n = 89), and BR (n = 21)] were used for variant detection analysis. With the Genome Analysis ToolKit (GATK) pipeline for variant detection, over five million (5,874,467) variants were detected from three chicken populations: LB 4,232,069,LSL 4,296,844,and BR 4,890,677. Overall, 80% of individuals shared 282,587 SNPs that were processed for pairwise F ST analysis. Additionally, a total of 292,485 commonly identified SNPs in the two-layer strains were analyzed to determine Allele-Specific Expression (ASE) within LB and LSL layer lines at each week's age (10, 16, 24, 30, and 60).

| Genetic differentiation of chicken breeds based on SNP genotype data
The breed structure was determined by using principal component analysis based on SNP genotype data. There was a strong separation between all three breeds with tight clustering within LB, LSL, and BR breeds revealing high levels of similarity among the individuals within the breed ( Figure 1). As expected, noticeable separation among the breeds was observed due to their genetic diversity.
F I G U R E 1 Principle Component Analysis (PCA) of three chicken breeds. PCA plot demonstrates the breed-specific clustering based on SNP genotype data. Individuals were colored based on their breed; the red circle represents LB, the blue circle indicates LSL, and the green circle depicts BR.

| SNPs annotation and SIFT prediction
A total of 282,587 SNPs were processed with Ensembl Variant Effect Predictor (VEP) tool, and annotation analysis of variants placed the identified SNPs in different genomic regions. Out of 282,587 variants, the majority (86,546, 30.63%) were located within predicated the coding sequence. SNPs within coding regions were classified as synonymous and nonsynonymous variants, while nonsynonymous variants were further classified as either missense or nonsense.
Of the 86,546 SNPs in the coding sequences, 70.67% were synonymous variants, 28.78% were missense variants, and only 0.55% were other variants, including stop-gained, stop-lost, start-lost, and stopretained variants as shown in Figure 2b.
Based on the fact that missense variants affect protein-coding sequences, these 24,916 variants were used for SIFT prediction analysis. According to our findings, 48.5% of SNPs were predicted to have deleterious effects (0.0-0.05), 39.9% had tolerated effects (0.05-1.0), while 6.1% and 5.5% of SNPs were classified as deleterious low confidence, and tolerated low confidence, respectively ( Figure 2c).

| Pairwise F ST analysis to detect differentiated genomic regions between LB and LSL
In the pairwise F ST analysis between LB and LSL, a total of 282,587 SNPs were processed. As part of the filtering process, 7490 A total of 8929 genes represented by more than or equal to three SNPs each were indicated to differentiate between LB and LSL. Of these, 17 genes passed the cutoff criteria of Z (F ST ) ≥ 4 and were considered candidate genes. These candidate genes were distributed on the following chromosomes: 1, 3, 5, 6, 7, 11, 12, 14, and 19. Our

| Pairwise F ST analysis to detect differentiated genomic regions between LB and BR
In the pairwise F ST analysis for the comparison between LB and BR, a total of 282,587 SNPs were processed. As part of the filtering process, 7490 nonautosomal SNPs and 40,824 monomorphic or low-frequency SNPs (MAF < 0.05) were excluded. Subsequently, 234,273 SNPs were used for gene-based aggregation of SNP-specific F ST values.
A total of 9056 genes represented by more than or equal to three SNPs were indicated to differentiate between LB and LSL. Of these, 16 genes passed the cutoff criteria of Z (F ST ) ≥ 4. These candidate genes were distributed on chromosomes: 3, 4, 5, 6, 7, 11, 15, 20, 24, 25, 27, 31, and 33, respectively. Our results revealed that the

| Pairwise F ST analysis to detect selection signature between LSL and BR
In the pairwise F ST analysis for the comparison between LSL and BR, a total of 282,587 SNPs were processed. Hence, 7490 nonautosomal SNPs and 53,654 SNPs that were monomorphic or had a MAF of <0.05 were filtered out. Subsequently, 221,443 SNPs were used for gene-based aggregation of SNP-specific F ST values.
A total of 8974 genes represented by more than or equal to three SNPs were indicated to differentiate between LB and LSL.

| Colocalized QTL regions with F ST genes
The QTL database for chickens (https://www.anima lgeno me.org/ cgi-bin/Q TLdb/ GG/brows e; chickenGRCg6a.bed) was used for mapping the F ST region. Corresponding traits for laying or broiler chickens, including egg production, egg quality, fat content, growth, meat quality, and feeding categories were considered. A total of 7052 QTL regions were used, corresponding to egg production (668), egg quality (1465) Table S2. vs. LSL were designated as clusters one, two, and three, respectively. Afterward, we used a term cluster for each comparison group.

| KEGG pathway enrichment analysis on top 1% of selected genes
Among the genes with top 1% Z (F ST ) values, cluster one had 906 genes, cluster two had 898 genes, and cluster three had 893 genes.
A total of 15 KEGG pathways (Table S3), including cell cycle, peroxisome, FoxO signaling pathway, NOD-like receptor signaling pathway, sphingolipid metabolism, necroptosis, ribosome, autophagy, endocytosis, mitophagy, glycolysis/gluconeogenesis, lysosome, inositol phosphate metabolism, phosphatidylinositol signaling system, and apoptosis were enriched by these selected genes with different proportions (Figure 7). The proportion of genes within the clusters and pathways were listed in Table 1.

| Cis-regulated genes in the jejunum mucosa of LB and LSL
Allele-specific expression analysis was performed on 292,485 identified cis-regulatory SNPs in jejunum mucosal transcripts of LB and LSL strains at five developmental stages (10, 16, 24, 30, and 60 weeks) to identify genes that are allele-specifically regulated using the gene-based ASEP approach. ASEP is based on a generalized linear mixed effects model with subject-specific random effects to account for the correlation of multiple SNPs within the same gene (Fan et al., 2020). ASEP can detect ASE across multiple individuals, in our case in each strain and aging group.
Our results showed 4945 ASE genes with a p-value of<0.01 (FDR < 5%) in at least one of the strains or aging groups (Table S4).
A number of ASE genes identified within strains and time points, as well as the strain-specific ASE genes, are listed in Table 2. A higher number of ASE genes was found in the LB strain than in the LSL strain at all time points (10, 16, 24, 30, and 60 weeks). We also found 7 ASE genes (FDR < 5%) common with fixation index analysis (Z (F ST ) ≥ 4) in the comparison between LB and LSL, including TERF, SAAL1, UBA32, F I G U R E 6 Candidate genes linked to chicken QTLs. Circosplot illustrates that genes selected as candidates in the pairwise F ST analysis were assigned to a chicken QTL database. Each comparison was designated as layer 1: LB vs. BR; layer 2: LB vs. LSL; and layer 3: LSL vs. BR. The red, green, blue, orange, aquamarine, and pink lines in each layer indicate SNPs in genes meeting the cutoff criteria of Z (F ST ) ≥ 4 (as shown in Figures 3, 4, and 5) and colocalized in QTL regions associated with egg production, egg quality, fat content, feeding, growth, and meat quality, respectively. The gray line indicates SNPs in genes that meet the cutoff criteria of Z (F ST ) ≥ 4 but are not colocalized with any of the QTLs considered here. The outer layer represents the genomic location of the QTL regions.
C1H12ORF57, KCNK17, SLC13A2, and RAC2. All ASE genes of each group with a p-value of <0.01 (Table S4) were submitted for IPA pathway analysis and we found that most ASE genes are involved in energy metabolism, including Sirtuin Signaling Pathways, Oxidative phosphorylation, and mitochondrial dysfunction (Figure 8). A high number of ASE genes was found at the 24th and 30thweeks which is the peak of laying and these ASE genes were also enriched particularly in Cholesterol Biosynthesis. While the prelaying period (10 and 16 weeks), ASE genes were enriched in Inositol Phosphate

Metabolism.
We also detected ASE genes that were strain-and age-specific ( Figure 9;  (Figure 9a). In addition, these strain-and age-specific ASEs were subjected to

| DISCUSS ION
During the last century, the success of the selective breeding of chickens specialized in laying eggs (layers) and fast-growing meat birds (broilers) became apparent. Using whole-genome sequence data, a recent study has elucidated selection forces in the genome of commercial chickens, particularly broilers and layers, compared with the ancestral red jungle fowl, and identified putative selection sweeps affecting loci likely to have been affected by the domestication and selection process associated with improved production (Qanbari et al., 2019). Our previous studies focused on the molecular and phenotypic characterization of the two laying hen strains LB and LSL. Interestingly, both lines have a similar performance on egg production but differ in body weight, immunity, bone metabolism, phytate degradation, and transcript abundance of genes that could be assigned to immune system processes and phosphorus metabolic processes (Habig et al., 2012(Habig et al., , 2014Omotoso et al., 2021;Ponsuksili et al., 2022;Sommerfeld, Omotoso, et al., 2020). Therefore, genetic variation between the two strains of laying hens (LB and LSL) and in contrast to the Cobb 500 broiler (BR) meat chickens was the subject of the present study. We used available RNA-Seq data sets of the three chicken populations to derive reliable SNPs and genotypes for analyzing the genetic differentiation between strains and ASE of genes across different developmental stages of the laying hen.

| Genes identified as candidates distinguishing between LB and LSL strains
By comparing LB vs. LSL, our pairwise fixation index analysis identified 17 genes that were selected as candidates and highly differentiated between the two-layer strains (Z (F ST ) ≥ 4). Several previous studies indicate that the immune and metabolic systems of both strains are acquired and maintained differently under high-performance conditions (Iqbal et al., , 2022Ponsuksili et al., 2022). In line with these results, our study shortlisted the CREBL2 on chromosome 1 as a candidate gene that was highly differentiated in the in this study as the ASE gene at the 24th week in LSL and is located in a QTL region for fatness and growth traits (Table S2). Consistent with a previous study that identified a QTL region on chromosome 11 associated with body weight traits (Jia et al., 2012). Furthermore, the gene SLC13A2 (Solute Carrier Family 13 Member 2), which codes for the sodium-coupled citrate transporter (Akhtar et al., 2022) is worth mentioning here. Our recent study showed that SLC13A2 is highly expressed in the intestinal mucosa, especially in the laying period (week 24), and was also identified as a biosignature associated with age and laying strain (LB vs. LSL; Ponsuksili et al., 2022). In the same study, we found that the expression of SLC13A2 in the intestinal mucosa was strongly negatively correlated with phosphorus utilization (r = 0.73; Ponsuksili et al., 2022).
Interestingly, based on the genetic variants, this study also selected SLC13A2 as a candidate gene that differs between laying lines and as an ASE gene for the LB laying strain at the 24th and 30th week; we speculate that SLC13A2 may be involved in the complex events that occur during the laying cycle and that it is strain specific during the selection process. This is consistent with a study showing that SLC13A2 and SLC35B4 are posttranscriptionally regulated in chickens via miRNAs, and the expression of SLC family genes associated with egg formation, oviposition, and embryonic development changes in chickens 3 and 20 h after ovulation (Lim et al., 2012). Interestingly, SLC13A2 and SLC35B4, which were identified as functional candidate genes in this study, are located in the QTL region that has been associated with egg number (Atzmon et al., 2008) and egg production (Goto et al., 2011), respectively. The LB and LSL strains have different immunological phenotypes, with LB having an innate immunological phenotype, whereas LSL chickens showed an adaptive immunological phenotype . Interestingly, this study also identified RAC2 on chromosome 1 as a functional candidate that is distinct between LB and LSL strains and as an ASE gene at week 16 of the LB strain. Consistent with this transcriptome analysis of the breast muscles of two native chickens with divergent feed efficiency revealed that Rac family small GTPase 2 (RAC2) is a critical signal transduction factor in immune cells and is associated with modulation of the various NADPH oxidase (NOXes) family members, responsible for the production of ROS in response to receptor activation such as inflammatory cytokines (Hordijk, 2006).
RAC2-deficient mice study shows pro-inflammatory cytokines and chemokines were inhibited (Zou et al., 2017). In chicken, ROS are produced in phagosomes of chicken heterophils by RAC2 to kill pathogens (Nambooppha et al., 2018). Together with the abundant evidence showing that RAC2 plays an immunomodulatory role, we proposed that this gene may be a factor responsible for the different immunological phenotypes of LB and LSL strains.
The role of Serum Amyloid A (SAA) family genes in the immune system of humans was well known, for instance, a variety of defense functions are induced by SAA, including neutrophil priming, leukocyte recruitment, and antiviral activity (Cai et al., 2007;Hatanaka et al., 2003;Lavie et al., 2006). Similarly, another study identified the SAAL1 gene playing a significant role in the immune defense of Oplegnathus fasciatus fish species (Revathy et al., 2012). Interestingly, our fixation index analysis shortlisted the SAAL1 as a functional candidate gene that was highly differentiated between the two strains [Z (F ST ) = 4.98] and also was identified as an ASE gene at the 24th week of the LB strain. Therefore, we suggest that further research is required to deeply understand the role of the SAAL1 gene in im-
Previously, a comparative transcriptomics study between two-layer chicken strains revealed that LSL hens had downregulated COMTD1, resulting in a reduced substrate conversion, and subsequently, LSL layers' humerus bone-breaking strengths were lower than those of LB layers (Habig et al., 2012). In the present study, we also identified the COMTD1 on chromosome 6 as a candidate gene providing additional evidence to the earlier finding.
F ST analysis revealed several genes that differ between LB and LSL, some of which show ASE. These genes play a role in metabolic and immune pathways that have been shown to be recruited differently and show a phenotypic difference between the lines against a background of consistent high egg production of both lines.

| Genes identified as candidates distinguishing between LB and BR strains
In the LB vs. BR comparison, our pairwise fixation index analysis identified 16 genes that were selected as candidates and highly differentiated between the LB and Cobb 500 broilers. According to a previous study, lipid metabolism and fat deposition processes significantly differed between broiler and layer chicken (Boschiero et al., 2018). Furthermore, during chicken embryogenesis, differences in lipid metabolism have been observed, for instance, broiler livers have higher triglyceride levels (Buzała et al., 2015). A study indicated that by activating CaM/CAMKK2/AMPK pathways, S100 Calcium Binding Protein A16 (S100A16) regulates lipid metabolism.
In another study, the S100A16 gene was shown to promote adipogenesis in the mouse 3 T3-L1 cell line (Liu et al., 2011). Interestingly, our study also identified S100A16 as a candidate gene that was differentiated between the brown layer and broiler, we assume that this gene may be responsible for higher fat accumulation in the liver of the broiler.
Together with S100A16, we shortlisted two functional candidate genes: alpha-1,2-Mannosyltransferase (ALG9/ENSGALG00000021174) is involved in lipid-linked oligosaccharide assembly and N-Acylsphingosine Amidohydrolase-1 (ASAH1) is involved in sphingolipid metabolism. These genes seem to play a vital role in lipid metabolism, but further studies in the chicken are required to explore their role in a broader spectrum. Interestingly, ALG9 is also located in the QTL region for abdominal fat weight reported previously (Gao et al., 2009). A comparative study of QTL for body weight and growth in broiler and layer found a QTL region for body weight on chromosome 4 (Podisi et al., 2013). In our study, ASAH1 was identified in this QTL region differing between layer and broiler lines.
In addition, this QTL region also covers other fatness QTL, including abdominal fat percentage and abdominal fat weight. ALG9 and ASAH1 are associated with chicken growth and development in the present study, providing further evidence to support the previous finding. However, no fatness QTL was found colocalized with SNPs in S100A16. Previous research demonstrated that a higher expression of the deubiquitination gene Zinc Finger RANBP2-Type Containing 1 (ZRANB1) boosted water-holding capacity via enhancing protein stability (Huynh et al., 2013;Ponsuksili et al., 2010;Wu et al., 2020).
As a comparison, an earlier study found that broiler meat had a much higher average water-holding capacity (13.50 ± 0.05%) than layer meat (11.82 ± 0.19%;Saleem et al., 2017). In accordance with these findings, the present study also identified the ZRANB1 as a functional candidate gene on chromosome 6. However, no meatquality QTL was colocalized with ZRANB1.
Recently, a study reported that Gastric Inhibitory Polypeptide (GIP) plays a key role in bone size, growth, and development in chickens Tan et al., 2021;Wang et al., 2020).
Similarly, research has discovered that GLP1 and GIP are related to synbiotic activity in broiler chickens (Seino et al., 2010). In growing broiler chickens, an injection of synbiotics could modify growth traits, immunological traits, and developmental traits on a broad spectrum (Madej et al., 2015;Madej & Bednarczyk, 2016). In this study, GIP and GLP1R were identified as candidate genes that differed strongly between BR and LB and were located within the body weight QTL regions on chromosome 27  and chromosome 3 , respectively. These QTL regions also overlapped with other QTL regions such as average daily gain, body weight, and body weight gain. We speculate that GIP and GLP1R are two essential functional candidates for studying differences in growth and immunological traits between broilers and laying hens.
A recent transcriptome study indicated the role of the Rac GTPase activating protein 1 (RACGAP1) gene in muscle growth and development by comparing two divergent broiler chicken groups, including a fast-growth group and slow-growth group and their results indicated that RACGAP1 was differentially expressed between the two groups and involved in many biological processes related to muscle growth and development (Chen et al., 2019). In the present study, we shortlisted RACGAP1 as a candidate gene that was differentiated between LB and BR and we propose that the RACGAP1 gene may play a role in the faster muscle growth performance of the broiler compared with the layer type. A study indicated that PHD Finger Protein 20 (PHF20) gene positively regulates bone formation, further PHF20-null mice showed a defective skeletal phenotype and complete knockout of PHF20 indicated spinal bone defects . Interestingly, in this study, PHF20 was found to be located in the QTL region for growth previously found on chromosome 20 (Wahlberg et al., 2009). Earlier, it has been shown that compared with layers and traditional lines, broiler tibiae were stronger, stiffer, and had lower stress values (Hocking et al., 2009). In this study, we identified PHF20 as a candidate gene that was differentiated between LB and BR and we suggest further studies to investigate its role in the chicken osteoblast differentiation. In addition, candidate genes such as mitogen-activated protein kinase 3 (MAPK3), which is involved in follicle development in chickens , or the expression of HspB1, which prolongs lifespan and increases resistance to heat stress (Alexander et al., 2022), have also been identified as functional candidate genes that differ between LB and BR.

| Genes identified as candidates distinguishing between LSL and BR strains
Comparing the LSL vs. BR, our pairwise fixation index analysis identified 18 genes that were selected as candidates and highly differentiated between the Lohmann Selected Leghorn (egg-type) and Cobb 500 Broiler (meat-type) chickens ( Figure 5). We found three candidate genes KCNK17, ENSGALG00000053358, and VPS9D1 that were common in comparisons between LSL vs BR and LB vs LSL and and were located in QTL regions for various traits. The eggshell formation is largely dependent on ion transport and Na+, Ca2+, and K+ ion channels are responsible for transporting Ca2+ from the plasma into the uterine lumen (Benos & Stanton, 1999). It has also been reported that KCNK17 is involved in transmembrane ion transport (Li et al., 2019 (Jing et al., 2014). Folic acid supplementation has been reported to improve egg quality and immunity in laying hens.
In the current study, we identified SLC19A1 as a functional candidate gene that was well known for folic acid transmembrane transporter activity. The Poly(A) Binding Protein Cytoplasmic 1 (PABPC1) gene has been found to have a significant and positive correlation with breast muscle weight (Kang et al., 2021). Another study reported that increased PABPC1 expression, as an active regulator of protein synthesis rates, in adult hearts increases the heart size and the ratio of heart to body weight (Chorghade et al., 2017). Similarly, our study also identified the PABPC1 as a candidate gene and differentiated between layer and broiler, therefore, we speculate that PABPC1 has a comparable impact on skeletal muscle growth and also indicates higher muscle growth in broiler than layer. PABPC1 was also found in the QTL region of average daily gain and breast muscle weight.
A study demonstrated that compared with wild-type mice, FABP5 knockout mice demonstrated a 24% reduction in body fat mass and suggested that FABP5 induced fat accumulation (Shibue et al., 2015).
FABP5 was found in QTL regions of fatness and growth, particularly, for average daily gain, body weight, breast muscle weight, and fat distribution traits. In the present study, we also identified FABP5 as a candidate gene that was differentiated between white layer and broiler and highlighted the role of FBAP5 in fat accumulation and body weight gain, we assume that this candidate gene may play a vital role in the higher growth rate in the broiler.

| Allele-specific expression (ASE) in the gut mucosa of laying strains
Allele-specific expression (ASE) in the gut mucosa of two highyielding layer lines was used to determine cis-regulatory divergence. In a previous study, it was shown that ASEP, a gene-based method for ASE detection, can be used to summarize information across individuals at the population level and SNPs within the same genes using RNA-seq data (Fan et al., 2020). ASE genes in each group were enriched in energy metabolism pathways, especially for the LB strain at the 24th week. It also emerges from this study that the LB strain not only has a higher body weight (Sommerfeld, Omotoso, et al., 2020) and a higher transcription level of genes Ponsuksili et al., 2021) located in energy metabolism pathways but also that the genes in this pathway are ASE genes. In addition, ASE genes were found to be particularly enriched in cholesterol biosynthesis at the 24th and 30th weeks, the peak of the laying period, which is consistent with previous findings (Ponsuksili et al., 2022). Other interesting ASE genes identified in the prelaying period were enriched in inositol phosphate metabolism. All these apparent ASE genes point to specific biological processes related to the metabolic and nutritional requirements, especially energy and fat metabolism and phosphorus homeostasis, of egg production in both strains.
The enrichment analysis of ASE genes that were strain-and agespecific ( Figure 8) revealed their involvement in specific metabolic pathways. As shown in a previous study, egg-laying activity alters the immune system toward a more pronounced humoral and innate immune response . On the one hand, cytotoxic T lymphocytes (CTL) and natural killer (NK) cells are important immune effectors related to granzyme B signaling, whereas tumor suppressor protein p53 (TP53) is an important transcriptional regulator that responds to a variety of cellular stresses and induces apoptosis or cell cycle arrest of damaged cells. ASE genes of LSL at 24 weeks were enriched in these signaling pathways. ASE genes in the LB 24- week group were enriched in the STAT3 pathway, which deals with signaling messengers and transcription factors and is involved in normal cellular responses to cytokines and growth factors. This additional evidence of ASE genes in such pathways supports the previous finding that the immune response differs in LB hens and LSL hens and that the cellular arm of the immune system is more prominent in LB . In the LSL strain, we found ALPL, encoding the alkaline phosphatase enzyme and involved in bone mineralization, as an ASE gene. This result is consistent with previous studies reporting that the two strains of laying hens differed in body weight, immunity, bone metabolism, and phytate usage (Habig et al., 2012(Habig et al., , 2014Sommerfeld, Omotoso, et al., 2020). In addition, higher ALP activity was found in the blood plasma of the LSL strain than in the LB strain in our previous study . Furthermore, at week 24, only the ASE genes within the LB strain were involved in the metabolic pathways, including glucose and glucose-1-phosphate degradation. Also, hexokinases (HK) belong to these pathways and catalyze the first step of glycolysis. HK2 is highly expressed in gut epithelium under microbial regulation via short-chain fatty acid (SCFA) and is involved in immune responses and inflammation (Hinrichsen et al., 2021). We found HK1 and HK2 as specific ASE genes in the gut of LB at the 24thweek.
Interestingly, 7 of 17 genes with Z (F ST ) ≥ 4 when comparing LB and LSL were ASE genes in at least one of the strain and age groups.
For example, we found KCNK17 and SLC13A2 (FDR < 3 x 10 −6 ) to be ASE genes at the laying peak (30 weeks) of the LB strain. As mentioned above, KCNK17 and SLC13A2 are involved in intestinal transport. These results suggest that the selection process is also associated with allelic imbalance at a particular time point of metabolic process and requirements, which is an essential step to decipher the genotype-phenotype map or functional diversity among chicken populations.

| CON CLUS ION
Our study provides several key candidate genes that were differentiated between three different chicken populations based on pairwise fixation index analysis. A number of these candidate genes are associated with characteristics such as growth, embryonic development, metabolism, immunology, and other traits of economic importance. Furthermore, Allele-Specific Expression (ASE), in the gut mucosa of two high-yielding layer chickens, suggests that the genetic architecture of different breeds may involve allelic heterogeneity, with multiple variants modifying the regulatory properties of different enhancers and regulating multiple genes and providing the baseline for further studies related to ASE for breed selection.
Genetic differences between laying hen strains LSL and LB, as well as ASE patterns observed in the laying hen population in the current study, provide insight into the complex processes underlying selection in modern commercial breeding programs.

ACK N OWLED G M ENTS
This work was financially supported by the Deutsche