A novel 2 bp deletion variant in Ovine-DRB1 gene is associated with increased Visna/maedi susceptibility in Turkish sheep

Visna/maedi (VM) is a multisystemic lentivirus infection of sheep that affecting sheep industry across the globe. TMEM154 gene has been identified to be a major VM-associated host gene, nevertheless, a recent study showed that the frequency of the VM-resistant TMEM154 haplotypes was very low or absent in indigenous sheep. Thus, the present study was designed to determine other possible co-receptors associated with VM. For this purpose, DRB1 gene, which is renowned for its role in host immune response against various diseases was targeted. A total number of 151 case–control matched pairs were constructed from 2266 serologically tested sheep. A broad range of DRB1 haplotype diversity was detected by sequence-based genotyping. Moreover, a novel 2 bp deletion (del) in the DRB1 intron 1 was identified. For the final statistic, the sheep carrying VM-resistant TMEM154 diplotypes were removed and a McNemar’s test with a matched pairs experimental design was conducted. Consequently, it was identified for the first time that the 2 bp del variant is a genetic risk factor for VM (p value 0.002; chi-square 8.31; odds ratio 2.9; statistical power 0.90) in the dominant model. Thus, negative selection for 2 bp del variant could decrease VM infection risk in Turkish sheep.

www.nature.com/scientificreports/ infection-free flocks as flocks will remain open to a new infection at any time. A typical example of such case is a field trial in 1979, where thirteen VM affected commercial flocks with 17% mean seroprevalence were subjected to screening and culling program for every 6 months. After two years, at the 5, 6, and 7th screening, all flocks were found seronegative. However, at the 8th test, three VM seropositive sheep were found where the exact source of the infection remained unclear 21 .
The search for ways to implement effective selective breeding strategies against lethal and chronic diseases like VM has been gaining increased attention worldwide. Several candidate loci such as ovar-DRB1 22,23 , CCR5 24 , DPPA2/DPPA4 and SYTL3 25 , ZNF389 26 , and TLR9 27 were reported to be associated with the VM serostatus and/ or VM virus proviral load. Moreover, a case-control matched pairs experimental design revealed two haplotypes (haplotypes 1 and 4) in Exon 2 region of TMEM154 gene having a major effect according to the recessive model on reducing susceptibility to VM in North American sheep 28 . Subsequent studies have confirmed this association in North American 29,30 , German, Iranian 31,32 , and Turkish sheep 16 . The frequencies of the protective TMEM154 haplotypes, however, were reported to be significantly low or absent in some indigenous Turkish and Iranian sheep 16,32 . TMEM154 encodes a transmembrane protein highly expressed in B lymphocytes. There are evidence that TMEM154 variants associated with type 2 diabetes susceptibility in human 33,34 . Despite the effect of TMEM154 on VM occurrence has been repeatedly reported, its biological role in the host response pathway against VM is not yet known. Although TMEM154 has been identified to be a major gene regarding resistance/ susceptibility to VM, there could be other possible co-receptors that affect VM occurrence. In case the possible co-receptors are detected; they could be used in selective breeding for VM in indigenous sheep breeds.
Earlier attempts by our team to investigate the association between some previously reported VM associated genes, i.e., CCR5, ZNF389, TLR9 and VM serostatus did not reveal any significant result when a case-control experimental design was implemented (unpublished data). Here we aimed to investigate the possible link between Ovar-DRB1 and VM serostatus. Ovar-DRB1 was reported to be associated with VM in two independent studies 22,23 . DRB1 is a MHC class II gene that encodes antigen-presenting receptor glycoproteins called histocompatibility molecules and plays a crucial role in recognizing peptides of pathogens and presenting them to the T-lymphocytes that eventually triggers host immune response. Because of its function in the immune system, DRB1 gene has been associated with various diseases in sheep as well as in other mammalians (reviewed in 35 ).
In the present study, a case-control matched pairs experiment was carried out to investigate the possible association between Ovar-DRB1 variants and serostatus of VM in Turkish sheep. Furthermore, we have performed computational functional analysis to predict peptide-protein affinity and the possible effect of amino acid substitutions.
A McNemar's test was performed over 142 matched pairs to investigate the correlation between the 2 bp del mutation or the most prevalent haplotypes (10:01, 13:01, 09:02, and 08:01) and VM serostatus (Supplemental table S2). Both recessive and dominant models were tested. In McNemar's test, the sum of the discordant pairs (1;0 and 0;1) is expected to be higher than 25 for statistical significance. Accordingly, all of the tested variants were observed to have higher than 25 discordant pairs for the dominant model but not for the recessive model. A statistically significant correlation was observed between VM serostatus and haplotype 13:01 and the 2 bp del mutation according to the dominant model (Table 2), among these two variants, however, the most significant association was detected between the 2 bp del mutation and VM serostatus (exact p value, 0.005; chi-square, 6.89; odds ratio, 2.36; CI 95) upon computing the dominant model. Statistical power analysis was performed using real sample size (142 matched pairs) and percentage of discordant pairs (33%). Hence, the statistical power of the first analysis was calculated to be 0.87 (p < 0.05).
Alongside to strong LD between the 2 bp del and haplotype 13:01, it was identified that Histidine (H) or Tyrosine (Y) substitution to Threonine (T) at codon 61 (H/Y61T) was also linked to both 2 bp del and 1301 haplotype (Fig. 4). Among 124 ovine DRB1 haplotypes deposited in Immuno Polymorphism Database (IPD), the 61T amino acid variant is specific to 13 (Table 2). Our case-control matched pairs panel was also available for TMEM154 genotypes which were identified to have a major effect for the recessive model on VM resistance/susceptibility 16,[28][29][30][31][32] . Finally, to test the relative protection of the wild type DRB1 genotype compared to the 2 bp del variant in the absence of the protective effect of TMEM154, these matched pairs were sorted and analyzed again. Briefly, if at least one member of a pair is a carrier of protective TMEM154 diplotypes (1;1, 1;4, or 4;4), these pairs were removed from the data set, thus, 92 matched pairs lacking the protective TMEM154 diplotypes remained (Supplemental table S3). Consequently, the third McNemar's test for correlated proportion was performed on this data set with 39 discordant pairs (1;0 and 0;1). It was determined that 2 bp del mutation was still significantly associated with increased susceptibility, and the wild type ones were associated with relative resistance to VM despite lacking the protective TMEM154 diplotypes (exact p value, 0.002; chi-square, 8.31; odds ratio, 2.9; CI 95). The statistical power of this analyzis was calculated to be 0.90 (odds ratio, 2.9; CI 95; p < 0.05). According to our results, the 2 bp del mutation in DRB1 Intron1 was identified as a genetic risk factor in dominant model, i.e., having one or two copies of 2 bp del mutation was found to increase the risk of contracting the VM virus by 2.9 fold (Table 2).
In silico peptide-protein docking analysis was performed to predict the docking affinity between VM virus Gag protein antigenic epitope and DRB1 high-frequency haplotypes i.e.13:01, 10:01, 08:01, and 09:02. The prediction result revealed docking scores were variable among these four haplotypes and both of the two tools (HDOCK and HEPDOCK) computed relatively lower docking score for VM associated haplotype 13:01 when compared to other common haplotypes (Table 3). Additionally, functional analysis of point mutation effect was predicted for H/Y61T substitutions using two popular web applications. Accordingly, PROVEAN predicted the H61T and Y61T substitutions have "deleterious" effect while, PANTHER predicted the same substitutions to exhibit "probably benign" effect (Table 4).

Discussion
As in many other mammalians, there are several identified ovine DRB1 gene variants. According to IPD records, 124 ovine DRB1 alleles or subtypes have been deposited so far. In the present study, a broad genetic diversity in Exon 2 region of DRB1 gene was identified in Turkish sheep. There were 44 different haplotypes detected in at least two times (Supplemental table S1). The highest relative allelic diversity (the number of the detected different alleles over total alleles in each breeds) was in Chios breed (0.625) and the lowest in Kivircik breed (0.176). Exon 2 region of the ovar DRB1 gene has gained great attention of researchers having a broad genetic diversity and its key role in immune defense. Previously, the association between various DRB1 alleles and different diseases such as Cystic Echinococcosis 36 , faecal egg count of gastrointestinal nematodes 37,38 have been reported. The association between DRB1 gene and resistance/susceptibility to VM have also been reported in two different studies. In the first study, it was found that DRB1 haplotypes 04:03 and 07:12 were significantly associated with lower provirus levels of the ovine progressive pneumonia 22 , which is the counterpart of VM in North America, and in the second study, it was reported that the haplotype 03:25 was associated with the susceptibility to VM 23 . However, these reported haplotypes were not detected in the present study.
According to present findings, strong LD was detected between 2 bp del variant and 61T amino acid substitution which was found in 13:01 and 13:03 haplotypes in our genotype panel. To further investigate, fifty Ovar-DRB1 sequences consisting of 2 bp del variant in intron 1 were obtained from GenBank (Accession numbers: MG000511.1, MG000512.1, MG000515.1, MG000516.1, MG000518.1, MG000538.1 to 552.1, and MG000555.1 to 586.1). All the sequences belonged to the Djallonke and Sahelian native sheep breeds of Ghana in West Africa. These sequences were searched using BLAST on IPD server, and the results revealed that all of these sequences were compatible with the haplotype 13:02 (D′ value, 1; r-squared, 1). Archaeological evidence and retrovirus integration analyses suggest that selection for desired traits common to modern sheep first began in Fertile Crescent including the Anatolia region of Turkey, and spread to Africa, Europa, and other parts of Asia 39 . All the native breeds in the present study were carriers of the 2 bp del mutation, but improved breeds by backcrossing (Karacabey merino and Ramlic) were not. Other two research flocks, Bandirma which was bred from native Kivircik breed and SBA crosses, were also carriers of this mutation. When the GenBank records and our findings evaluated together; it can be inferred that having been found in African and Turkish native breeds, 2 bp deletion is an ancient mutation and perfectly linked to both haplogroup 13 (13:01, 13:02, and 13:03) and 61T amino acid variation. Furthermore, it could be speculated that 61T amino acid substitution, either alone or together with other codons, may be the causative variant for increased susceptibility against VM infection. However, case control studies are required with other 61T carrier haplotypes i.e. 05:01:01, 05:02:01, 05.03:01, 26:01, and 29:01 to strengthen this hypothesis. Furthermore, the LD status between 2 bp del and 61T amino acid variant on these haplotypes remained ambiguous. Despite these limitations, it was clearly demonstrated that 2 bp del mutation in the intron 1 of the ovine DRB1 is a strong predictor for the 61T amino acid variant in haplogroup 13, and significantly associated with increased susceptibility against VM. www.nature.com/scientificreports/ In a recent study, the VM resistant TMEM154 haplotypes in Turkish native sheep breeds were detected either at a very low frequency or complately absent (0 to 0.12) when compared to the improved breeds by backcrossing 16 . A similar observation was made in Iranian native sheep breeds 32 . Hence, selective breeding regarding protective TMEM154 genotypes might be a long lasting process to improve herd level genetic resistance to VM for native breeds. Alternatively, introgression or gene editing technology, i.e., CRISPR might be required for the native breeds lacking resistant TMEM154 haplotypes. Otherwise, 2 bp del mutation and linked H/Y61T amino acid substitution was found a range of 3-25% indicating that negative selection for 2 bp del mutation could provide a chance for relatively rapid genetic improvement against VM disease for some native sheep.
From the first matched pairs panel, 21 ewes were carriers of susceptible DRB1 2 bp del variant and protective TMEM154 diplotypes (1;1, 1;4, or 4;4) at the same time. According to the previous reports, TMEM154 protective diplotypes provide approximately threefold protection to VM, thus, it is expected that only 1/4 of these 21 ewes could be VM positive. Ten of these ewes, however, were VM positive which might indicate the 2 bp del variants potentially limits the protective effect of the TMEM154. Thus, further genetic improvement could be achieved using combination of these two genetic markers (DRB1 and TMEM154) for selective breeding of both native and improved sheep breeds.
Genotyping of MHC II locus, including DRB genes is highly problematic due to the excessive variation in a relatively short exon often causes sequence reading errors. Therefore, bidirectional sequencing is needed to reduce Identifying peptides of antigen is a key step in the development of therapeutic options for many chronic infectious diseases 40 . Predicting viral epitope especially is cost effective and time saving 41 . Despite bioinformatics prediction does not mimics experimental results, predicting epitopes plays important roles. The identified Gag protein antigenic epitope was relatively weakly bind a 2 bp del associated DRB1 receptor allele i.e. 13:01 compared to the other high-frequency alleles i.e.10:01, 08:01 and 09:02. The variation in binding affinity and peptide core could be due to the highly polymorphic nature of DRB1 especially the specific binding cores 42 . This finding possibly strengthen the susceptibility feature of 13:01 because of the poorly docking affinity between the antigenic epitope of VM virus and this haplotype which potentially interferes the immune response upon infection.
As observed in this study, H/Y61T substitution was strongly associated with 13 haplogroup. We hypothesized that the variant may maneuver the receptor core affinity to the viral antigenic epitope. This assumption is supported by the proclivity of peptides to bind to DRB1 of a particular allele in human 43 . Similarly, variants of peptides could define ligand specific binding capacity 44 which in turn is defined by the docking energy score. Further to that, according to PROVEAN web server result, the new variant (T) was predicted to have deleterious effect  www.nature.com/scientificreports/ than the commonly observed variants at position 61. The results obtained from PANTHER and PROVEAN servers were slightly inconsistent. This could be due to the assumptions that the databases are based on. PROVEAN takes the 75% sequence homology modeling to predict the effect of amino acid substitutions while PANTHER is based on the evolutionary history of a domain assuming a preserved region has functional importance. Additionally, PANTHER predicts the function of a query against the already existing 100 species in the databases which does not include Ovis aries. The sum of the results in terms of docking affinity and predictive effect of point mutation could give us a clue on how different alleles and variants that are harbored in the haplotypes affect the cellular immune reaction in fighting against the viral invasion.
In conclusion, this is a pioneering study in the identification of the association between the 2 bp del variation in ovine DRB1 intron 1 and VM serostatus in the presence and/or absence of the protective effect of the major gene TMEM154. The protective haplotypes of TMEM154 were previously detected at a high frequency in improved breeds such as Karacabey merino and Ramlic, and for native breeds they were in very low frequencies or absent. Conversely, 2 bp del variant of the DRB1 gene having a strong association with increased susceptibility Table 2. McNemar's statistics of common DRB1 haplotypes and 2 bp del mutation. Table abbreviations: a, indicate the model that is "at least one copy" is dominant and "two copy" is recessive model; b, total number of ewes that are carriers of interested haplotype as one and/or two copies; c, sum of discordant pairs; OR, Odds ratio; χ 2 , the McNemar's test statistic with continuity correction; d, representing haplogroup 13 (*1301 and *1303).   Host genotype and disease resistance/susceptibility association study requires the control of confounding factors such as variations in exposure intensity, exposure duration and other environmental factors such as climate, herd management, herd density, access to the pasture. Furthermore, it is required to increase statistical power to reduce false positive results potentially caused by the population stratification phenomenon. Therefore, association analysis was performed according to case-control matched pairs experimental design. A total of 151 matched pairs (151 cases and 151 controls) were constructed according to age, breed type, and flock. Briefly, a seropositive ewe was randomly matched with a seronegative one from the same flock, same breed type, and same age group. This matched pairs panel was used for haplotype phasing and LD analiysis. Due to limited animal numbers in SBA (6) and Hampshire crosses (6), Ramlic (2), and native Chios (4), these breeds were ignored, and statistical analysis was performed on the remaining 142 matched pairs ( Table 5). The age information was not available for a small proportion of matched pairs panel (4.2%). Thus, these individuals were matched according to breed type and flock. The pairs were assigned to be (1;1), (1;0), (0;1), and (0;0) i.e. both case and control are carriers of genetic risk/protective factors (1;1), the case is the carrier of genetic risk/protective factors, but control is not (1;0), the case is not a carrier of genetic risk/protective factors, but control is (0;1), and both case and control are not carriers of genetic risk/protective factors (0;0). In this experimental design, the higher rate of the discordant pairs (1;0 and 0;1) provides higher statistical power. A McNemar's test 50 for correlated proportions was performed using these matched pairs. To determine the statistical power of the experimental design, a power analysis was conducted using the G*Power software 51 . In silico analysis. The binding affinity of the DRB core to the peptide is a potential factor in presenting antigen to immune cells 52 . Identifying the viral epitope is the first step to calculate the affinity of the given antigenic peptide to the DRB cores. Gag protein epitope of VM virus was extracted from Immune Epitope Database and Analysis Resource (IEDB, http:// www. immun eepit ope. org). Gag protein epitope was used for further ligand-protein affinity analysis using HEPDOCK and HDOCK web servers 53 . Besides, the effect of the H/Y61T