**A Genome-Wide Association Study Identifying Genetic Variants Associated with Growth, Carcass and Meat Quality Traits in Rabbits**

**Xue Yang 1,2,**†**, Feilong Deng 1,3,**†**, Zhoulin Wu 1, Shi-Yi Chen 1, Yu Shi 1, Xianbo Jia 1, Shenqiang Hu 1, Jie Wang 1, Wei Cao <sup>1</sup> and Song-Jia Lai 1,\***


Received: 21 May 2020; Accepted: 17 June 2020; Published: 20 June 2020

**Simple Summary:** Rabbit meat has been widely consumed in China and is considered as an ideal food source due to its high protein, low fat, low cholesterol and low sodium contents. The growth rate, carcass characteristics and meat quality are considered economically important traits in the rabbit industry. Genomic selection (GS) could facilitate genetic selection for important economic traits, however, the lack of molecular markers for these traits limits the application of GS in rabbits. Genome-wide association study (GWAS) has the potential to comprehensively identify trait-associated molecular markers and has been applied in animal and plant research. In this study, GWAS was used to examine growth, carcass and meat quality traits of meat rabbits based on specific-locus amplified fragment sequencing (SLAF-seq) technology to identify significantly associated SNPs and functional genes, to be used as a basis for prompting the application of GS in rabbits.

**Abstract:** Growth, carcass characteristics and meat quality are the most important traits used in the rabbit industry. Identification of the candidate markers and genes significantly associated with these traits will be beneficial in rabbit breeding. In this study, we enrolled 465 rabbits, including 16 male Californian rabbits and 17 female Kangda5 line rabbits as the parental generation, along with their offspring (232 male and 200 female), in a genome-wide association study (GWAS) based on SLAF-seq technology. Bodyweight at 35, 42, 49, 56, 63 and 70 d was recorded for growth traits; and slaughter liveweight (84 d) and dressing out percentage were measured as carcass traits; and cooking loss and drip loss were measured as meat quality traits. A total of 5,223,720 SLAF markers were obtained by digesting the rabbit genome using RsaI + EcoRV-HF® restriction enzymes. After quality control, a subset of 317,503 annotated single-nucleotide polymorphisms (SNPs) was retained for subsequent analysis. A total of 28, 81 and 10 SNPs for growth, carcass and meat quality traits, respectively, were identified based on genome-wide significance (*<sup>p</sup>* <sup>&</sup>lt; 3.16 <sup>×</sup> <sup>10</sup><sup>−</sup>7). Additionally, 16, 71 and 9 candidate genes were identified within 100 kb upstream or downstream of these SNPs. Further analysis is required to determine the biological roles of these candidate genes in determining rabbit growth, carcass traits and meat quality.

**Keywords:** *Oryctolagus cuniculus*; SNPs; SLAF-seq; genome-wide association study; growth trait

#### **1. Introduction**

Rabbit meat has a long history of consumption starting from around 1100 BC [1]. It has high nutritional value and is considered a healthy food because of its high protein content and low fat, cholesterol and sodium [2]. Growth rate, carcass characteristics and meat quality are considered important economic traits in rabbit breeding. In the past few years, researchers have worked on improving growth performance and meat quality by advanced molecular breeding methods. Dozens of single-nucleotide polymorphisms (SNPs) have been identified by resequencing gene regions and found to be associated with growth traits in rabbit. Melanocortin receptor 4 (*MC4R*) [3], fat mass and obesity-associated (*FTO*) [4,5], *LEP* [6], *TBC1D1* [4] and *GHR* [7] genes are correlated with growth and carcass traits in rabbit. However, these researchers mainly investigated the correlation between a single SNP present in a specific DNA fragment with a given trait using low-throughput methods [4–6]. For complex traits such as growth performance and meat quality, large-scale analysis is necessary to detect trait-associated SNPs.

Genome-wide association study (GWAS) [8] represents a powerful approach to correlating SNPs and functional genes with quantitative traits. SNPs associated with a specific trait can be considered as molecular markers for application in genomic selection (GS) [9] and as genetic markers [10]. The most important step in GWAS is to acquire high-quality SNPs at the genome-wide level. A high-density SNP array is a high-throughput, cost-effective genotyping assay and is the most widely used genotyping method in GWAS [11,12]. Although there are still disadvantages, for example, that only known SNPs can be detected, there are high costs and great effort involved in establishing an array and that marker distribution is biased [13], it has become possible for researchers to perform GWAS using 10,000 individuals. Whole-genome resequencing is another major genotyping method that has been used over the last 10 years. It is a powerful method for whole-genome SNP discovery [14,15]. However, whole-genome resequencing can be prohibitively expensive in GWAS. Therefore, specific-locus amplified fragment sequencing (SLAF-seq), a high-resolution strategy for large-scale genotyping at the genome level, is a great alternative approach to SNP genotyping in non-popular research species such as rabbits [16]. Compared with high-density SNP arrays and whole-genome resequencing, SLAF-seq is an efficient method for de novo SNP discovery with such advantages as high genotyping accuracy, relatively low cost and a high capacity for large sample sizes. SLAF-seq has been successfully applied in chicken [17,18] and pig [19].

Despite the great success of GWAS in animal science [20], including a recent GWAS study successfully performed by Sosa-Madrid and colleagues to identify genomic regions associated with the intramuscular fat of rabbits based on a high-density SNP array [21], there is still a lack of large-scale research studies linking important economic traits to candidate genes in rabbit. Identification of SNPs associated with economically important traits via GWAS, as a first step, would provide a basis for further improving the breeding efficiency of rabbits. Here, we performed a GWAS study of the growth, carcass and meat quality traits of meat rabbits based on SLAF-seq technology to identify the associated SNPs and to predict functional genes. This study will provide a molecular basis for marker-assisted selection and gene-based selection to improve those traits.

#### **2. Materials and Methods**

#### *2.1. Ethics Statement*

All animal experiments were approved by the Institutional Animal Care and Use Committee of Sichuan Agricultural University (Permit Number: No. DKY-B20141401). The experiments were performed in accordance with the institutional regulations (no public availability).

#### *2.2. Animals and Phenotypes*

For the experiment, 53 Californian female rabbits and 22 Kangda5 line rabbits, a commercial strain of meat rabbit from Qingdao Kangda Foodstuffs Co., Ltd., were selected as the parental generation to generate crossbred F1 offspring. Rabbits with the same day of mating (2 May 2017) and the same day of delivery (2 June 2017), including 16 male Californian rabbits and 17 female Kangda5 line rabbits, were enrolled in the study. All rabbits were fed a pellet diet (10.9 MJ/kg DE, 16.5% CP, 13.3% CF) and housed in cages of 50 cm × 40 cm × 40 cm in size. During day 35–70, two rabbits were placed into one cage block and one rabbit was raised in a separate cage block from 70 days until slaughter. An environmental conditioning control system was used when the conditions were out of the normal range with respect to the ventilation (0.5–4 m3/h/kg), temperature (18–25 ◦C) and humidity (20%–70%). A 12 h period of artificial normal environmental lighting was provided. All rabbits of the F1 generation were weaned at the age of 35 days.

Bodyweights were recorded on days 35, 42, 49, 56, 63 and 70. In China, meat rabbits are usually slaughtered at 70 or 84 days old. Therefore, we recorded bodyweight at day 70 and slaughtered at day 84. A total of 432 F1 generation rabbits, including 232 male and 200 female rabbits, was slaughtered by electrical shock after fasting for 24 h. The skin and head were separated from the body by cutting at the level of the third caudal vertebra and of the distal epiphyses of radius–ulna and tibia bones. After bleeding, the hot carcasses were properly handled following procedures described in Reference [22]. The slaughter liveweight (84 d, SLW) and hot carcass weight (HCW) were measured. The dressing out percentage (DoP) was defined as HCW divided by SLW.

Two meat quality traits include cooking loss (CL) and drip loss (DL). The drip loss was quantified by the method of Reference [23]. The cooking loss was measured throughout using the following method—approximately 20 g samples of cube-like raw meat from the biceps femoris muscle of the hind leg were weighed (W1) and steamed for 30 min. Cooked samples were cooled down to room temperature and re-weighed (W2). CL was calculated as follows:

$$\text{CL (\%)}=100 \times \text{(W2/W1)}\tag{1}$$

#### *2.3. SLAF-Seq Design*

The rabbit genome (GenBank assembly accession: GCA\_000003625.1) was used as the reference genome to predict and characterize the presence of putative restriction endonuclease sites. Genome assembly revealed a genome size of 2.74 Gbp and a GC (guanine-cytosine) content of 43.75%. Restriction endonuclease profiles of the reference genome were determined using self-developed software [16]. The most efficient enzyme digestion scheme was selected based on the following criteria—(1) the proportion of restriction fragments located in the repetitive region is as low as possible; (2) zymogenic fragments are distributed evenly in the genome; (3) the number of observed restriction fragments (SLAF tags) meets the expected number of tags according to the primary in silico investigation based on reference genome sequences.

Genomic DNA was extracted from whole blood samples using the phenol–chloroform method and was digested with the selected optimal enzymes. A single "A" nucleotide was added to the 3 end of the digested fragment (SLAF tag). Dual-index [24] sequencing adapters were ligated to the A-tailed fragments using T4 DNA ligase. Further, PCR was conducted using diluted digested DNA samples, with forward primer (5 -AATGATACGGCGACCACCGA-3 ) and reverse primer (5 -CAAGCAGAAGACGGCATACG-3 ). PCR products were purified, pooled and separated by electrophoresis using a 2% agarose gel. Fragments of 300–500 bp size were excised and purified using a QIA Gel Extraction Kit (Qiagene, Germany). The library was sequenced on an Illumina HiSeq 2500 (Illumina, Inc., San Diego, CA, USA) platform following the manufacturer's instructions. Raw sequencing reads were identified by dual-indexing and classified to each sample. Clean reads were mapped to the reference genome (GenBank assembly accession: GCA\_000003625.1) using SOAP software (http://soap.genomics.org.cn) [25]. The sequences mapped to the reference genome were retained for further analysis.

#### *2.4. Genotyping and Statistical Analysis*

Efficiency of RsaI EcoRV-HF® digestion was evaluated using a positive control sample (*Oryza sativa* ssp. *japonica*). SLAF tags were mapped to the reference genome using BWA software (http://biobwa.sourceforge.net/) [26] and SNPs were identified using two methods, GATK [27] and Samtools [28]. SNPs identified using both these methods with integrity > 0.8 and MAF > 0.05 were retained for GWAS analysis.

The population structure of the rabbits was evaluated using the ADMIXTURE program [29]. The association analysis between traits and SNPs was performed according to a general linear mixed (GLM) model using PLINK2 software (http://www.cog-genomics.org/plink/2.0/) (–glm) [30,31]. The SNP effects were estimated using the following model:

$$y = Xb + Zu + \mathfrak{e} \tag{2}$$

where *y* is the vector of phenotypes; *b* is the fixed effect of sex; *u* is the SNP effect with *u* ∼ *N* -0, *I*σ<sup>2</sup> *u* ; *X* and *Z* are each the incidence matrix for *b* and *u*, respectively; and *e* is a vector of residual effects with *e* ∼ *N* - 0, *I*σ<sup>2</sup> *e* . *I* is an identity matrix; σ<sup>2</sup> *<sup>u</sup>* is the variance of SNP effects; and σ<sup>2</sup> *<sup>e</sup>* is the residual variance.

Furthermore, Bonferroni correction was applied to determine significance at the genome-wide level. SNPs with an adjusted *p*-value less than the 10%genome-wide Bonferroni-corrected threshold were annotated using ANNOVAR software [32]. Genome-wide linkage disequilibrium (LD) blocks were estimated using PLINK2 (–ld) and the LD decay is shown in Supplementary Figure S1. Based on the LD decay, the genes within 100 kb of significant associated SNPs were considered as trait-associated candidate genes. We extracted the candidate genes within the ±100 kb region of the associated SNPs according to the genome annotation information in GFF format using a script written in-house in Python3 programming language.

#### **3. Results and Discussion**

#### *3.1. Phenotype*

Bodyweight at six time points, two carcass traits and two meat quality traits were measured. The descriptive statistics for the measured phenotypic traits are listed in Table 1.


**Table 1.** Descriptive statistics of phenotypic traits.

BW: bodyweight; SD: standard deviation.

#### *3.2. SLAF-seq*

A total of 3,127,736,005 paired-end reads was generated in this study. Almost 97.07% of the raw reads were successfully mapped to the rabbit genome. The effectiveness of digestion is the key indicator of a successful SLAF-seq. In this study, 96.05% of the sequences were digested normally, which is an ideal digestion rate. A total of 5,223,720 SLAF tags was identified across the whole genome, with sequencing to a 16.57× average depth. A total of 2,498,314 polymorphic SLAF tags was identified

and 8,838,009 SNPs were selected for genome-wide association analysis based on the selection criteria (integrity > 0.8; minor allele frequency > 0.05). The density distribution of SNPs was calculated throughout the rabbit genome and is shown in Figure 1. Almost all of the genome's non-overlapped 1 Mbp region contained SNPs, which indicates that the data are reliable.

**Figure 1.** Single-nucleotide polymorphism (SNP) density distribution on chromosomes of the rabbit genome. The horizontal axis (*X*-axis) shows the chromosome length (Mbp). SNP density was calculated per 1 Mbp window. Different colors represent different SNP density levels.

#### *3.3. Population Structure Analysis*

The population structure of the rabbits was analyzed using the ADMIXTURE program. The population was first divided into 1–20 subgroups (K). The cross-validation (CV) error of populations was calculated under different k numbers. The number of the Kvalue with the lowest CV is the most suitable. Our analysis revealed that a K value of 18 was most optimal. The samples were divided into 18 subgroups (Figure 2).

**Figure 2.** The coefficient of variation for each K value. The point of the lowest cross-validation (CV) error is indicated with red color.

#### *3.4. Genome-Wide Association Analysis*

As population stratification might affect GWAS, quantile–quantile (Q-Q) plots of all traits were drawn. The observed *p* value calculated by the association study fit the expected ones, which suggests that the population stratification was well-corrected and the association analysis using GLM was reliable. The Q-Q plot of each trait is shown following the Manhattan plot of the corresponding trait.

GWAS based on SLAF-seq were successfully used to identify SNPs associated with important economic traits in chickens and pigs [18,19]. In this study, we performed SLAF-seq-based GWAS of six growth, two carcass and two meat quality traits using a General Linear Mixed (GLM) Model in meat rabbits. Table 2 shows the number of SNPs associated with each trait and the number of genes located within 100 kb upstream and downstream. Table S1 shows the location on the genome, the *p* value and the annotation of all SNPs linked with each trait.


**Table 2.** Number of SNPs associated with each trait and the number of genes around.

Based on the GLM and Bonferroni correction, 28 SNPs exhibited genome-wide association with the bodyweight trait at 35, 49, 56, 63 and 70 d (Table 2) and 16 genes near or within those SNPs were identified as associated with bodyweight. Manhattan plots for the growth traits are shown in Figures 3 and 4. Interestingly, one SNP on chromosome 8 (Chr8) was found to be associated with bodyweight at both day 49 and 56, which is located within *FGD4* (FYVE, RhoGEF and PH domain containing 4). This gene encodes a protein involved in the regulation of actin cytoskeleton and cell shape [33], which is important for cell growth. The *DNM1L* (dynamin 1 like) gene, located within 100 kb of an SNP significantly associated with BW49 and BW56, encodes a member of the dynamin GTPase superfamily, which is involved in regulating mitochondrial metabolism. The gene expression of *DNM1L* is reported to increase in skeletal muscle following exercise [34]. We thus propose these genes as potential targets for molecular breeding in meat rabbits.

**Figure 3.** Manhattan plots and quantile–quantile (Q-Q) plots of the general linear mixed (GLM)-based genome-wide association study (GWAS) for bodyweight at day 35 (**A**,**B**), 42 (**C**,**D**) and 49 (**E**,**F**). Negative log10 *p* values of the filtered high-quality SNPs were plotted against their genomic positions. The dashed lines of orange and blue indicate a 10% and 1% genome-wide Bonferroni-corrected threshold, respectively.

One SNP at position 61,268,427 of Chr9 was identified to be associated with bodyweight at day 49, 56 and 63. The SNP at position 8,681,239 in Chr8 was identified to be associated with bodyweight at day 49 and 56. SNPs at 83,815,488 and 83,815,516 of Chr13 were found to be associated with bodyweight at day 56 and 63. However, no genes were found within 100 kb upstream or downstream of these SNPs. We speculate these SNPs to have an association at two or more time points of bodyweight and be important for the growth of the rabbit and are potential targets for further molecular breeding research.

**Figure 4.** Manhattan plots and Q-Q plots of the GLM-based GWAS for bodyweight at day 56 (**A**,**B**), day 63 (**C**,**D**) and day 70 (**E**,**F**). Negative log10 *p* values of the filtered high-quality SNPs were plotted against their genomic positions. The dashed lines of orange and blue indicate 10% and 1% genome-wide Bonferroni-corrected threshold, respectively.

A total of 81 SNPs was found to be associated with two carcass traits and 71 genes were found near those SNPs. Manhattan plots for the growth traits are shown in Figure 5. Among them, one SNP located at 37,208,481 of Chr18 was identified to be associated with slaughter liveweight (84 d), with six genes located around the SNP. We identified *LIPA* (lipase A, lysosomal acid type) as the nearest gene, which encodes lipase A, the lysosomal acid lipase (cholesterol ester hydrolase). This enzyme catalyzes the hydrolysis of cholesteryl esters and triglycerides in the lysosome [35]. In addition, 80 SNPs were found associated with carcass rate and 65 genes near these SNPs encoded proteins with various functions (Supplementary Material Table S1). However, there have been no reports of the association between these candidate genes and animal production traits. Further study is required to validate the links between candidate genes and traits.

**Figure 5.** Manhattan plots and Q-Q plots of GLM-based GWAS for the dressing out percentage (**A**,**B**) and slaughter liveweight (**C**,**D**). Negative log10 *p* values of the filtered high-quality SNPs were plotted against their genomic positions. The dashed lines of orange and blue indicate a 10% and 1% genome-wide Bonferroni-corrected threshold, respectively. Red ellipses indicate that SNPs within the same ellipse share the nearest candidate gene.

For the meat quality phenotype, 15 SNPs were associated with cooking loss and drip loss. Manhattan plots for the growth traits are shown in Figure 6. Five SNPs were associated with cooking loss, with 2 genes located nearby. The nearest gene *OSBPL11* (oxysterol-binding protein like 11) encodes a member of the oxysterol-binding protein (*OSBP*) family, a group of intracellular lipid receptors [36]. Ten SNPs were associated with drip loss, with 9 genes around those SNPs. Two main types of proteins in the human body involved in the maintenance of zinc ion homeostasis—zinc binding protein, which acts as a buffer substance or as a donor of intracellular zinc and zinc transporters, which are involved in the intake and excretion of zinc in cells. Among the five genes found in close proximity to the associated SNPs of the drip loss trait in this study, two encode proteins involved in the maintenance of zinc homeostasis. *ADAM7* (ADAM metallopeptidase domain 7) on chromosome 2 encodes a member of the *ADAMs* family of zinc proteases and *SLC39A6* (solute carrier family 39 member 6) on chromosome 9 encodes a protein with structural characteristics of zinc transporters [37]. These findings indicate that the potential targets for meat quality are the genes relevant to zinc homeostasis.

**Figure 6.** Manhattan plots and Q-Q plots of the GLM-based GWAS for the meat quality traits of cooking loss (**A**,**B**) and drip loss (**C**,**D**) in rabbits. Negative log10 *p* values of the filtered high-quality SNPs were plotted against their genomic positions. The dashed lines of orange and blue indicate a 10% and 1% genome-wide Bonferroni-corrected threshold, respectively. The red ellipse indicates that SNPs in the same ellipse share the nearest candidate gene.

#### **4. Conclusions**

In this study, we performed GWAS of 432 F1 meat rabbits. All samples were genotyped using SLAF-seq technology. The GWAS of 10 economic traits revealed 111 significantly associated SNPs. A total of 98 putative candidate genes were located within 100 kb of significantly associated SNPs. Among them, *FGD4* and *DNM1L* genes, linked with BW49 and BW56 in our study, are reported to be involved in cell growth and mitochondrial metabolism. *ADAM7* on chromosome 2 and *SLC39A6* on chromosome 9 may be associated with drip loss due to its effect on zinc homeostasis that functions in determining meat quality. These findings provide novel insights into the genetic basis of growth, carcass and meat quality traits in rabbits and may contribute to the application of GS in rabbits. However, because only a draft version of the rabbit reference genome is available, there are numerous sequence gaps and missing annotations. These gaps could influence the accurate estimation of LD patterns. Moreover, the missing annotations would inevitably lead to false-negative results when attempting to find causative genes or variants. In addition, relatively high numbers of false positives were observed in the GWAS study. Therefore, a future validation study should be conducted using an independent population.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-2615/10/6/1068/s1. Figure S1. A schematic representation of genome-wide linkage disequilibrium (LD) decay. Table S1. Significant SNPs for growth, carcass and meat quality traits in rabbits.

**Author Contributions:** Conceptualization, S.-J.L. and S.-Y.C.; methodology, X.Y. and F.D.; validation, X.Y., and S.H.; formal analysis, X.Y., Z.W., Y.S. and F.D.; investigation, X.Y., X.J., J.W. and W.C.; writing—original draft preparation, X.Y and F.D.; writing—review and editing, S.-J.L.; visualization, X.Y.; project administration, S.-J.L.; funding acquisition, S.-J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the earmarked fund for China Agriculture Research System, grant number CARS-43-A-2.

**Conflicts of Interest:** All authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Myeloperoxidase and Lysozymes as a Pivotal Hallmark of Immunity Status in Rabbits**

#### **Rafał Hrynkiewicz, Dominika B ˛ebnowska and Paulina Nied ´zwiedzka-Rystwej \***

Institute of Biology, University of Szczecin, Felczaka 3c, 71-412 Szczecin, Poland; rafal.hrynkiewicz@gmail.com (R.H.); bebnowska.d@wp.pl (D.B.)

**\*** Correspondence: paulina.niedzwiedzka-rystwej@usz.edu.pl

Received: 15 August 2020; Accepted: 3 September 2020; Published: 4 September 2020

**Simple Summary:** Rabbit breeding is a very important element in the context of broadly understood industrial breeding, as rabbits are one of the main and most frequently chosen economic directions. Effective rabbit breeding, however, requires full control over the health of these animals, which is particularly related to the orientation regarding their immune status. There are many indicators that can be used to assess the immune system, but the greatest attention should be paid to those that change rapidly over time and reflect the body's first line of defense. Peripheral blood granulocytes contain enzymes with strong antimicrobial properties, the level of which changes as a result of various external factors, e.g., viral infection, which was assessed in this study. The aim of the study was to evaluate the dynamics of myeloperoxidase (MPO) and lysozyme (LZM) in the experimental infection of rabbits with the *Lagovirus europaeus*/GI.1a virus, which is a pathogen causing high mortality, decimating rabbit farms all over the world in a short time. The results obtained in the dynamic system show that the levels of assessed enzymes significantly change in the blood during infection. Assessing the immune system using these indicators could therefore be a potential biomarker for the immune status of rabbits.

**Abstract:** Infectious diseases, due to their massive scale, are the greatest pain for all rabbit breeders. Viral infections cause enormous economic losses in farms. Treating sick rabbits is very difficult and expensive, so it is very important to prevent disease by vaccinating. In order to successfully fight viral infections, it is important to know about the immune response of an infected animal. The aim of this study was to analyze the immune response mediated by antimicrobial peptides (myeloperoxidase (MPO) and lysozyme (LZM)) in peripheral blood neutrophils and rabbit serum by non-invasive immunological methods. The study was carried out on mixed breed rabbits that were experimentally infected with two strains (Erfurt and Rossi) of the *Lagovirus europaeus*/GI.1a virus. It has been observed that virus infection causes changes in the form of statistically significant increases in the activity of MPO and LZM concentration, while in the case of LZM activity only statistically significant decreases were noted. Additionally, clinical symptoms typical for the course of the disease were noted, and the probability of survival of the animals at 60 h p.i. (post infection) was 30% for the Erfurt strain, and −60% for the Rossi strain. The obtained results of MPO and LZMs suggest that these enzymes, especially MPO, may serve as a prognostic marker of the state of the immune system of rabbits.

**Keywords:** myeloperoxidase; lysozyme; rabbits; viral infection; rabbit hemorrhagic disease

#### **1. Introduction**

Rabbit breeding is an important and future-oriented direction of animal production, especially from an economic point of view. Apart from being a key element of the economy of many countries, rabbits also exist as an important link in the trophic chains of the Mediterranean ecosystems [1–5]. As laboratory animals, they are widely used in various types of diagnostic tests and analyses [2,6–8]. Currently, rabbit breeding is becoming increasingly popular. The data concerning the world production of these animals indicate that annually about 2 million of them are slaughtered exclusively for meat [9].

The loss of rabbits, especially due to various diseases, is very problematic for breeders, because even the best farms, by not following the rules of disease prevention and improper treatment, face huge problems. The treatment of sick rabbits is very difficult, as it requires cumbersome, long-lasting and systematic procedures [10,11], and the causes of rabbits' diseases can be very different. The most common diseases of rabbits are infectious diseases, which can occur en masse, thus leading to serious losses [10,11].

One of the most dangerous diseases that can decimate rabbit breeding is rabbit haemorrhagic disease (RHD) [10–12]. RHD is a highly infectious and deadly viral disease manifested by acute viral hepatitis in rabbits, caused by single-stranded RNA (ssRNA) virus ssRNA virus *Lagovirus (L.) europaeus* [13–15]. The genus *L. europaeus* is divided into two main gene groups related to rabbit hemorrhagic disease virus (RHDV): GI or European hare syndrome virus (EBHSV): GII [14]. Within the GI gene group, four genotypes—GI.1, GI.2, GI.3, and GI.4—were distinguished. The GI.1 genotype includes classical RHDV strains (*L. europaeus*/GI.1), while the GI.2 genotype was classified as RHDV2, discovered in France in 2010 (*L. europaeus*/GI.2) [16]. Genotype GI.3 is represented by RCV-E1 (*L. europaeus*/GI.3), and GI.4 by RCV-A1 (*L. europaeus*/GI.4) and RCV-E2 (*L. europaeus*/GI.4d). The strains present in these genotypes are called nonpathogenic rabbit calicivirus (RCV). The strains of the GI.1 and GI.2 genotypes are known infectious agents responsible for the development of RHD. According to the recently proposed nomenclature [14], GI.1 is further subdivided into different antigenic variants (GI.1a-d) based on phylogeny and genetic distance [14].

Due to its features, RHD has been included in the list of notifiable diseases to the World Organisation for Animal Health (OIE). Every animal in which RHD is suspected or confirmed is immediately euthanized and cremated [10,11]. The first data of RHD were reported in Wuxi, Jiangsu Province in China in 1984 [17]. The disease was observed in the population of European rabbits (*Oryctolagus cuniculus*) of the Angora breed, which were imported from the former German Democratic Republic to China for breeding purposes [17,18]. In China, within one year, the disease contributed to the deaths of over 140 million domestic rabbits and spread over about 50,000 km2 [18]. The Chinese epidemic of RHD caused enormous losses in the economic sphere of the country [18,19]. The strong expansion of the virus soon led to its rapid spread throughout Europe [2]. Since 2010, the RHDV2 virus (*L. europaeus*/GI.2) has also spread to countries in Europe, Australia, America, and Africa [20]. The appearance of RHD on the Iberian Peninsula led to huge losses in the wild rabbit population [1–3,16,21]. There was a great interest in this problem due to the rabbit's key contribution to the Mediterranean ecosystems, as it provides food for several endangered endemic species, such as the Iberian lynx (*Lynx pardinus*) and the Iberian imperial eagle (*Aquila heliaca*) [3]. The reduction of the wild rabbit population therefore has a negative impact on the persistence of these predators [3,22,23].

Currently, the disease is present in Europe, Asia, Africa, and Australia. Single epizooties are recorded from time to time in North America (United States, Canada, Cuba). According to the current information in the WAHID (World Animal Health Information Database) system [11], in the years 2005–2018, the occurrence of RHD was reported or suspected in 50 countries, of which more than half of the reports were reported in European countries [11]. Unfortunately, so far there is no effective cure that can save infected rabbits. The only correct method of fighting RHD is its prevention through prophylactic vaccinations [1,3,20,24].

Moreover, the key to administering proper vaccination, as well as efficient fighting against the virus, is the knowledge of the immune system response of the infected animal. In light of the above, the aim of the study was to analyze with a noninvasive immunological method, the immunological response mediated by antimicrobial peptides (myeloperoxidase and lysozyme) in neutrophils of peripheral blood and serum of rabbits infected experimentally with *L. europaeus*/GI.1a.

#### **2. Materials and Methods**

#### *2.1. Animals*

The study was conducted on a group of 25 Polish mixed-breed rabbits of various sexes, of the *Oryctolagus cuniculus* species, weighing in the range of 2.50–3.80 kg, marked as conventional animals, from licensed breeding under constant veterinary and zoo-hygienic supervision [25]. The rabbits were not vaccinated against *Lagovirus europaeus*, were in good health, and did not show any disease symptoms. Throughout the experiment, the animals stayed in a vivarium belonging to the Institute of Biology, the Faculty of Mathematical, Physical and Natural Sciences, the University of Szczecin. During the tests, appropriate zoo-technical conditions were ensured, in accordance with the recommended Polish standards developed in line with the European Union Directive with regards temperature and humidity, as well as lighting and size of cages for animals [26]. After transporting to the vivarium, the animals were subjected to a two-week adaptation period and tested for the presence of anti-RHDV antibodies using a ready-made ELISA (Institutio Zooprofilattico Sperimentale, Italy). Each of the rabbits had a separate metal cage, was fed with complete feed for rabbits (16% Rabbit, Motycz, Poland) at an amount of 0.15–0.20 kg per day, and had unlimited access to fresh drinking water. All studies were conducted with the approval of the Local Ethics Committee for Experiments on Animals in Pozna ´n (license no. 1/2009).

#### *2.2. The Scheme of the Experiment*

Two antigenic variants of *Lagovirus europaeus*/GI.1a with positive hemagglutination (HA+) were selected for the study:


The rabbits that were used in the experiment were divided into three groups. The animals were classified into groups based on randomness. The first group (*n* = 10) were rabbits (60% female; 40% male) infected with *L. europaeus*/GI.1a/Erfurt. The second group (*n* = 10) were rabbits (50% female; 50% male) infected with *L. europaeus*/GI.1a/Rossi (*n* = 10). The third group (*n* = 5) were rabbits (40% female; 60% male) that were classified as the control group.

The virus infection occurred after the first blood sampling (0 hour of the experiment). The virus was infected by intramuscular injection of *L. europaeus*/GI.1a antigens into the lower limb muscle. Control group (*n* = 5) rabbits were given a placebo in the form of PBS (phosphate-buffered saline) in the same way. Subsequent blood samples were taken in hours: 8, 12, 24, 36, 48, 52, 56, and 60 h post infection (p.i.). Myeloperoxidase (MPO) activity was determined in peripheral blood samples, whereas in blood serum the concentration and activity of lysozyme (LZM) was determined.

#### *2.3. Virus Preparation and Administration*

Viruses were obtained from animals that died under natural conditions in different parts of Germany (strains were obtained in lyophilized form from Dr. Horst Schirrmeier, Friedrich Loeffler Institute Greifswald, Germany). Virus identification was performed using the real-time PCR method under the conditions previously described by Nied ´zwiedzka-Rystwej et al. [27]. Liver homogenate was prepared from the recovered livers, and the animals were experimentally infected in order to increase the amount of virus. After their death, the liver was prepared as a 20% homogenate, which was purified by centrifugation (3000 rpm), treatment with 10% chloroform for 60 min, and centrifugation again. Then, a suspension in glycerol was prepared in a 1:1 ratio [28,29]. All the antigens prepared in this way had the same number of virus particles, with the density from 1.310 to 1.340 g/cm3. Additionally, the virus titer was determined using the HA hemagglutination test, which was 1:1280 for both strains.

#### *2.4. Myeloperoxidase (MPO) Activity*

Determination of myeloperoxidase (MPO) activity in neutrophils was performed according to Graham's method, as described by Zawistowski [30]. This method is based on making a smear of edetate blood, collected in a sample with an anticoagulant in the form of heparin, and the use of a color reaction with myeloperoxidase (MPO) in PMN (polymorphonuclear) cells, based on the reaction of benzidine, which in the presence of hydrogen peroxide and peroxidase passes into brown oxybenzidine. The reagent used in the color reaction with MPO was made from 250 mL of 40% ethyl alcohol, in which 350 mg of benzidine was dissolved; then the reagent was filtered and 0.35 mL of hydrogen peroxide was added. All blood smears were fixed with a mixture of formalin and ethyl alcohol, and then, for at least 3 min, were stained with a benzidine reagent. At a later stage, the preparations were stained with a 7% Giemsa solution for 7 min to stain the PMN cell nucleus. Then, the preparations prepared in this way were viewed under an optical microscope at 1000-fold magnification using immersion. In the preparations, the number of granulocytes was calculated, and the degree of color of the granules was determined according to the adopted scale (Figure 1). MPO activity was expressed as a factor, according to the intensity of the granule color in 100 PMN cells, and calculated on the basis of this formula according to Afanasyev and Kolot [31]:

**Figure 1.** Scale showing the intensity of granule staining in PMN cells.

#### *2.5. Concentration and Activity of Lysozyme (LZM)*

The concentration of LZM was determined in the blood serum towards *Micrococcus (M.) lysodeikticus* using the platelet diffusion method, according to the method of Hankiewicz [32].

The substrate was a 1% agarose gel prepared from 1 g of agarose dissolved in 100 mL of 0.067 M phosphate buffer at pH 6.2, with the addition of 1 g of NaCl, in which *M. lysodeikticus* (Sigma, Saint Louis, MO, USA) was suspended (150 mg of bacteria in 15 mL phosphate buffer, at 0.067 M pH

6.2). The agarose gel prepared in this way was poured between glass plates. The agarose gel plates were allowed to set, and then 17 μL wells were carved into the gel.

The LZM standard with a concentration of 100 mg/L was obtained by dissolving 10 mg of hen egg white (Sigma, Saint Louis, MO, USA) with activity of 50,000 IU/mg in 100 mL of 0.067 M phosphate buffer at pH 6.2. A subsequent dilution at a concentration of 64 mg/L was obtained by transferring 6.4 mL of LZM standard solution (100 mg/L) into 3.6 mL of 0.067 M phosphate buffer at pH 6.2. Subsequent concentrations of 32, 16, 8, 4, 2, 1, 0.5, 0.25, and 0.125 mg/L were obtained by serial dilution by transferring 1 mL of standard from each higher concentration to 1 mL of buffer. The thus-prepared LZM standard (17 μL) and the test serum (17 μL) were spotted into the grooved wells in an agar medium, and then the plates with the medium were incubated in a humid chamber and placed in an incubator at 37 ◦C for 18 h. After incubation, the diameter of radiolucency zones around individual wells was measured; during this process, the size of the diameter depends on the concentration of LZM in the blood serum. The results of the LZM concentration are read from the calibration curve made on the basis of standard solutions.

In turn, the activity of LZM was calculated on the basis of the formula presented by Szmigielski [33] and presented as an indicator:

$$\text{LZM activity ratio} = \frac{\text{LZM concentration in serum (mg/L)}}{\text{absolute number of neutrons} \text{ of neutrons}} \tag{2}$$

The absolute neutrophil count (\*) was calculated by multiplying the total leukocyte count (in 1000) by the percentage of neutrophils from the blood quality picture (leukogram).

#### *2.6. Clinical Studies of Experimental Animals*

Throughout the experiment, clinical symptoms and the survival of rabbitsinfectedwith*L. europaeus*/GI.1a were recorded. Rabbit survival is presented as the percentage of survival analyzed using the Kaplan– Meier method.

#### *2.7. Statistical Analysis*

All results were statistically analyzed using the Student's *t*-test with Cochran–Cox correction, assuming *p* = 0.05. The analysis was performed using in the *Statistica* software, ver. 13.1 (Statsoft, Poland) and Microsoft Excel (Microsoft 365, Redmond, WA, USA).

#### **3. Results**

#### *3.1. Myeloperoxidase (MPO) Activity*

The values of MPO activity in the group of rabbits infected with *L. europaeus*/GI.1a/Erfurt ranged from 1.15 to 2.13, with the standard deviation (SD ±) within the range of 0.08–0.29 (Figure 2). On the other hand, in the case of animals infected with *L. europaeus*/GI.1a/Rossi, MPO activity ranged from 1.07 to 2.34, with the standard deviation (SD ±) from 0.16 to 0.65 (Figure 2). In the group of control rabbits, MPO activity values ranged from 1.01 to 1.27, with the standard deviation (SD ±) from 0.02 to 0.15 (Figure 2).

Both results of MPO activity in rabbits infected with *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/ GI.1a/Rossi throughout the entire digestion of the experiment showed a clear upward trend of the parameter under study, compared to the results of MPO activity obtained in the studies of rabbits from the control group. The obtained results were subjected to statistical analysis, which showed that both *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi strains caused eight statistically significant changes in the form of increases in the parameter under study. These changes were observed after 8, 12, 24, 36, 48, 52, 56, and 60 h p.i. (Figure 2).

**Figure 2.** Values of myeloperoxidase (MPO) activity in rabbits infected with *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi, as well as control rabbits. \* statistically significant change with respect to the control group (*p* < 0.05).

#### *3.2. Lysozyme (LZM) Concentration*

The parameters of LZM concentration for *L. europaeus*/GI.1a/Erfurt infected rabbits ranged from 4.22 mg/L to 6.40 mg/L, with standard deviation (SD ±) from 0.78 to 1.91 (Figure 3). On the other hand, the results recorded for *L. europaeus*/GI.1a/Rossi infection ranged from 4.03 mg/L to 5.49 mg/L, with the standard deviation (SD ±) from 0.54 to 2.12. (Figure 3). The concentration of LZM in control rabbits ranged from 3.60 mg/L to 4.50 mg/L, with the standard deviation (SD ±) ranging from 0.20 to 0.90 (Figure 3).

**Figure 3.** Values of lysozyme (LZM) concentration in rabbits infected with *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi, as well as control rabbits. \* statistically significant change respect to the control group (*p* < 0.05).

Both of the parameters for LZM concentration in rabbits infected with *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi, in relation to the results obtained in the study of rabbits from the control group, showed an upward trend throughout the duration of the study. Nevertheless, they were accompanied by significant standard deviation. The obtained results were analyzed statistically, which showed that in the case of infection with *L. europaeus*/GI.1a/Erfurt, three statistically significant changes were noted, which fell at 12, 24, and 56 h p.i., while the *L. europaeus*/GI.1a/Rossi strain in within the examined parameter showed only one change of statistical significance; this change was observed at 12 h p.i. (Figure 3). High standard deviations and relatively low numbers of statistically relevant changes may be caused by the rapid losses of rabbits and the methodical character of the analysis, which depends on many outside factors.

#### *3.3. Lysozyme (LZM) Activity*

The results of studies on LZM activity in rabbits infected with *L. europaeus*/GI.1a/Erfurt ranged from 0.0015 to 0.0057, with the standard deviation (SD ±) from 0.0002 to 0.0047 (Figure 4). In *L. europaeus*/GI.1a/ Rossi-infected animals, LZM activities ranged from 0.0004 to 0.0014, with a standard deviation (SD ±) of 0.0001 to 0.0008 (Figure 4). In the group of control rabbits, LZM activity values ranged from 0.0013 to 0.0023, with standard deviation (SD ±) from 0.0002 to 0.0008 (Figure 4).

**Figure 4.** Values of lysozyme (LZM) activity in rabbits infected with *L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi, as well as control rabbits. \* statistically significant change respect to the control group (*p* < 0.05).

In the case of rabbits infected with *L. europaeus*/GI.1a/Erfurt, the parameters of LZM activity, in relation to the results obtained in the study of rabbits from the control group, showed a slight upward trend at 0 hours. Then at 12 and 24 h p.i., there was a decrease in the tested parameter; from 24 to 60 h p.i., the tested parameter showed only an upward trend. Analyzing the values of LZM activity in rabbits infected with *L. europaeus*/GI.1a/Rossi, it was observed that the parameter values, in relation to the control group, showed only a downward trend throughout the duration of the experiment. The obtained result was subjected to statistical analysis, which showed that in the case of infection with *L. europaeus*/GI.1a/Erfurt, all changes within the value of the examined parameter are not statistically significant, while in the group of rabbits infected with *L. europaeus*/GI.1a/Rossi, there were five statistically significant changes in the form of decreases in the activity of LZM. These changes were observed at 8, 24, 36, 48, and 60 p.i. (Figure 4). It needs to be added that similar to LZM concentration, LZM activity is characterized by noticeably high standard deviations, which was probably caused by individual animals with relatively higher or lower results. This may lead to an overall conclusion on LZM that it is not enough stable to serve as a prognostic factor.

#### *3.4. Clinical Signs Infection of Lagovirus Europaeus*/*GI.1a.*

During the experiment, many clinical symptoms were noted, including sneezing, thickening of the blood, no response to external stimuli, rapid breathing, general body stiffness, lethargy, increased thirst, increased blood clotting, and lack of appetite. All observed clinical symptoms were characteristic of the course of RHD. The first symptoms in animals were observed as early as 8 h p.i. In the following hours of the experiment, symptoms began to worsen, and continued until the animal died or until the end of the experiment. Using the Kaplan–Meier method, we calculated the probability of survival of animals from both groups of infected rabbits and plotted them in the graph (Figure 5). In the group of animals infected with *L. europaeus*/GI.1a/Erfurt, there were seven deaths from 10 tested animals during the experiment.

**Figure 5.** Percentage of rabbit survival recorded for the two tested strains (*L. europaeus*/GI.1a/Erfurt and *L. europaeus*/GI.1a/Rossi), analyzed by Kaplan–Meier method.

The first death was recorded between 12 and 24 h p.i., then between 24 and 36 h p.i. another four deaths were recorded. The last fall of the animals was observed between 56 and 60 h after the start of the experiment. The probability of survival 60 h after infection with *L. europaeus*/GI.1a/Erfurt was 30% (Figure 5).

In the case of the group of animals infected with *L. europaeus*/GI.1a/Rossi, there were four deaths from 10 infected animals. The first two deaths, as in the case of *L. europaeus*/GI.1a/Erfurt, were recorded between 12 and 24 from the beginning of the experiment. Another two deaths were observed between 24 and 36 h p.i. The probability of survival for 60 h after infection in the group of animals infected with *L. europaeus*/GI.1a/Rossi was 60% (Figure 5).

#### **4. Discussion**

Results show that MPO seems to be a sensitive marker of the condition of immune system activity of infected rabbits, as it by definition increases over the course of infection. It was previously described that MPO in viral infections in increasing, and the elevated level of MPO results from the possible lysis of leukocytes, being the first line of antimicrobial defense [34,35]. After entering of the virus into the host cells, phagocytosis is performed, and the leukocytes lead to increased reactive oxygen species through MPO and NADPH oxidase, meaning that the MPO level and activity is directly correlated with the activation of phagocytes [36]. Moreover, the level of MPO enzyme is also discussed as being a prompt indicator of endothelial dysfunction, inflammation, atherosclerosis, and oxidative stress [37]. Our results show that MPO activity levels seems to be stable and behave similarly in both studied

strains of *Lagovirus europaeus*/GI.1a levels of MPO—in the Erfurt strain levels were between 1.15 and 2.13, and in case of the Rossi strain the range was 1.07–2.34. The trend is similarly increasing from 8 to 60 h p.i. Also, no impact was noted and considered as significant as far as sex of the animals is concerned. Taking all of the above into consideration, it may be stated that MPO activity may be a more useful prognostic tool, compared to LZM, for defining the status of the immune system in the breed animals, in order to avoid the unexpected losses of the farmed rabbits.

Among animal models, the role of myeloperoxidase has been emphasized mainly in mice, where experiments with MPO-deficient mice showed their susceptibility to *Candida albicans* [38] and *Klebsiella pneumoniae* [39] infection. Also, in rats, infection causes decreases in MPO levels [40]. Keeping in mind that animal models are not ideally reflecting the human condition, it is worth noticing that in humans, several types of tissue injuries and pathogenesis of many diseases, like rheumatoid arthritis, cardiovascular diseases, liver disease, diabetes, and cancer, have been linked with MPO-derived oxidants [41]. Due to the above, it can be stated that elevated levels of MPO activity may be one of the best diagnostic tools of inflammatory and oxidative stress biomarkers. Finally, latest research shows that MPO may not only be pivotal for innate immunity, but seems to have a potential role in adaptive immunity [42] and autoimmune diseases [43].

Lysozyme is an enzyme hydrolyzing glycosidic linkages in bacterial peptidoglycan (PGN); due to the low toxicity of lysozyme, it is used as a natural preservative to control bacteria in meat products [44]. What is more, this enzyme may act as chitinase and activate bacterial autolysins, and it also exhibits antiviral activity against several human and animal viruses [45]. Also, since the natural substrate for lysozyme is PGN, it has similar functions as PGN recognition molecules, such as CD14, TLR2, or NOD-like proteins [46]. The antibacterial, antiviral, and antifungal roles of LZM has been claimed mainly in human-orientated studies [47]; its role in animal models is modestly represented. Yet, an extremely interesting study that corresponds to ours was found related to influenza, in which it was registered that the level of LZM is inhibited by this virus, whereas simultaneously, the virus does not influence the level of MPO [48]. So far, apart from ours, this is the only noted correlation between those potential biomarking substances of the immune system's condition.

The mortality of the infected animals showed that the probability of survival differs in the studied viruses of *L. europaeus*/GI.1a; however, it is important to add that this condition was probably affected by several other factors impacting immune system status.

#### **5. Conclusions**

The results on the roles of MPO and LZM as important hallmark and prognostic factors of survival of rabbits suffering from infection with *L. europaeus*/GI.1a suggests that those antimicrobial enzymes located in neutrophils not only play a significant role in defending the host from the infection, by activating the phagocytes to actively fight with the virus, but also may serve as a prognostic marker of immune system status. Basing on the results, MPO may be a more reliable indicator of inflammatory response than LZM. Further, more extensive studies on the subject are required, especially to check the correlation between the enzymes' activity and viral loads, in order to better understand the mechanism and use it in the protection of rabbits. Also, considering the fact that currently, *Lagovirus europaeus*/GII is widespread all around the world, it would be interesting to check if the role of MPO and LZM is repeatable.

**Author Contributions:** Conceptualization, P.N.-R. and R.H.; methodology, P.N.-R and R.H.; software, R.H.; validation, R.H. and D.B.; formal analysis, R.H., D.B., and P.N.-R.; investigation, R.H. and D.B.; resources, P.N.-R.; data curation, R.H.; writing—original draft preparation, R.H., D.B., and P.N.-R.; writing—review and editing, R.H., D.B., and P.N.-R.; visualization, R.H.; supervision, P.N.-R.; project administration, P.N.-R.; funding acquisition, P.N.-R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Animals* Editorial Office E-mail: animals@mdpi.com www.mdpi.com/journal/animals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-4498-4