Genes involved in immune, gene translation and chromatin organization pathways associated with Mycoplasma ovipneumoniae presence in nasal secretions of domestic sheep

Mycoplasma ovipneumoniae contributes to polymicrobial pneumonia in domestic sheep. Elucidation of host genetic influences of M. ovipneumoniae nasal detection has the potential to reduce the incidence of polymicrobial pneumonia in sheep through implementation of selective breeding strategies. Nasal mucosal secretions were collected from 647 sheep from a large US sheep flock. Ewes of three breeds (Polypay n = 222, Rambouillet n = 321, and Suffolk n = 104) ranging in age from one to seven years, were sampled at three different times in the production cycle (February, April, and September/October) over four years (2015 to 2018). The presence and DNA copy number of M. ovipneumoniae was determined using a newly developed species-specific qPCR. Breed (P<0.001), age (P<0.024), sampling time (P<0.001), and year (P<0.001) of collection affected log10 transformed M. ovipneumoniae DNA copy number, where Rambouillet had the lowest (P<0.0001) compared with both Polypay and Suffolk demonstrating a possible genetic component to detection. Samples from yearlings, April, and 2018 had the highest (P<0.046) detected DNA copy number mean. Sheep genomic DNA was genotyped with the Illumina OvineHD BeadChip. Principal component analysis identified most of the variation in the dataset was associated with breed. Therefore, genome wide association analysis was conducted with a mixed model (EMMAX), with principal components 1 to 6 as fixed and a kinship matrix as random effects. Genome-wide significant (P<9x10-8) SNPs were identified on chromosomes 6 and 7 in the all-breed analysis. Individual breed analysis had genome-wide significant (P<9x10-8) SNPs on chromosomes 3, 4, 7, 9, 10, 15, 17, and 22. Annotated genes near these SNPs are part of immune (ANAPC7, CUL5, TMEM229B, PTPN13), gene translation (PIWIL4), and chromatin organization (KDM2B) pathways. Immune genes are expected to have increased expression when leukocytes encounter M. ovipneumoniae which would lead to chromatin reorganization. Work is underway to narrow the range of these associated regions to identify the underlying causal mutations.


Introduction
Mycoplasma ovipneumoniae is considered a commensal bacterium [1] found in domesticated sheep and goats, but it can at times play a role in the aetiology of chronic, non-progressive pneumonia [2,3]. The bacteria adhere to and infect airway epithelial cells and can cause oxidative damage and apoptosis through an extracellular-signal-regulated kinase signaling-mediated mitochondrial pathway [4]. In addition, M. ovipneumoniae can moderate the host immune system by eliciting alterations in macrophage activity [5], suppress lymphocytes [6], and induce production of autoantibodies to host ciliary antigen [7].
Little is known about the lifecycle of M. ovipneumoniae in domestic sheep. At slaughter, 83 to 90% of lambs' lungs with chronic bronchopneumonia were positive for M. ovipneumoniae [8,9]. Nasal and lung isolates may not be independent populations, as a longitudinal study of New Zealand lambs found the same isolate of M. ovipneumoniae in lungs at slaughter as well as in the nasal cavity about three weeks before slaughter using older DNA technology [10]. In the U.S., 88.5% of flocks had individual sheep with M. ovipneumoniae detected in their nasal secretions (APHIS Info Sheet #708.0615). Using National Animal Health Monitoring System data, comparison of positive with negative M. ovipneumoniae flocks revealed a slight decrease in productivity but known effects that could contribute to this small difference (~4% [11]), such as litter size and ewe age, were not considered in the comparison. Seroprevalence of M. ovipneumoniae was found to be highly variable (9.76 to 30.61%) among Kazak, Hu, Merino, and Duolong sheep breeds [12]. This breed variability suggested there may be a host genetic component that is impacting bacterial shedding [12,13]. Genome-wide association studies (GWAS) can identify associated genomic regions of interest with underlying causal mutations that affect infectious disease susceptibility.
Previous exposure to M. ovipneumoniae can be assessed by the presence of antibodies [12], though this does not enable quantification of bacteria present. There are several molecular amplification assays used to detect this bacterium [26][27][28] with varying specificity depending on the genomic region of the bacterium targeted. In order to reduce false positives and ensure bacterium-specific quantification of M. ovipneumoniae, a species-specific quantitative PCR was developed for use with DNA from nasal secretions. Previous work utilizing a modified PCR [29] for screening of domestic animals and wildlife resulted in unintended or non-specific amplification of other bacteria [30]. The assay used in this study permitted a more precise measurement of M. ovipneumoniae.
The first objective of this study was to quantify longitudinal detection of M. ovipneumoniae DNA copy number from nasal secretions in a large, extensively managed flock composed of three prevalent US sheep breeds to determine patterns in detection per ewe and/or per breed. Our second objective was to identify genomic regions associated with detection of M.
ovipneumoniae DNA copy number given patterns of detection and/or breed differences identified in objective one. The breeds evaluated were developed and continue to be selected for high quality wool, fast growing, muscular lambs, or larger litters with good mothering ability to ensure wide applicability of the findings for the sheep industry. This study elucidated M. ovipneumoniae DNA copy number detection patterns over time and identified genomic regions associated with detected M. ovipneumoniae DNA copy number. Validation studies and fine mapping of markers from these genomic regions will allow for the possibility of applying selective breeding methods to reduce susceptibility to nasal shedding of M. ovipneumoniae in domestic sheep.

Ethics statement
All animal care and handling procedures were reviewed and approved by the Washington State University Institutional Animal Care and Use Committee (Permit Number: 4885 and 4594) and/or by the U.S. Sheep Experiment Station Animal Care and Use Committee (Permit Numbers:15-04, 15-05).

Sheep
Whole blood (jugular venipuncture) and nasal mucosal secretions were collected from ewes of Rambouillet (N = 321), Polypay (N = 222), and Suffolk (N = 104) breeds, ranging in age from 1 to 7 years. Repeated sampling occurred three times per year (February, April, and September/October) for 4 years (2015-2018). Ewes were managed as a single flock in an extensive rangeland production system. Sampling corresponded with mid-pregnancy (February), earlylactation (~35 days post lambing, April), and post-weaning (~30 days post-weaning, September/October) production stages. At mid-pregnancy and early-lactation sampling events, ewes had been managed in a penned, close-contact environment (e.g., dry-lot) and fed harvested feeds for at least 45 and 120 days, respectively. At the post-weaning sampling event, ewes had been managed in an extensive rangeland grazing environment. Preliminary (unpublished) data of a very small number of sheep housed under the same intensive managment conditions year-round, suggested that although some variation in the DNA copy number of M. ovipnumoniae existed, three sample time points were sufficient to detect a pattern. Therefore, a minimum of three time points was collected from each ewe but we were unable to sample every ewe at each timepoint. The number of sheep available over the four-year period was variable, as some animals were remove from the flock for other health or production reasons, each ewe was sampled 3 to 10 times, for a total of 2,876 nasal swabs collected for this study.

M. ovipneumoniae detection
Secretions were collected with sterile cotton swabs that were inserted approximately 10 cm into each nostril, consecutively, and rotated briefly. Swabs were stored without media at -20˚C until DNA was extracted with QiaAmp DNA Mini kit (Qiagen, Germantown, MD). Nasal shedding status of M. ovipneumoniae was determined by qPCR with newly designed speciesspecific primers 439F: TGATGGAACATTATTGCGCT and 439R: TGCCATTATTTGAAACRAGA. Reactions consisted of 1uL each (10uM) primer, 10uL of SsoFast EvaGreen Super mix (BioRad), 6uL of distilled water and 2uL of DNA template diluted 1:10. Each sample was run in triplicate on a CFX96 Touch system (BioRad, Hercules, CA) and quantified using a 10-fold dilution series standard curve ranging from 0 to 10 6 copies (S1 Fig). The standard curve was generated with DNA from the M. ovipneumoniae Y98 type-strain cultured in SP4 medium, pelleted by centrifugation and extracted with the QiaAmp DNA Mini kit.

Genotyping
Blood was collected into EDTA-coated vacutainer tubes. Genomic DNA was isolated using the Invitrogen GeneCatcher gDNA 3-10 ml Blood Kit (Life Technologies, Carlsbad, CA) or Gentra PureGene (Qiagen, Germantown, MD), with minor modifications of manufacturers' instructions. Quality and quantity of DNA were quantified using an ND-1000 spectrophotometer (Nanodrop, Wilmington, DE) and 300ng per ewe was sent for genotyping. Genotyping services were provided by Geneseek Inc. (Lincoln, NE) using the OvineHD BeadChip (Illumina Inc., San Diego, CA) with a set of greater than 600,000 SNP designed by the International Sheep Genome Consortium.

Phenotype analysis
Univariate analysis was conducted (SAS v9.4; SAS Inst. Inc., Cary, NC) to identify outliers and test for normality. Mycoplasma ovipneumoniae DNA copy number data were not normal, and therefore, values greater than zero were transformed to log 10 to reduce heterogenous variances. To test if there were shedding differences a reduced mixed model in SAS that accounted for fixed effects of breed, age at sample collection, year of sample collection, and approximate month of sample collection, ewe as a repeated effect, and sire as a random effect was conducted. Tukey-Kramer adjustment was implemented to identify pair-wise differences.
Quality control was performed on 647 individual samples and 606,006 SNP sites using Plink1.90 software. The following SNPs were excluded: 1) 203 for duplicate base pair location on the same chromosome; (2) 14,953 not in accordance with the Hardy-Weinberg equilibrium; 3) 31,770 with a minor allele frequency less than 1%; 4) 26,412 with genotyping rate less than 95%. In addition, multidimensional scaling was conducted to identify breed outliers. Three individuals were excluded because of failure to cluster with breed in multidimensional scaling (Fig 1). Genotyping rate was 99.8% for 644 sheep and 532,668 SNPs were used in the analyses.
Data were reformatted in PLINK1.9 for use with Efficient Mixed-Model Association eXpedited (EMMAX) [32]. The full model evaluated included fixed effects of breed or principal components 1 to 6 (eigenvalues greater than 4.9; generated with PLINK) and random effect of relationship matrix (Balding-Nichols estimate). For individual breed analysis, principal components with eigenvalues greater than 4.9 were Polypay PC1-2, Rambouilet PC1-4, and Suffolk PC1 were considered fixed effects and random effect of relationship matrix (Balding-Nichols estimate). Phenotypes considered were 1) binomial where all samples from an ewe were either greater or were less than or equal to100 DNA copy numbers, 2) categorical where a mixed detection among sampling times was added to the binomial data, and 3) untransformed mean DNA copy number. Genome-wide suggestive significance was set at a nominal P<9x10 -8 consistent with other sheep GWAS using the OvineHD Beadchip [33]. In R, qqman [34] was used for visualization of results in Manhattan and quantile-quantile plots. Further, the significant SNPs were interrogated using a mixed model similar to those described above, accounting for log 10 transformed M. ovipneumoniae DNA copy number, breed, age, year, sample time, and genotype with SAS.
Ensembl [35] and NCBI [36] were used to determine the location of the SNP within genome assembly OAR_Rambouillet_v1.0 as well as identify annotated genes within 100 kb of the SNP. The Illumina designated name of the SNP with the dbSNP rs# cluster id, the dataset where the significant SNP was identified, minor allele, and M. ovipneumoniae mean DNA copy trend associated with the minor allele can be found in S1 Table. Genotypes and phenotypes by ewe for the significant SNPs may be found at https://data.nal.usda.gov/dataset/datagenomic-regions-associated-mycoplasma-ovipneumoniae-presence-nasal-secretionsdomestic-sheep. GeneCards [37] and a search of Pubmed [36] was used to provide known information about genes within 100kb of the associated SNPs. Reactome [38] was used to determine the pathway(s) associated with the annotated genes.

Mycoplasma ovipneumoniae DNA copy number
Nasal mucosal secretions were collected from three to 10 times from each ewe, median number of collections was four. Untransformed M. ovipneumoniae DNA copy number ranged from 0 to 109,000 with a mean of 1,461 and median of 22. Transformed (log 10 ) M. ovipneumoniae DNA copy number ranged from 0 to 5.03 with mean of 1.52 and median of 1.34. There were 224 ewes with DNA copy number categorized as � 100 (zero to low detection, considered not detected) for all samplings; specifically, these ewes were sampled three to eight times with a mode of four samplings. Ninety-four ewes consistently had DNA copy number greater than 100 (considered positive for each sample time) for all sampling times; these ewes were sampled three to ten times with a mode of four samplings. Finally, there were 326 ewes with mixed levels of detection (above and below 100) when sampled over time where each ewe was sampled three to six times with a mode of four samplings. There were only 28 ewes with zero DNA copy numbers detected at all time points tested and these ewes were sampled three to six times with a mode of four samplings.
There were significant (P<0.0001) differences among breeds of log 10 transformed M. ovipneumoniae DNA copy number (Table 1). Rambouillet ewes had lower (P<0.0001) detected DNA copy numbers than Polypay and Suffolk ewes which were not different (P>0.21). Detected DNA copy number was greatest for yearlings and declined with increasing age until a slight increase at 7 years old (Table 1). Yearling ewes had significantly (P<0.046) higher detected DNA copy numbers than ewes 2 through 6 years of age. No other pair-wise age comparisons were different (P>0.08). Time during the year samples were collected significantly (P<0.0001) impacted DNA copy number (Table 1). All three sample periods were different (P<0.003) from one another where April had the highest and February had the lowest mean DNA copy number. Year of sample collection significantly (P<0.0001) impacted detected DNA copy number, with means increasing with time ( Table 1). All years were different (P<0.004) in pair-wise comparisons except for 2016 and 2017 (P>0.87).

Association analysis
No significant genomic associations were identified when M. ovipneumoniae DNA copy number detection was considered binomial, where all samples from an ewe were either positive or undetected or when categorization included mixed detection. Seven SNPs were identified as significantly (P<10 −8 ) associated with untransformed mean M. ovipneumoniae nasal shedding (Fig 2). These SNPs represented 3 regions on chromosomes 6 and 7 (Table 2) where the least frequent genotype for all seven SNPs had the highest log 10 DNA copy number and the minor allele was present in all breeds. There were 4 closely spaced SNPs on chromosome 7, approximately 21Kb apart, that were highly significant (2.72 to 6.60 Analysis of untransformed M. ovipneumoniae DNA copy number was conducted within each breed with 320, 220, and 104 ewes classified as Rambouillet, Polypay, and Suffolk, respectively (Table 3). Within Rambouillet, there were nine SNPs on five chromosomes (3, 4, 9, 10, and 22) that were significantly (<10 −8 ) associated with M. ovipneumoniae mean DNA copy number (Fig 4). Three significant SNPs on chromosome 4 were within Thyrotropin Releasing Hormone Degrading Enzyme (TRHDE) and one significant SNP was within 25Kb of TRHDE. The two significant SNPs on chromosome 22 were within Solute Carrier Family 16 Member 12 (SLC16A12). Analysis of only Polypay identified three SNPs on two chromosomes (7 and 17) that were significantly (P<4.63x10 -8 and 9.34x10 -8 , respectively) associated with M. ovipneumoniae mean DNA copy number (Fig 5). The significant SNPs on chromosome 17 were within Lysine-specific demethylase 2b (KDM2B) and Anaphase Promoting Complex Subunit 7 (ANAPC7). Nine SNPs on two chromosomes (7 and 15) were significantly associated with M. ovipneumoniae mean detection in the Suffolk-only analysis (Fig 6). Eight of these nine significant SNPs were on chromosome 15, with three of the SNPs within 31Kb of one another and two SNPs within piwi (P element-induced Wimpy)-like RNA-mediated gene silencing 4 (PIWIL4). In addition, one SNP was within Unc-13 Homolog C (UNC13C) on chromosome 15 and one on chromosome 7 was within 10Kb of Cullin 5 (CUL5).  We identified a shedding pattern among the sheep, where a subset consistently did not have detectible M. ovipneumoniae DNA, was always detected, or had mixed detection. Understanding this variability is important for conducting genomic studies and improves the possibility for selection of non-shedding sheep. Unfortunately, there were too few animals in the always detected category to identify genomic associations with the categorical analysis in this study. In other livestock, approximately 4.8% of cattle were positive (nasal secretions) for Mycoplasma bovis testing [41]. A longitudinal study of dairy cattle also found a low positivity rate of about 4% but identified one heifer that was shedding M. bovis for each time tested over 2 years  September to October had the lowest detection and April through August had the greatest. The Spanish and this study together demonstrate that transmission of M. ovipneumoniae among sheep is not necessarily just a function of proximity, but that climate may also be a factor.
Genomic associations with M. ovipneumoniae detection were identified in all-breed and individual-breed analyses. Within the all-breed analysis, significant SNPs were close to TMEM229B on chromosome 7 and PTPN13 on chromosome 6. Both these genes are found in immune pathways, and impact T cells and innate lymphoid cell type 2. Specifically, TMEM229B was found to be down regulated in memory CCR6+ T cells in an HIV study [45]. The CCR6+ T cells direct Th17 cells to effector mucosal sites [46,47], where they secrete cytokines. Some of these cytokines stimulate neutrophils and other immune cells to promote clearance of extracellular bacteria, such as Streptococcus pyogenes (Sp), Klebsiella pneumoniae, and Bordetella pertussis [48][49][50][51][52][53]. Interestingly, TMEM229B was identified as under extreme selection in chickens from tropical and arid climates [54]. The immune system of these chickens would be hypothesized to respond differently to pathogens as the two environments in the signatures of selection study were vastly different. While PTPN13 has been shown to play a role in T cell activation [55], it is specifically expressed in Th17 T cells [56] and has variable expression in T cells of people with Rett Syndrome [57]. In addition, PTPN13 is a conserved transcript in innate lymphoid cell type 2 (ILC2), which are considered helper-like and found in upper and lower airway and blood of humans [58] and mice [59]. In mouse lung, ILC2s were found to proliferate upon intranasal application of allergens, and memory-like lung ILC2s displayed upregulated expression of cytokines and their respective receptors, but also receptors for alarmins (IL25R, ST2), and genes associated with cell activation and proliferation (BCL2A1B, IER3) [60]. Also, PTPN13 (PTP-BL) interacts with the Rho effector kinase protein kinase C-related kinase 2 (PRK2) [61] and PRK2 has been demonstrated to have a role in cilia signaling [62] as well as regulation of the formation of apical junctions in bronchial epithelial cells [63]. Infection of a murine macrophage cell line with M. ovipneumoniae increases levels of nucleotide-binding oligomerization domain-containing 2 (NOD2) and NOD2 along with activation of c-Jun NH2-terminal protein kinase (JNK) induces autophagy [64]. NOD2 interacts with RAC family small GPTase 1 (Rac1) in membrane ruffles [65] and it has been suggested that Rho GTPases such as Rac1 regulate PRK2 activity during control of the cell cycle [66]. In individual-breed analyses, 21 SNPs on 8 chromosomes were identified as associated with mean M. ovipneumoniae detection. One of these SNPs was in a similar chromosomal region as the all-breed analysis, but the majority were on different chromosomes. Specifically, eight SNPs were identified on chromosome 15 within three regions in Suffolk. PIWIL4, on chromosome 15, has been identified in the chromatin pathway [67,68], where human PIWIL4 has a role in regulating gene expression in immune cells with gene silencing [69]. In vitro experiments showed that PIWIL4 knockdown in primary human cells from people with HIV-1 on suppressive combined antiretroviral therapy had increased HIV-1 transcription and reversed HIV-1 latency in HIV-1 latently infected Jurkat T cells and primary CD4+ T lymphocytes, as well as resting CD4+ T lymphocytes [70]. Mutations in CUL5, also located on chromosome 15, have been associated with more rapid loss of CD4+ T cell in HIV infected patients [71].
Two SNPs on chromosome 17, near two different genes, were significantly associated in the Polypay analysis. The expression of KDM2B has been shown to help maintain ILC3 cells, while reduction in ILC3s leads to susceptibility to bacterial infection [72]. Interestingly, optimal Th17 cell formation in mice after Streptococcus pyogenes infection requires polycomb repressive complex 1.1, which includes KDM2B [73]. The other significant gene, ANAPC7, is a part of the anaphase promoting complex/cyclosome (APC/C), which is a substantial E3 ubiquitin ligase that controls mitosis by catalyzing ubiquitination of key cell cycle regulatory proteins [74]. In human lungs, an overexpression of ANAPC7 protein has been identified and differential expression was found between smokers and non-smokers [75].
In the Rambouillet analysis, four SNPs were within or very near TRHDE on chromosome 3. This is an important regulator of thyroid hormones [76]. Although there does not appear to be a direct involvement of this gene with the immune system, a study found T-and B-cells from hyperthyroidic mice had higher mitogen-induced proliferation compared with the same cells from hypothyroidic mice [77]. Two SNPs were near an annotated gene, SLC16A12, a transmembrane transporter on chromosome 22. This gene was found to be induced by budesonide (a corticosteroid to treat asthma by reducing inflammation) in A549, BEAS-2B and primary human bronchial epithelial cells [78]. Interestingly, another GWAS examining susceptibility of mice to Streptococcal pneumoniae infection found an association with Slc16a12 as well and the authors made the correlation with Slc11a1, another known solute carrier transporter, that has been found to affect susceptibility to bacterial infections [79].
As indicated above there are other bacteria and viruses that are associated with changes in immune function genes and future studies may further support findings here. It is not unexpected that a host's immune system reacts similarly to certain organisms. Additional studies will require the qualification and quantification of the repertoire of viruses and bacteria present in nasal mucosal secretions to improve our understanding of which organisms contribute to polymicrobial pneumonia and potentially confirm the findings of this study. Further studies are also needed to quantify the strains of M. ovipneumoniae in domestic sheep and if different strains contribute to polymicrobial pneumonia.

Conclusion
Domestic sheep nasal mucosal secretions had variable mean M. ovipneumoniae DNA copy number as detected with a new species-specific qPCR. Rambouillet had the lowest mean while Suffolk and yearlings had the greatest. Time of year samples were collected affected detected DNA copy number. Twenty-eight SNPs on 9 different chromosomes were associated with mean DNA copy number. Annotated genes near these SNPs have been identified in immune function pathways (TMEM229B, PTPN13, CUL5, and ANAPC7), chromatin organization (KDM2B), and gene silencing (PIWIL4). Nasal shedding of M. ovipneumoniae may affect mucosal immunity, specifically Th17 cells. These data will be used to identify causal mutations, validate an association with improved control of M. ovipneumoniae, and causal mutations then may be developed into commercial markers for the domestic sheep industry.