Genomic association with pathogen carriage in bighorn sheep (Ovis canadensis)

Abstract Genetic composition can influence host susceptibility to, and transmission of, pathogens, with potential population‐level consequences. In bighorn sheep (Ovis canadensis), pneumonia epidemics caused by Mycoplasma ovipneumoniae have been associated with severe population declines and limited recovery across North America. Adult survivors either clear the infection or act as carriers that continually shed M. ovipneumoniae and expose their susceptible offspring, resulting in high rates of lamb mortality for years following the outbreak event. Here, we investigated the influence of genomic composition on persistent carriage of M. ovipneumoniae in a well‐studied bighorn sheep herd in the Wallowa Mountains of Oregon, USA. Using 10,605 SNPs generated using RADseq technology for 25 female bighorn sheep, we assessed genomic diversity metrics and employed family‐based genome‐wide association methodologies to understand variant association and genetic architecture underlying chronic carriage. We observed no differences among genome‐wide diversity metrics (heterozygosity and allelic richness) between groups. However, we identified two variant loci of interest and seven associated candidate genes, which may influence carriage status. Further, we found that the SNP panel explained ~55% of the phenotypic variance (SNP‐based heritability) for M. ovipneumoniae carriage, though there was considerable uncertainty in these estimates. While small sample sizes limit conclusions drawn here, our study represents one of the first to assess the genomic factors influencing chronic carriage of a pathogen in a wild population and lays a foundation for understanding genomic influence on pathogen persistence in bighorn sheep and other wildlife populations. Future research should incorporate additional individuals as well as distinct herds to further explore the genomic basis of chronic carriage.


| INTRODUC TI ON
Host susceptibility (vulnerability to infection), tolerance (ability to reduce the cost of infection), infectivity (ability to transmit infection), and resistance (ability to clear infection) influence transmission and persistence of pathogens in populations. Heterogeneity in these characteristics can produce a subset of individuals that contribute disproportionately to disease dynamics. For example, highly infectious individuals may act as superspreaders, resulting in disproportionately more new infections relative to other less infectious individuals (Lloyd-Smith et al., 2005;Paull et al., 2012). Individuals with low resistance and high tolerance, or those with incomplete resistance, may act as chronic carriers (individuals that host and transmit pathogens for long durations) allowing infection to persist in populations. There is a growing body of literature assessing pathogen infection and host genomics in wildlife (e.g., DeCandia et al., 2020;Donaldson et al., 2017), and the influence of host genetics on host susceptibility, resistance, and tolerance has been widely documented (recent examples include Batley et al., 2019;Bishop & Woolliams, 2014;Oleński et al., 2019); however, few studies have investigated host-genomic influence on chronic carriage.
Understanding the role of host genomics on pathogen carriage risk can provide insights for management of disease in humans, domestic animals, and wildlife. For example, in animals, genomic data can inform selective breeding for resistance to pathogens, which has been broadly discussed in agriculture and aquaculture (Bishop & Woolliams, 2014;Houston et al., 2020). These same approaches may be applied to selective breeding for individuals with low risk of becoming a chronic carrier. In free-living populations, identifying individuals genetically predisposed to maintaining infection could inform future research as well as management actions such as culling, captive breeding, and translocations.
One approach to understanding the impact of host genomics on carrier status is through genome-wide association (GWA) analyses.
Genome-wide association methods have been used to better understand the underlying architecture of disease (i.e., are few or many loci responsible for disease-related phenotypes?) and to identify single nucleotide polymorphisms (SNPs) that may be associated with disease phenotypes (Bush & Moore, 2012). These methods have been successfully utilized to understand immunity and resistance in wildlife disease systems (Elbers et al., 2018;Kosch et al., 2019;Margres et al., 2018). While GWA methods are typically restrained to large sample sizes, they can be adapted for individuals with shared ancestry, including parent-offspring pairs and other relatives, through family-based association designs (Ott et al., 2011), which can be useful when sample sizes are limited. Additional metrics assessing genome diversity, such as heterozygosity and allelic richness, can complement association analyses, identifying broader trends in diversity that may play a role in disease phenotypes.
Introduced disease has played a major role in reducing abundance and inhibiting recovery in bighorn sheep (Ovis canadensis) populations across western North America. Specifically, pneumonia epizootics (primary etiological agent Mycoplasma ovipneumoniae) have caused all-age die-off events of varying severities (48% mortality on average), with the most extreme resulting in 100% mortality and localized extirpation . Following outbreak events, M. ovipneumoniae may persist within a population via persistently infected and shedding individuals (often, asymptomatic; Plowright et al., 2017). Continued shedding by these chronic carriers often results in seasonal lamb exposure and mortality (20%-100%), suppressing recruitment and recovery in years following epizootics Cassirer & Sinclair, 2007). These dynamics can result in ongoing population declines and stagnation with occasional herd extinctions . It is recognized that M.
ovipneumoniae carriers have a significant impact on herd recovery (Garwood et al., 2020) Fox et al., 2015); however, the causative agent of tumors is unknown (Fox et al., 2016).
The outcome of M. ovipneumoniae infection in bighorn sheepwhether mortality, recovery, or transition to chronic carrier-is likely influenced by a multitude of factors, including pathogen strain, host immunity, host demographics (specifically, age, see Plowright et al., 2017), dosage, and environmental conditions. Here, we explore the association of host-genomic variation and chronic carrier phenotype in a well-studied bighorn sheep herd (Lostine herd, Hells Canyon) using family-based genome-wide association methodology and genome-wide diversity metrics. Long-term monitoring of the Lostine herd has identified chronic carriers, intermittent carriers, and noncarriers within the herd, and we hypothesize that host-genomic variation plays a role in carrier phenotype. A previous study of the Lostine herd found a strong association between carrier status and heterozygosity at a major histocompatibility complex (MHC) locus, suggesting genetic diversity at the MHC gene complex may influence carrier status (Plowright et al., 2017). This prior study examined relationships between carrier status and diversity at four microsatellite markers associated with four genes with immune system functions, as well as eleven microsatellite markers in "neutral" genomic regions (i.e., not associated with genes); here, we expand on that study by investigating associations between carrier status and allelic composition at thousands of SNPs across the genome. We also estimate genome-wide and MHC region diversity (heterozygosity and allelic richness) between phenotypes. Pathogen exposure, strain, and habitat conditions are similar across individuals within the Lostine herd, and therefore, this study system provides an opportunity to investigate genome-wide associations without the confounding influence of these variables. Due to the difficulty in collecting longitudinal, invasive data in free-ranging wildlife, our sample sizes are necessarily small; however, individuals within the Lostine herd have shared ancestry, and thus, we employ a family-based GWA approach. This study represents one of the first to explore associations between genomic variation and carrier status in a free-ranging species. Furthermore, our study demonstrates the potential power of family-based genome-wide association analyses to aid in management decision-making for bighorn sheep and other wildlife species threatened by disease.

| Study population: Lostine herd
The Lostine River bighorn sheep herd occupies an approximately 230 km 2 range in the Wallowa Mountains of northeastern Oregon and is part of the Hells Canyon metapopulation, which includes 18 sheep populations and spans the states of Idaho, Oregon, and Washington (~2.3 million hectares, Figure 1; Fortin & Cassirer, 2015). Bighorn sheep are native to this region but were locally extirpated by 1945.
The Lostine herd was re-established through the translocation of 20 sheep from Alberta, Canada in 1971. Since the reintroduction event in 1971, the population has been monitored annually. Following population growth post-introduction, management efforts were adopted to maintain the population at approximately 80 individuals through translocations of sheep out of the Lostine herd (Coggins, 2006). In 1986-87, a pneumonia outbreak resulted in the loss of two-thirds of the population (Cassirer et al., 1996;Coggins, 1988). By 1997, the population had made a near-full recovery reaching 92% of its pre-epidemic size. In the early 2000s, a different M. ovipneumoniae strain that was associated with pneumonia epizootics throughout the Hells Canyon metapopulation was detected in the Lostine population , and since then, the population has been stagnant to declining, estimated at 60 -90 sheep with no translocations out of the population (Oregon Dept. of Fish and Wildlife, unpublished data).
The Lostine herd is a well-mixed population-with many related individuals-that winters on a range where the herd is supplementally fed as part of management efforts and where our genetic samples were collected. Disease exposure was homogeneous throughout the population as was confirmed by M. ovipneumoniae antibody presence in all sheep included in this study (Plowright et al., 2017). This is an important advantage in our association analysis, as results can be confounded by individuals that are phenotypically disease-free due F I G U R E 1 Distribution of the Hells Canyon bighorn sheep metapopulation across Idaho, Washington, and Oregon. This metapopulation consists of 18 herds including the Lostine herd (red). Herd range geospatial data were freely available from state wildlife agencies' websites. The base map was sourced from ESRI (Esri & Nov., 2020) to lack of exposure, but exhibit "case" genotypes (i.e., individuals that have carrier genotypes but are classified as controls due to not being exposed). Further, a single M. ovipneumoniae strain was detected during the sampling period in the Lostine population, eliminating strain variation complexities (the same strain has been present since at least 2007; Cassirer et al., 2016). From 2011 to 2016, sheep (lambs, yearlings, andadults) in the Lostine herd were baited with alfalfa pellets and captured in a corral trap or chemically immobilized using a jab stick or dart gun.

| Sample collection and infection status
Individuals were aged and given a unique marker for recapture, and samples were collected including skin tissue, blood, and nasal swabs (Plowright et al., 2017). A subset (n = 52) of individuals was periodically recaptured, allowing for assessment of longitudinal infection status. The sampling regime focused on ewes, given the research interests of the study for which these data were collected (Plowright et al., 2017). Mycoplasma ovipneumoniae infection status was determined using conventional and real-time PCR on nasal swabs (see Plowright et al., 2017). Age was documented for each sheep at every capture event and is reported as a range spanning age at first capture to age at last capture (Table 1). A subset of individuals that died was necropsied postmortem, and sinus tumor presence was documented.
Infection phenotype was assigned based on nasal swab PCR results with the following classifications: (i) individuals that consistently tested negative for M. ovipneumoniae through time were given the status "negative"; (ii) individuals that tested positive and negative in consecutive captures, that is, infected and recovered, were given the status "intermittent"; and (iii) individuals that tested positive in two consecutive sampling events across two or more years were given the status of "chronic carrier" (Table 1) (Plowright et al., 2017).
Individuals had to be tested ≥2 times for inclusion in this study (see Table 1). Our main objective was to understand drivers for chronic (i.e., persistent) carrier status; thus, we combined the negative and intermittent classes into one "negative" control class, as this grouping allowed us to specifically address carriage.

| DNA extraction and RADseq library preparation
Restriction site-associated DNA sequencing (RADseq) data were available from a previous parentage study of the Lostine herd, including 52 individuals with longitudinal disease data ; NCBI BioProject ID PRJNA454718, SRA accession SRP144608). RADseq generates sequence data from loci distributed relatively randomly across the genome, enabling a survey of a large portion of the genome (Andrews et al., 2016). Samples used to generate RADseq data consisted of blood and skin tissue samples collected as described above. DNA extraction and RADseq library prep are described in . Briefly, DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen Inc., Germantown, MD, USA). RADseq libraries were constructed using 50 ng of high-molecular-weight genomic DNA. DNA was digested using the Sbfl restriction enzyme, and cut sites were ligated with biotinylated RADseq adapters containing 8 bp barcodes unique to each sample. Ligated products from all samples were multiplexed and sheared to 400 bp using a Covaris M220 Focused-ultrasonicator. Streptavidin bead washes were used to remove DNA fragments that did not have ligated adapters, and remaining fragments were prepped for sequencing using NEBNext
Parameters for alignment included the following: -sensitive (

| Observed heterozygosity and allelic diversity
Several genomic diversity metrics were estimated to investigate the relationship between diversity and carrier status. Standard heterozygosity (proportion of heterozygous loci divided by the mean expected heterozygosity of typed loci) was calculated for each individual using the genhet function in R v3.5.3 (Coulon, 2010; R Development Core Team, 2016). The standard heterozygosity estimate accounts for differences in the numbers of genotyped loci among individuals. In addition, several population-level diversity metrics were calculated to compare carrier and control phenotypes. Observed heterozygosity was estimated for each phenotype in R v3.5.3 using the basic.stats function in package hierfstat v0.5-7 (which accounts for group sample size; Goudet & Jombart, 2015).
We also employed a rarefaction approach to estimate metrics of overall allelic richness (the average number of alleles per locus) and private allelic richness (number of unique alleles in a population) for both carrier and control phenotypes using program ADZE v1.0 (Szpiech et al., 2008). This rarefaction approach allows for comparisons among groups where sample sizes differ, as is the case for our carrier and control groups. The maximum standardized sample size was set to 6, which is the sample size of the carrier group. For each diversity metric, we used Welch's two-sample t test in program R v3.5.3 to assess if mean diversity differed between control and carrier groups across all loci. We also implemented each of these diversity analyses to compare heterozygosity, allelic richness, and private allelic richness between phenotypes for loci that fell within the Ovar-Mhc region-where the major histocompatibility complex (MHC) class I through III genes are located-on chromosome 20.
Ovar-Mhc SNPs were identified from base 7,166,117 to 27,665,927 spanning the locations of the DYA to OLA genes based on our current understanding of the Ovar-Mhc region structure and location in domestic sheep (Dukkipati et al., 2006). Mycoplasma positive, n = 6) and analyzed using the following model (Zhou & Stephens, 2012:

| Association test
where y is a vector of binary disease labels (0 for negative, 1 for carrier), W is a matrix of covariates including a column of ones (no additional variables included in this analysis), α is a vector of corresponding coefficients (including intercept), x is a vector of SNP genotypes, β is the corresponding effect of SNPs, u is a vector of random effects, and ϵ is a vector of errors. The random effect term (u) can incorporate a relatedness matrix (centered relatedness matrix generated in GEMMA; Data S2), which we included to account for relatedness as we observed many closer-than third-degree relatives in this population (see also the kinship coefficient matrix from PLINK v2.0; Data S2). The relatedness matrix was estimated using all sheep with RADseq data in the Lostine population (n = 82; , not just the individuals with longitudinal disease information used in the GWA. The association analysis was performed using GEMMA v0.98.1 (using -miss 0.1 -hwe (1) y = W + x + u + 0.000001 -maf 0.01-km 1 -lmm 4; Zhou & Stephens, 2012), and the likelihood ratio test p-values were used to test the null hypothesis (β = 0) for each SNP. The p-value significance threshold was adjusted using Bonferroni correction resulting in a threshold of 4.71 × 10 -6 .
Further, as the Bonferroni method can be overly conservative, we controlled for type I errors by using false discovery rate (FDR) methodology, whereby new significance values (q-values) and threshold are defined: in this case, we use α = 0.15 (15% of significant values are expected to be false positives) (Benjamini & Hochberg, 1995). Manhattan and quantile-quantile plots (Q-Q plot) were created using the qqman package v0.1.4 in R (Turner, 2014).

| Identification of candidate genes
Genes with potential for influencing carrier status were identified as any genes located within a 500 kb window (up-or downstream) of SNPs associated with persistent carriage. This window is based on linkage decay reported for this RADseq dataset . Candidate genes were visualized using the NCBI Genome Data Viewer (https://www.ncbi.nlm.nih.gov/genom e/gdv) and the domestic sheep genome (Ovis aries, Oar v4.0; GenBank accession GCA_000298735.2; annotation release 102).

| Gene ontology enrichment analysis
We implemented an overrepresentation test using PANTHER (http://geneo ntolo gy.org/; Mi et al., 2018), which identifies over-or under-represented gene ontology terms for a gene subset in relation to a reference gene set. The subset of genes associated with the significant SNPs identified by the family-based GWA analysis (the candidate genes identified by the previous section) was compared to a reference gene set comprised of all genes within a 500 kb window of all SNPs used in the analysis. To determine the list of genes in the reference set, we first used the UCSC Genome Table Browser database to obtain a full list of annotated genes with their positions on the O. aries v4.0 genome (https://genome.ucsc.edu/cgi-bin/hgTables; Karolchik et al., 2004). We then used BEDTools v2.29 (Quinlan & Hall, 2010) to identify genes positioned within a 500 kb window of all loci used in the family-based GWA. The PANTHER analysis utilized biological processes, cellular components, and molecular functions identified for genes in the human genome and applied Fisher's exact test, and p-values were corrected using false discovery rate (FDR) calculations.

| Genetic architecture
To understand if pneumonia carrier status is polygenic (many loci influencing the phenotype, each having little effect) or oligogenic (few loci affecting the phenotype, each having a large effect), a Bayesian sparse linear mixed model (BSLMM) was implemented in GEMMA v0.98.1 (Zhou et al., 2013). The structure of the model was as follows: where 1 n is n-vectors of 1s, μ is a scalar for phenotype mean, x is a vector of SNP genotypes, β is the corresponding genetic marker effects, u is a vector of random effects (relatedness matrix), and ϵ is a vector of errors. The BSLMM estimates the proportion of variance in phenotypes explained by the SNPs using xβ and u (PVE; "chip heritability") and estimates a parameter rho whereby values close to 0 imply a polygenic architecture and values close to 1 imply oligogenic structure (Zhou, 2016). The probit BSLMM was run for 50 million sampling iterations following a 500k iteration burn-in, and parameter estimates were recorded every 10 iterations for a total of 500,000 sampled values.
While this method can be limited by small sample size, we report the results and acknowledge a level of uncertainty.

| Sample size and diversity metrics
A total of 98,307 variant loci were identified using genomic data from 52 bighorn sheep (19 carriers and 33 control genotypes).
After quality filtering of the sequence data (see Data S1), a total (2) y = 1 n + x + u +

| Single SNP association
The Manhattan and Q-Q plots for the carrier trait association analysis are shown in Figure 2. After implementing false discovery rate methods, the univariate linear mixed model identified two SNPs that were associated with the carrier phenotype (Figure 2c Table 2). The functions of these potential candidate genes are outlined in Table 3 and were determined using the National Center for Biotechnology Information and UniProt databases, https://www.ncbi.nlm.nih.gov/gene and https:// www.unipr ot.org, respectively (Hutchins, 2014), as well as published literature. Additional information regarding candidate gene function, processes, and disease associations are outlined in Data S4.

| Enrichment analysis
The gene ontology biological process enrichment analysis compared seven candidate genes to 1,003 reference genes identified by the F I G U R E 2 Q-Q (a), Manhattan (b), and false discovery rate (c) plots for the results of the family-based genome-wide association of the M. ovipneumoniae chronic carrier trait. (b) The Manhattan plot displays the chromosomal location and -log 10 (P-value) for each variant included in the association test (represented as points). The Bonferronicorrected threshold of 4.71 × 10 -6 (5.32 on -log 10 scale) is designated by the blue horizontal line. (c) False discovery rate corrected q-values are plotted with a cutoff of 0.15 (blue horizontal line). Two SNPs are identified as being significant (see Table 2) TA B L E 2 Single nucleotide polymorphism (SNP) loci identified as significant in the family-based genome-wide association analysis (n = 2). Allele frequencies were calculated using PLINK v1.9 (--assoc command) Note: Chr, chromosome; ID, variant identification; UTR, untranslated region; A1, allele 1 (minor allele); A2, allele 2 (major allele); f(A), frequency of allele 1 among individuals with the case phenotype (chronic carrier); f(U), frequency of allele 1 observed in individuals with the control phenotype (noncarrier); H e , observed heterozygosity at locus; LD region, genes, pseudogenes, or protein-coding regions within ± 500 k base pairs up-or downstream of the SNP location (within the linkage disequilibrium decay threshold).
USCS table browser to be ± 500 kb from any SNPs in the dataset.
The reference genes related to 8,884 biological processes (GO terms), 1,709 molecular functions, and 879 cellular mechanisms in the human genome. After FDR corrections, no gene ontology terms were found to be overrepresented in any of the three categories.

| Genetic architecture
The BSLMM revealed that 10,605 SNPs explain approximately 55.9% (PVE median, mean 54.5% ±32 S.D.) of the phenotypic variance, with a select few (median n = 16 SNPs) having a relatively large effect, explaining 45.2% (PGE median, mean 46.2% ± 31.7 S.D.) of phenotypic variance (Table 4). The value for rho (0.50, median) does not clearly identify the genetic architecture as either polygenic or oligogenic (Table 4). These results yield high uncertainty (see Table 4 and sampling distributions in Data S5), which is likely the result of small phenotypic sample size. While convergence was difficult to attain for some parameters (n[gamma] and pi), conservative inference may be made for others (PVE, PGE, and rho).

| Genetic diversity and M. ovipneumoniae carriage
Loss of genetic diversity and inbreeding are expected to leave host populations vulnerable to disease risk (Spielman et al., 2004).
Further, polymorphisms near-to or within genes involved in immune function may play an important role in an individual's susceptibility, resistance, and tolerance to disease (O'Brien & Evermann, 1988).
Here, we compared heterozygosity and allelic richness, across the entire SNP panel and for SNPs located within the MHC region, for M.
ovipneumoniae carrier and control groups in bighorn sheep.
Host heterozygosity has often been associated with genetic health and the ability to resist pathogen infection (e.g., heterozygous advantage; Spurgin & Richardson, 2010). This idea has been previously demonstrated in free-living sheep. For example, inbred Soay sheep (Ovis aries)-as assessed by heterozygosity at several microsatellite loci-were found to be more susceptible to gastrointestinal parasites (Coltman et al., 1999). In the bighorn sheep population studied here, differences in heterozygosity between M. ovipneumoniae carriers and noncarriers were previously assessed at eleven neutral and four putatively adaptive microsatellite loci (Plowright et al., 2017). Plowright et al. (2017) found that one of the four putatively adaptive Allelic richness-or, the average number of alleles per locus-is another metric by which genetic diversity can be assessed and is a reliable predictor of a population's potential adaptability (Caballero & García-Dorado, 2013). A related metric, private allelic richness (the number of unique alleles in a group) is used to genetically distinguish groups from one another (Kalinowski, 2004). We assessed allelic richness in our dataset using rarefied methods to account for disparities in sample sizes and found no difference in overall allelic richness or private allelic richness between phenotypes. The lack of significant differences for genetic diversity observed between TA B L E 3 Description of candidate gene function and tissues where they are expressed in O. aries. Information regarding the tissues where the highest expression of each gene is observed in O. aries was obtained from the National Center for Biotechnology Information website (https://www.ncbi.nlm.nih.gov/gene/). Additional information regarding biological processes for each gene can be found in Data S4

Highest expression in O. aries Citation
GHSR: growth hormone secretagogue receptor (1) Stimulation and release of hormones, (2) regulation of energy metabolism and food intake, (3) regulation of cell growth and survival, (4) pancreatic function, (5) regulation of immune functions related to aging and gastrointestinal homeostasis, (6) mitigation of inflammatory processes, and (7)  phenotypes for all diversity measures suggests that genome-wide diversity is likely not predictive of carrier status.

| Candidate genes and sinus tumors
Seven candidate genes were identified nearby to the SNPs associated with chronic carrier status, using the annotation of the domestic sheep  (Fox et al., 2015). Specifically, 36% of sheep with sinus lining thickened >5 mm were PCR positive for M. ovipneumoniae (as opposed to 13% PCR positive with no tumor presence). These findings suggest that sinus tumors may disrupt the host's ability to clear pathogens, allowing for maintenance of chronic infections in the upper respiratory tract and potential prolonged shedding (Fox et al., 2015). In our study, six of the seven candidate genes associated with persistent pathogen carriage have functions in tumor regulation, potentially indicating that a genetic predisposition to tumor growth influences pathogen carriage. However, no gene ontology terms were found to be significantly overrepresented for these seven candidate genes.
Additional research would help assess the potential for a genetic pre-

| Mutation in the growth hormone secretagogue receptor (GHSR) gene
Two SNP variants were identified as being associated with M. ovipneumoniae carrier status in bighorn sheep. One variant was located in the 5' untranslated region (5'UTR; mRNA leader sequence) of the growth hormone secretagogue receptor (GHSR) gene. The GHSR gene encodes a G-protein-coupled receptor that binds ghrelin, a peptide hormone, and is most widely recognized for its role in energy metabolism and growth (Müller et al., 2015). However, ghrelin and its receptor GHSR may play a complex role in the regulation of several other physiological functions, including immune function (Dixit & Taub, 2005;Yin et al., 2014). The variant identified by the genome-wide association analysis falls within the 5' untranslated region of the GHSR, a region that does not directly translate into proteins; however, both 5' and 3' untranslated regions are transcribed and function in post-transcriptional regulation (Hubé & Francastel, 2018). In instances when a start codon is located in the 5' untranslated region, peptide (short chain of amino acids) translation can occur, which may affect expression of down-stream proteins (Calvo et al., 2009) In domestic bovids, mutations in the GHSR gene have been of interest for selective breeding potential and quality assessment.
Polymorphisms identified in the 5' untranslated region of the GHSR were associated with increased carcass weight and average daily gain in domestic cattle (Komatsu et al., 2011). In caprine, investigations of polymorphisms in the protein-coding region of the GHSR gene have been associated with changes in body metrics, including body weight, body length, blood cholesterol, and abdominal fat in domestic sheep (Bahrami et al., 2012), and body weight and chest depth in domestic goats (Da et al., 2013). Past research has indicated that ghrelin and GHSR may play an important role in inflammation regulation (suppression; Dixit & Taub, 2005), and GHSR mRNA expression in human leukocytes (T and B cells and monocytes) suggests the ghrelin/GHSR pathway may have a role in generating or controlling immune response (Dixit & Taub, 2005;Taub, 2007). In the Lostine bighorn sheep herd, neither carcass trait data nor immune function data were available, limiting our understanding of the impact that the GHSR 5' untranslated region mutation may have had on bighorn sheep body metrics and immune response.

| Genetic architecture of carriage-polygenic?
Genetic architecture refers to the genetic contribution to-or heritability of-a given phenotype, specifically whether the phenotype is influenced by few or many variants and the respective contribution of each variant to the phenotype. For example, in wild bighorn sheep, narrow-sense heritability of fitness-related traits, specifically horn size, has been investigated using quantitative trait loci (QTL) (Johnston et al., 2010;Poissant et al., 2012). However, heredity of disease-related traits-such as susceptibility, resistance, and infectivity-can be difficult to quantify in wildlife. An understanding of genetic architecture for disease phenotypes can provide important insight into effective mitigation practices, including risk assessment, drug development, and screening (Timpson et al., 2018

| CON CLUS IONS
Here, we investigated the influence of genomic composition on persistent carriage of M. ovipneumoniae in bighorn sheep. To our knowledge, this is the first study to investigate the connection between genomic composition and persistent carriage in a free-living species using restriction site-associated DNA sequencing (RADseq) technology. We observed no difference in genome-wide diversity (allelic richness and heterozygosity) between phenotypes. However, we did identify two variant loci and seven genes associated with carrier status, and six of these genes have functions related to tumor regulation. This result suggests potential causative links between genomic composition, tumor susceptibility, and pneumonia carriage. Support

ACK N OWLED G M ENTS
We thank Digpal Gour and Jennifer Adams for assistance with laboratory work in generating the RADseq data in the previous study that were used here. We thank Oregon Department of Fish and Wildlife, Pat Matthews, Holly E. Tuers-Lance, and Chad Dotson for field as-

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.
Sequences are demultiplexed and quality filtered (using process_ radtags in STACKS).