Identification of SNPs Associated with Stress Response Traits within High Stress and Low Stress Lines of Japanese Quail

Mitigation of stress is of great importance in poultry production, as chronic stress can affect the efficiency of production traits. Selective breeding with a focus on stress responses can be used to combat the effects of stress. To better understand the genetic mechanisms driving differences in stress responses of a selectively bred population of Japanese quail, we performed genomic resequencing on 24 birds from High Stress (HS) and Low Stress (LS) lines of Japanese quail using Illumina HiSeq 2 × 150 bp paired end read technology in order to analyze Single Nucleotide Polymorphisms (SNPs) within the genome of each line. SNPs are common mutations that can lead to genotypic and phenotypic variations in animals. Following alignment of the sequencing data to the quail genome, 6,364,907 SNPs were found across both lines of quail. 10,364 of these SNPs occurred in coding regions, from which 2886 unique, non-synonymous SNPs with a SNP% ≥ 0.90 and a read depth ≥ 10 were identified. Using Ingenuity Pathway Analysis, we identified genes affected by SNPs in pathways tied to immune responses, DNA repair, and neurological signaling. Our findings support the idea that the SNPs found within HS and LS lines of quail could direct the observed changes in phenotype.


Introduction
Phenotypic diversities occur due to genetic variations within and between populations. Many genotypic variations arise as the result of genetic mutations. One of these mutations is known as a Single Nucleotide Polymorphism (SNP). SNPs are a common mutation and can be considered synonymous or non-synonymous [1]. Synonymous SNPs will either manifest within a noncoding region of a chromosome or silently within a coding region. Non-synonymous SNPs occur within coding regions of chromosomes and may be either missense or nonsense mutations. Both alter the amino acid at the point of mutation, though in differing ways. Missense mutations substitute one amino acid for another, while nonsense mutations stop translation entirely. Each type of mutation alters the protein produced [2]. These protein alterations may lead to alterations in phenotype. Due to this simple mechanism, SNPs are employed as biomarkers in studies surrounding population genetics [3].
Mitigation of stress is one of the primary goals of poultry producers. Stress has been shown to suppress the immune system's function, decrease egg production, and slow body growth in chickens [4,5]. Hormones released during stress can affect the production of cytokines and chemokines by leukocytes, which have proinflammatory functions. Immunosuppression resulting from chronic stress has also been shown to decrease the effectiveness of vaccines, which would increase the occurrence of disease with a flock. Two different approaches are used for mitigating stress in chickens: improving the production environment and genetic selection for healthier stress responses [5].
This study aimed to identify genetic variations that could influence differential stress responses using a whole genome resequencing method. Two different lines of Japanese quail bred for differential stress responses were used. SNPs resulting in amino acid changes (nonsynonymous) were identified and their potential effects on functionalities were determined.

Ethics Statement
This study was conducted following the recommended guidelines for the care and use of laboratory animals for the National Institutes of Health. All procedures for animal care were performed according to the animal use protocols that was reviewed and approved by the University of Arkansas Institutional Animal Care and Use Committee (IACUC Protocol #14012).

Birds and DNA Sequencing
For the purposes of this study, two different lines of Japanese quail-a high stress (HS) line and a low stress (LS) line-have been utilized. These two lines have been selected since the 1980s for divergent blood corticosterone response to restraint-induced stress in order to serve as models for stress responses in poultry [6]. The differences between the two lines include higher sociability, lower fear, higher body weight, increased egg production, and a reduced heterophil/lymphocyte ratio in the LS line compared to the HS line. Additionally, the LS line had decreased stress-induced osteoporosis, accelerated onset of puberty, heightened male sexual activity and efficiency, and a decrease in corticosterone levels of approximately one-third when compared to the HS line [6].
The development of HS and LS lines of Japanese quail up to the 12th generation is documented by Satterlee and Johnson (1988) [7]. Each line was selected for its plasma corticosterone response to immobilization [6,7]. An independent random mating condition has since been used for their maintenance [8][9][10]. These research lines were shipped to the University of Arkansas at generation 44 from Louisiana State University and maintained at the Arkansas Agricultural Experimentation Station, Fayetteville, AR [11].
Bird selection and DNA sequencing was performed in direct accordance with Khatri et al. (2019) [6]. Male HS and LS birds were chosen for this study as they have a more stable physiology as a result of having fewer hormonal changes when compared to females. Blood samples (3 ml) from 20 birds each from HS and LS lines were collected. Genomic DNA was purified from each sample using QiaAmp DNA mini kit (Qiagen, Hilden, Germany) following manufacturer's method. DNA quality was assessed using NanoDrop 1000 (Thermo Scientific, Waltham, MA, USA) and agarose gel electrophoresis. Twelve samples showing highest quality per line were pooled to represent each line. Library preparation and Illumina sequencing for the pooled DNA samples were performed by the Research Technology Support Facility at Michigan State University (East Lansing, MI, USA) using Illumina HiSeq 2 × 150 bp paired end read technology [6].

SNP Detection and Analysis
Detection and analysis of SNPs followed the methods of Khatri et al. (2018) [14]. The JMP genomics (SAS Institute Inc., Cary, NC, USA) program was used for filtering unique SNPs between the HS and LS lines of Japanese quail. Single nucleotide polymorphism (SNPs) occurring in both HS and LS lines were removed, leaving behind the unique SNPs for each line. To identify highly fixed and homozygous SNPs, the SNPs were filtered based on SNP percentages (SNP%) of SNP% ≥0.90 (for example, number of SNP = 9 of read depth = 10) were chosen. SNPs that induce non-synonymous changes in CDS (coding DNA sequences; protein coding) regions were chosen and unique SNPs in either HS or LS showing ≥10 read depths were selected as reliable SNPs. To reduce false positives, reliable SNPs chosen by criteria described above were confirmed by double-checking the initial assembly results with alignment view in SeqMan Pro program of Lasergene package (DNAStar) [14].

Bioinformatic Pathway Analysis
Genes containing unique SNPs from the HS and LS lines were entered into Ingenuity Pathway Analysis (version 60467501; Ingenuity Systems; Qiagen, Redwood City, CA, USA) for variant effect analyses for functional interpretation of genes including SNPs. The top 5 canonical pathways were chosen based on p-value calculated by the Fisher's exact test, which was provided by the IPA algorithm.

Sorting Intolerant from Tolerant (SIFT) Prediction
The SIFT prediction for amino acid changes to be deleterious was conducted following the instruction (https://sift.bii.a-star.edu.sg/; accessed on 1 March 2021). Briefly, since SIFT database for Japanese quail (Coturnix japonica) is not available, a new database for Japanese quail genome was generated using sequences in FASTA format for genome and GFF file retrieved from NCBI (https://www.ncbi.nlm.nih.gov/genome/?term=txid93934[orgn]; accessed on 1 March 2021) following instruction (https://sift.bii.a-star.edu.sg/sift4g/ SIFT4G_codes.html; accessed on 1 March 2021). The newly generated database was, then, annotated with .VCF file containing SNP list causing non-synonymous mutations. The SIFT score ≤ 0.05 is considered as deleterious amino acid substitution according to the instruction (https://sift.bii.a-star.edu.sg/sift4g/AnnotateVariants.html; accessed on 3 March 2021).

Genome Re-Sequencing and Distribution of SNPs
In accordance with Khatri et al. (2019) [6], we performed whole genome resequencing of pooled DNA samples from 12 birds each from HS and LS lines of the quail and produced 250 and~257 million reads of 150 bp respectively. Of those,~85 and~84 million reads were mapped to the reference genome (NCBI/C. japonica) and their respective depth of coverage reached to~41× and~42× for HS and LS ( Totally, 10,364 SNPs were uniquely found across the HS and LS line of Japanese quail with SNP% ≥ 0.90 and present within the CDS region. The HS line contained 6551 SNPs, and the LS line contained 3813 SNPs. These SNPs were further filtered to remove synonymous SNPs, leaving only those within the CDS region that would induce changes in amino acid sequences. This brought the number of SNPs to 1831 within the HS line and 1055 within the LS line (Supplementary Materials S1 and S2). This list of SNPs was further filtered to include genes that occurred uniquely within each line. There were 980 and 539 uniquely affected genes (which were used for Ingenuity Pathway Analysis) in the HS and the LS line, respectively. The SIFT prediction tool predicted 707 and 374 amino acid changes to be deleterious substitutions in HS and LS lines, respectively (Supplementary Materials S1 and S2). The 20 SNPs from each line with the highest read depths can be found in Tables 2 and 3, while the occurrence of each type of impact of the unique SNPs can be found in Table 4. A previous study of SNPs within a chicken autoimmune vitiligo model identified similar mutations in coding regions, causing changes in amino acid sequences [13].

Bioinformatics and Pathway Analyses
IPA was used to analyze connections of genes containing SNPs to canonical pathways and changes in molecular interactions.
Glycoprotein VI (GP6) Signaling, Signaling by Rho Family GTPases, Rac Signaling, Amyotrophic Lateral Sclerosis Signaling, and Synaptogenesis Signaling were the top 5 canonical pathways found to be significantly affected within the HS line (p-value < 0.01) ( The Synaptogenesis signaling and Amyotrophic Lateral Sclerosis signaling pathways found in the HS line are both tied to nervous system development and signaling (Supplementary Materials S5 and S6). Affected genes within the Amyotrophic Lateral Sclerosis signaling pathway are involved in regulation of neuronal damage and degeneration, as well as motor neuron apoptosis, while affected genes within the synaptogenesis signaling pathway are tied primarily to cytoskeletal development mechanisms, such as microtubule stabilization and synaptic spine organization and development [15][16][17][18][19].
The GP6 Signaling pathway is found to contain SNPs in both lines of quail (Figure 3a,b). Each line shows SNPs within collagen, laminin, phosphoinositide 3-kinase (PI3K), and protein kinase C (PKC) gene families. In the HS line, Protein kinase C delta (PKCδ) is also affected. PKCδ is part of a signaling pathway with proinflammatory functions within the nervous system [20]. In the LS line, talin 1 (TLN1) and inositol 1,4,5-triphosphate receptor type 1 (ITPR1) are affected. TLN1 expression mediates leukocyte adhesion to platelets [21]. ITPR1, also known as IP3R1, is a regulator of B cell apoptosis [22]. GP6 is a member of the immunoglobulin superfamily. It is expressed within platelets and their precursor megakaryocytes, and it serves as the major signaling receptor for collagen, which leads to platelet activation and thrombus formation [23,24]. GP6 is activated by a few different factors, one of which being low shear stress [25][26][27]. Shear stress is one of many factors that can drive gene expression within endothelial cells of the circulatory system, and low shear stress is often an indication of improper blood flow [28].   ber of the immunoglobulin superfamily. It is expressed within platelets and their precursor megakaryocytes, and it serves as the major signaling receptor for collagen, which leads to platelet activation and thrombus formation [23,24]. GP6 is activated by a few different factors, one of which being low shear stress [25][26][27]. Shear stress is one of many factors that can drive gene expression within endothelial cells of the circulatory system, and low shear stress is often an indication of improper blood flow [28].   ber of the immunoglobulin superfamily. It is expressed within platelets and their precursor megakaryocytes, and it serves as the major signaling receptor for collagen, which leads to platelet activation and thrombus formation [23,24]. GP6 is activated by a few different factors, one of which being low shear stress [25][26][27]. Shear stress is one of many factors that can drive gene expression within endothelial cells of the circulatory system, and low shear stress is often an indication of improper blood flow [28].     The Th17 Activation pathway, which includes genes containing SNPs in the LS line, is an integral part of cellular immune response (Figure 4). In addition to controlling differentiation of Th17 cells, this pathway produces IL-17A, IL-17F, IL-21, and IL-22 cytokines [29,30]. Both IL-17 cytokines have pro-inflammatory functions, are important regulators of neutrophil activity, and help mount antifungal immune responses alongside IL-22 [31,32]. IL-21 induces differentiation of B cells and T cells [33]. Similarly, the activation of IRF pathway is also associated with an immune response. This pathway is primarily involved in a host's response to viral infection, but its activation is nonetheless a result of a stimulated innate immune response (Supplementary Material S7) [34].
The canonical pathways observed within the HS line of quail show significant overlap, most of which is attributed to the presence of phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit delta (PIK3CD), phosphatidylinositol-4,5-bisphosphate 3-kinase regulatory subunit 3 (PIK3R3), and phosphatidylinositol-4,5-bisphosphate 3-kinase regulatory subunit 6 (PIK3R6), which are genes found to be active in all five pathways listed in Figure 1. These three genes code for a family of molecules known as phosphoinositide 3-kinases (PI3Ks). PI3Ks are part of the immune response. These kinases carry out phosphorylation within leukocytes and regulate the activity of Rho GTPases, which are important controllers of leukocyte motility [35].
Found within the Rac signaling, ALS signaling, signaling by Rho GTPases, and synaptogenesis signaling pathways is the gene p21-activated kinase 1 (PAK1), which encodes a p21-activated kinase (PAK) protein. PAK proteins form the bridge between Rho GTPases and the cytoskeleton, leading to changes in cell motility and growth. This bridge is necessary for successful thrombin-mediated platelet activation [36].
Mitogen-activated protein kinase 8 (MAPK8), neutrophil cytosolic factor 1 (NCF1), and nuclear factor kappa B subunit 1 (NFKB1) genes are found within both the signaling by Rho GTPases and Rac signaling pathways. MAPK8, also known as c-Jun N-terminal protein kinase (JNK1), belongs to the MAP kinase family and reacts to cellular stimuli to trigger early gene expression responses ( Figure 5). In mice, this specific kinase family has been linked to T cell proliferation, differentiation, and apoptosis [37]. NCF1 is a subunit of NADPH oxidase, an enzyme complex that drives production of superoxide during an immune response [38]. NFKB1 is part of a protein complex known as nuclear factor kappa B (NF-κB) that acts as a transcriptional regulator integral in both innate and adaptive immune response. Issues with NF-κB activation can lead to inflammatory diseases [39].
Within the LS line, only one gene, Janus kinase 2 (JAK2), is found within multiple pathways. JAK2 is found within both the Th17 Activation and IL-15 production pathways ( Figure 4; Supplementary Material S8). JAK2 codes for a kinase that is part of the JAK family of kinases. This kinase is integral for cell responses to interferon γ (IFNγ), a cytokine that is an important part of both innate and adaptive immune responses, thus cementing this gene within the immune response [40].
The LS line of quail contained multiple SNPs at various stages of a DNA damage response that responds to damage caused by ionizing radiation (Figure 6). Both a gene involved in the initial activation complex, Fanconi anemia complementation group F (FANCF), and genes involved in each of the end points of the pathway were found to contain SNPs. In the initial activation, FANCF was identified. BRCA2, Bloom syndrome RecQ-like helicase (BLM), SWItch/Sucrose Non-Fermentable (SWI/SNF), and E2F were genes identified at end points of the repair process. This DNA damage response mediates homologous recombination and chromatin remodeling, while also playing a role in tumor suppression [41][42][43].      Diseases and bio functions that were found to be significantly affected within the HS line (p-value < 0.05) include cancers, neurological disease, organismal injury and abnormalities, psychological disorders, and developmental disorders. In the LS line, the significantly affected (p-value < 0.05) diseases and biofunctions include cancers, and organismal injury and abnormalities, neurological disease, hematological disease, and metabolic disease (Supplementary Material S9).
The top physiological systems found to be significantly affected (p-value < 0.05) in development and function within the HS line include embryonic development, nervous system development and function, organismal development, tissue development, and cardiovascular system development and function. The top significantly affected (p-value < 0.05) systems within the LS line include nervous system development and function, tissue development, embryonic development, organismal development, and organ development (Supplementary Material S10). As the nervous system plays an integral role in stress response signaling, further analysis of its affected developmental functions was performed. The factors involved in nervous system development and functions were listed in Supplementary Materials S11 and S12. The top 10 functions and the related genes affected by SNPs in each line are listed in Supplementary Materials S13 and S14. Significantly affected (p-value < 0.01) nervous system developmental functions in the HS line include branching, growth and outgrowth of neurites, neuritogenesis, dendritic growth and branching, proliferation of neuronal cells, polarization of hippocampal neurons, transport of synaptic vesicles, migration of cerebellar granule cells, and development of neurons (Supplementary Material S13). Significantly affected (p-value < 0.02) nervous system developmental functions in the LS line include function, maturation, and development of neurons, development of central nervous system, myelination of nerves and the sciatic nerve, formation of inhibitory synapse and the brain, long-term potentiation of mossy fibers, and quantity of dendrites (Supplementary Material S14). Diseases and bio functions that were found to be significantly affected within the HS line (p-value < 0.05) include cancers, neurological disease, organismal injury and abnormalities, psychological disorders, and developmental disorders. In the LS line, the significantly affected (p-value < 0.05) diseases and biofunctions include cancers, and organismal injury and abnormalities, neurological disease, hematological disease, and metabolic dis-

Discussion
Our analysis found unique SNPs in pathways tied to a shear stress induced immune response, a DNA damage repair response network, and several other immune response pathways, including cytokine production and regulation and differentiation of immune cells. Additionally, pathways related to nervous system development were found to be significantly affected by SNPs.
Differing mutations within GP6 Signaling pathway could be an indication that the two lines of quail have different immune responses to changes in blood pressure [44]. It is possible these two lines may experience changes in blood pressure when exposed to stressors, but the LS line could more efficiently minimize the effects that would have affected endothelial cells. The location of SNPs in the BRCA2 DNA Damage Response is evidence that the LS line may have become better adapted at repairing cellular damage caused by stress. These two factors taken in tandem indicate that the two lines of quail primarily differ in how their bodies have adapted to control the impact of stress at the cellular level. The LS line, when exposed to the same stressors as the HS line, would not show the same negative impacts to productive traits should its physiology sufficiently compensate.
The LS line was found to have SNPs affecting multiple different pathways tied to the immune response. In particular, the presence of SNPs in the Th17 activation pathway and IL-15 production indicate a possibly direct effect on cytokine production. IL-17A, IL-17F, IL-21, and IL-22 are considered signature cytokines within the immune system with various functions each, including regulation and differentiation of immune cells. Changes in regulatory T-cell activity may protect the cells of the LS line of quail, as these cells suppress many harmful effects of immune responses [5].
The HS line also contained SNPs in genes related to immune function across multiple canonical pathways. GP6 signaling, Rho GTPase signaling, and Rac signaling pathways were found to be directly connected transitively through SNP mutations of PI3K genes and Rho GTPases. The connections of these kinases to leukocytes, as well as each pathway's connection to inflammatory responses, once again show an association between the immune system and the SNP mutations in these birds.
Both lines of quail were selected based on their response to restraint stress, which is a psychological stressor. Psychological stressors are detected by the nervous system, which can, in turn, activate hormonal responses within the endocrine system. The hypothalamopituitary-adrenal (HPA) axis is one example of this type of response. The HPA axis is a driving force behind the endocrine response to stress, as it initiates a hormone cascade leading to the release of corticosterone hormone [45]. The HS and LS lines of Japanese quail show divergent blood corticosterone levels, indicating a lasting effect of genetic selection on a significant change in the response of the neuroendocrine system to stress. Our findings show that SNPs are present in genes throughout several nervous system developmental networks, though the extent to which these SNPs could alter nervous system function still remains unclear. Further studies are needed to understand the relationship of these genes containing SNPs to nervous system development associated with the differential stress responses observed in the HS and LS lines of Japanese quail. Comparison of SNPs from the brain transcriptome could also provide more insight into the effect of stress on neuronal gene expression.

Conclusions
Overall, these findings seem to support the idea that the higher body weight and increased egg production observed in the LS line of Japanese quail may be attributed to differing neuronal development and a stronger protection from stress induced immune responses. When paired with recent findings that have found CNVs affected genes within nervous and endocrine system development, humoral and cell-mediated immune responses, and various metabolisms, our results support the theory that genotypic variation leads to the observed phenotypic diversity between the HS and LS lines of Japanese quail [6]. Additionally, more research directed towards the immune responses and cellular activity of each line of quail could deepen our understanding of the effects chronic stress may have on differently adapted production birds.