Skip to main content
  • Research article
  • Open access
  • Published:

Marek’s disease vaccines-induced differential expression of known and novel microRNAs in primary lymphoid organ bursae of White Leghorn

Abstract

Marek’s disease (MD) is a contagious disease of domestic chickens caused by MD viruses. MD has been controlled primarily by vaccinations, yet sporadic outbreaks of MD take place worldwide. Commonly used MD vaccines include HVT, SB-1 and CVI988/Rispens and their efficacies are reportedly dependent of multiple factors including host genetics. Our previous studies showed protective efficacy of a MD vaccine can differ drastically from one chicken line to the next. Advanced understanding on the underlying genetic and epigenetic factors that modulate vaccine efficacy would greatly improve the strategy in design and development of more potent vaccines. Two highly inbred lines of White Leghorn were inoculated with HVT and CVI988/Rispens. Bursa samples were taken 26 days post-vaccination and subjected to small RNA sequencing analysis to profile microRNAs (miRNA). A total of 589 and 519 miRNAs was identified in one line, known as line 63, 490 and 630 miRNAs were identified in the other, known as line 72, in response to HVT or CVI988/Rispens inoculation, respectively. HVT and CVI988/Rispens induced mutually exclusive 4 and 13 differentially expressed (DE) miRNAs in line 63 birds in contrast to a non-vaccinated group of the same line. HVT failed to induce any DE miRNA and CVI988/Rispens induced a single DE miRNA in line 72 birds. Thousands of target genes for the DE miRNAs were predicted, which were enriched in a variety of gene ontology terms and pathways. This finding suggests the epigenetic factor, microRNA, is highly likely involved in modulating vaccine protective efficacy in chicken.

Introduction

Marek’s disease (MD) is a proliferative disease of domestic chickens caused by Gallid alphaherpesvirus 2 (Herpesviridae), commonly known as MD virus (MDV). MD was first described by József Marek, a Hungarian veterinarian [1]. Clinical symptoms of MD vary among different lines of chickens and are dependent of exposure to different strains of MDV [2, 3], but commonly include polyneuritis, visceral lymphoma, acute transient paralysis, immunosuppression, brain edema, and acute rash [4, 5]. The morbidity and mortality could range from 11 [6] up to 100% [3]. Although MD has been well under control most of the time in most regions of the world by wide use of MD vaccines and improved management measures since the 1970s [7,8,9], sporadic outbreaks still occur from time to time all over the world. MD continues to pose a real threat and reportedly costs more than 2 billion U.S. dollars annually to the poultry industry [9, 10].

There are three commonly used MD vaccines commercially available up to date in most parts of the world. They are CVI988/Rispens, SB-1, and HVT, which were derived from MDV-1, MDV-2, and MDV-3 (also known as serotype 1, 2, and 3 MDV), respectively [11,12,13]. The SB-1 and HVT are naturally non-oncogenic. The CVI988/Rispens is a product of serial passaging in tissue culture and is considered an attenuated strain of live MDV-1. Like SB-1 and HVT, CVI988/Rispens is infectious but incapable of inducing lymphoid organ atrophy or tumor formation [14]. CVI988/Rispens and HVT have been used individually or in different combination of bivalent or trivalent format along with SB-1 on commercial farms worldwide [14]. Most users and researchers believe CVI988/Rispens is the most protective one of all commercially available vaccines against MD in chickens [15, 16], and HVT is no longer sufficiently protective or far less protective than CVI988/Rispens in preventing tumor formation against very virulent strains of MDV induction in chickens [17,18,19,20].

Protective efficacy of a vaccine is not only determined by the potency of the vaccine itself, but also host genetics of the recipients. By using B-congenic lines of experimental chickens [21] in trials of vaccination followed by MDV challenge, it was clearly shown that chicken major histocompatibility complex (MHC) B-haplotype poses significant influence on vaccinal immunity against MD [22]. This important finding was also subsequently demonstrated to hold true in commercial chickens [23]. Based on the data from those studies, chickens of B*2/B*2, B*13/B*13, B*15/B*15 and B*21/B*21 haplotypes gain the best protection against MD by MDV-1 vaccines like CVI988/Rispens; B*5/B*5 haplotype chickens are better protected by MDV-2 vaccines like SB-1; and in those genetic lines, B*2/B*2 and B*13/B*13 chickens were not better protected against MD by none of the MDV-1, MDV-2, and MDV-3 vaccines in contrast to other haplotypes of chickens.

In addition to B-haplotypes, non-MHC host genetic components also play a significant role modulating MD vaccine efficacy against MD [24]. The two highly inbred genetic lines of chickens, lines 63 and 72, maintained at USDA, Agricultural Research Service, Avian Disease and Oncology Laboratory (East Lansing, Michigan, USA), share a common B*2 haplotype, but line 63 is relatively resistant and line 72 is highly susceptible to MD [3]. Data from vaccination and challenge trials conducted in recent years showed that line 63 chickens can convey better or equivalent protective efficacy in response to MDV-3 HVT vaccination than MDV-1 CVI988/Rispens. In a sharp contrast, the line 72 birds respond to HVT very poorly and can convey some protection in response to CVI988/Rispens [25].

As early as in the 1980s, researchers realized that new routes in development of safer and more effective vaccines could be opened upon advances in molecular biology, genetics, and epigenetics [26, 27]. New evidences in the molecular, genomic and epigenomic levels in response to vaccination should highly likely lead to better understanding on protective efficacy of vaccines against infectious diseases including virus-induced cancers and offer new insights of host immune defenses as well as interaction between host immune system and vaccines. This study was designed to profile microRNAs (miRNA) in two genetically divergent lines of chicken post MD vaccine inoculation and to explore differentially expressed miRNAs in response to vaccination by CVI988/Rispens or HVT in contrast to non-vaccinated birds.

Materials and methods

Experimental animals

Day-old chicks were randomly sampled from two highly inbred lines of White Leghorns, known as the line 63 and line 72 chickens, maintained on the Avian Disease and Oncology Laboratory farm at East Lansing, Michigan, USA. Both lines are B*2 haplotype homozygous but drastically differ in genetic resistance to MD. Line 63 is resistant and line 72 is highly susceptible [21]. The two lines of chickens also drastically differ in response to MD vaccines [28].

Challenge trial

The sampled 1-day old chicks from line 63 and line 72 were randomly divided into two groups per line. One group was inoculated with HVT intraperitoneally (IP) at a dose of 2000 plaque-forming units (PFU) each; and the other group was inoculated with CVI988/Rispens (IP) with the same dose. In addition, a control group was also included with the same sample size for each of the chicken lines under the same conditions in conjunction with a joint project simultaneously (to minimize the use of animals per experiment). The corresponding control group datasets were used in this study only to examine the differential expression of miRNAs in response to HVT and CVI988/Rispens inoculation. All chickens used in this study were housed in a BSL-2 experimental facility during the trial. Feed and water were supplied ad libitum. The chickens were observed daily throughout the entire duration of the experiment. This animal experiment was approved by USDA, Avian Disease and Oncology Laboratory Institutional Animal Care and Use Committee (IACUC). The IACUC guidelines established and approved by the ADOL IACUC (April 2005) and the Guide for the care and use of Laboratory Animals by Institute for Laboratory Animal Research (2011) were closely followed throughout the experiment.

Total RNA extraction

Three chickens from each group were randomly euthanized at 26 days post-inoculation. Bursa samples were individually collected, and immediately placed into RNAlater solution (Qiagen, Valencia, CA, USA). The collected samples were stored at −20 °C until extractions of the total RNA samples. Total RNA samples were extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions.

Small RNA sequencing

Total RNA samples were quantitatively and qualitatively checked with a NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), respectively. Good quality RNA samples were chosen to construct standard cDNA libraries using Illumina TruSeq Small RNA Library Preparation kits following the manufacturer’s recommendations. Completed libraries were subjected to routine quality control (QC) checks and quantified using a combination of Qubit dsDNA HS and Agilent Bioanalyzer High Sensitivity DNA assays. The libraries were sequenced on an Illumina HiSeq 4000 sequencer using SBS (Sequencing by Synthesis) reagents. Base calling was accomplished by use of Illuima Real Time Analysis (RTA) v2.7.7 and the output of RTA was demultiplexed and then converted to FastQ format data with Ilunima Bcl2fastq v2.19.1 (Illumina, San Diego, CA, USA). The small RNA sequencing datasets supporting the results and conclusions of this article are available in the NCBI SRA repository [29]. The sequence datasets (accession numbers: SAMN11674924 to SAMN11674929) of the unvaccinated line 63 and 72 control groups used in the comparison analyses of miRNA differential expression were also deposited to NCBI SRA repository. The small RNA sequencing operations, including library preparation and preliminary reads quality control, were performed at the Research Technology Support Facility, Michigan State University (East Lansing, MI, USA).

Small RNA_Seq reads data analyses

Small RNA_Seq reads data files that passed QC were analyzed one at a time with the miRDeep* software v3.8 [30] using the default parameters except the adapter sequence and the chicken genome build index files. The adapter sequence used in the analysis is TGG AAT TCT CGG GTG CCA AGG AAC TCC AGT CAC (Illumina); and the chicken genome build index (build_bwt_idx) files were constructed based on the chromosome information of the galGal 5.0 genome build. In addition, the “knownMiR.gff” file used in miRDeep* analysis of this study was the “gga.gff3” file at the miRbase download website [31, 32], which was constructed in accordance to galGal 5.0 assembly. Target genes of differentially expressed miRNAs were predicted using the built-in target gene prediction function in miRDeep*, which employees the most commonly used target gene prediction tool, TargetScan, in predicting the target genes of known and novel miRNAs [30, 33].

Droplet digital™ PCR analysis

To validate the miRNA expressions derived from small RNA reads data, identified miRNAs were randomly elected to subject to droplet digital PCR (ddPCR) analysis to collect the absolute quantification reads of each miRNA from each of the total RNA samples of the four treatment groups (line 63 HVT, line 63 CVI988/Rispens, line 72 HVT, and line 72 CVI988/Rispens). Primers were designed for each of the selected miRNAs based on its mature sequence as described by Balcells et al. [34], which were used in the ddPCR (QX200™ ddPCR system; Bio-Rad Laboratories, Inc., Hercules, CA, USA) analysis. The cDNA samples used in ddPCR validation were reversely transcribed from individual RNA samples using the iScript™ RT Supermix Kit (Cat No. 170-8841) and following the manufacturer’s instructions (Bio-Rad). A ddPCR reaction of 25 μL in final volume was initially prepared per miRNA per biological sample containing 2 μL of cDNA, 12.5 μL of EvaGreen Supermix (Cat No. 1864034), 0.5 μL of each forward and reverse primers (200 nM; synthesized by Eurofins Genomics, Huntsville, AL), and 9.5 μL of nuclease-free water. Of which, 20 μL were loaded into one of 8 sample channels of a DG8™ cartridge (Cat No. 1864008, Bio-Rad). Each oil well was loaded with 70 μL of droplet generating oil (Cat No. 1864006, Bio-Rad). The loaded DG8™ cartridges were placed on a QX200™ droplet generator (Bio-Rad) to generate the digital droplets. Forty μL of the generated droplet emulsion for each sample were transferred to a well in a 96-well PCR plate followed by polymerase chain reaction with EvaGreen on a C1000™ Thermal Cycler (Bio-Rad). The cycling conditions were 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s, 58 °C for 60 s, and a final extension step of 98 °C for 10 min. The droplets post PCR were read well by well on a QX200™ droplet reader (Bio-Rad). PCR-positive and PCR-negative droplets in each of the wells were counted and analyzed with the QuantaSoft™ Software (Version 1.7, Bio-Rad).

Differentially expressed miRNA identification and GO terms enrichment analysis

The number of reads per miRNA for each biological sample were counted using HTSeq [35]. In each of the pairwise comparisons (between MD vaccine inoculated group and the control group for each chicken line, between the two MD vaccinated groups for each chicken line, and between the two chicken lines for each vaccinated group), differentially expressed miRNAs were identified by use of a custom R script encompassing the DESeq R package (2.1.0). A filter criterion of FDR < 0.05 and FC > 2 was enforced. For some of the differentially expressed miRNAs that ended up with a zero-statistic estimate for a normalized average TPM (baseMeanA or baseMeanB) in a contrast, an arbitrary small value of 1 was assigned to substitute the zero in order to computer a numeric fold change, and then a log2 fold change value for easier comprehension of the estimates. To better understand the functional involvements of the identified miRNAs differentially expressed in response to the vaccination, predicted target gene lists of differentially expressed miRNAs for each of the contrasts were subjected to GO terms and pathway analysis using the g:Proflier [36] online tools with the following options: Organism: Gallus gallus; Statistical domain scope: All known genes; Significance threshold: Bonferroni correction; User threshold: 0.01 [37].

Results

Small RNA sequencing

Small RNA sequencing generated an average of 35.7 million PF reads (the number of clusters that passed Illumina’s “Chastity filter”) per biological sample, with a range of 23 to 50.8 million PF reads for all 12 bursal samples of the lines 63 and 72 chickens inoculated with either HVT or CVI988/Rispens. The raw sequence datasets (accession numbers: SAMN11675491 - SMAN11675502) are available at the Sequence Read Archive (SRA) under the National Centre for Biotechnology Information (NCBI) website with an assigned BioProject SRA accession number: PRJNA543526 [38].

MicroRNA profiles

A total of 693 miRNAs was identified in this study, which were mapped to 31 pairs of chromosomes (gga chr1-28, chr32-33, chrW and Z), and was identified in no fewer than three biosamples with reads counts of five and above. Four hundred and ninety-five out of the 693 miRNAs were novel miRNAs (see Additional file 1 for a complete list of the identified miRNAs), which have been deposited to miRbase Registry for processing in validation and finalizing name assignment. The total numbers of identified miRNAs were 589 and 519 in bursae of the line 63 birds, and 490 and 630 in bursae of the line 72 birds inoculated with HVT or CVI988/Rispens, respectively (Table 1). The numbers of exclusively identified miRNAs and the numbers of identified miRNAs in common between vaccine treatment groups of the line 63 and 72 birds are depicted in a Venn diagram (Figure 1). The details of the identified miRNAs, including miRNA name (miRNA_ID), average normalized transcripts per million (Mean (TPM)), mature loci, chromosome (chr), hairpin and mature sequences, are given in Additional files 2, 3, 4, 5.

Table 1 Numbers of microRNAs identified in bursae of line 63and 72chickens 26-days post HVT or CVI988/Rispens inoculation after hatch
Figure 1
figure 1

A venn diagram depicting numbers of identified miRNAs. The numbers of identified miRNAs exclusively or commonly expressed in one of or between the vaccinated treatment groups of the line 63 and 72 birds are graphically illustrated.

Differentially expressed microRNAs in response to HVT vaccination

Four out of the 589 miRNAs identified in bursae of the line 63 birds were differentially expressed in response to HVT vaccination in contrast to an unvaccinated group of line 63 birds (1.33E−6 < p < 2.58E−4, 9.62E−4 < FDR < 4.68E−2). Two of the miRNAs were significantly upregulated (log2 Fold Change: 3.27 and 7.63), and the other two were significantly downregulated (log2 FC: −8.68 and −6.09). None of the identified miRNAs in bursae of the line 72 birds was differentially expressed in response to HVT vaccination in contrast to its counterpart of an unvaccinated group (Table 2).

Table 2 HVT and CVI988/Rispens-induced differentially expressed microRNAs in bursae within each line (63and 72) of chickens 26-days post HVT or CVI988/Rispens inoculation after hatch

Differentially expressed microRNAs in response to CVI988/Rispens vaccination

One of the novel miRNAs, novelMiR_215, was significantly downregulated (p = 2.14E−10, FDR = 1.72E−08) with a log2 FC = −4.69, and 12 miRNAs were significantly upregulated 1.60E−134 < p < 7.65E−4, 1.16E−131 < FDR < 4.27E−2) with a range of log2 FC from 3.06 to 10.77 in bursae of line 63 birds in response to CVI988/Rispens vaccination. In bursae of line 72 birds, CVI988/Rispens induced a single novel miRNA, novelMiR_1251, with significantly downregulated expression (log2 FC = −10.01, p < 4.53E−5, FDR = 3.43E−2) in contrast to line 72 unvaccinated group (Table 2).

Differentially expressed microRNAs in response to CVI988/Rispens in contrast to HVT vaccination within each of the chicken lines

Twenty-eight identified miRNAs were differentially expressed in bursae of the line 63 birds in response to CVI988/Rispens inoculation in contrast to HVT inoculation. Four of the miRNAs were significantly downregulated with log2 FC ranged from −7.67 to −3.75 (4.76E−13 < p < 1.32E−5; 3.25E−11 < FDR < 5.22E−4). The other 24 miRNAs were significantly upregulated. The log2 FC ranged from 1.51 up to 11.26 (3.05E−192 < p < 1.73E−3; 2.29E−189 < FDR < 4.49E−2). In contrast to HVT inoculation, CVI988/Rispens-induced a single novel microRNA, novelMiR_1251, with significantly downregulated expression (log2 FC = −11.45; p < 1.80E−7 and FDR < 1.37E−4) and another novel microRNA, novelMiR_215, with significantly upregulated expression (log2 FC = 3.21; p < 1.70E−4 and FDR < 4.29E−2) in line 72 birds (Table 2). The numbers of differentially expressed miRNAs in response to HVT or CVI988/Rispens inoculation and between the CVI988/Rispens and HVT vaccination within each of the two chicken lines are graphically illustrated in a Venn diagram (Figure 2).

Figure 2
figure 2

A venn diagram graphically illustrating the numbers of differentially expressed miRNAs. The numbers of differentially expressed miRNAs exclusively in one or commonly between vaccine treatment groups of a chicken line and between the lines as well as the comparisons between vaccine treatments within each of the lines are depicted.

Differentially expressed microRNAs in response to HVT or CVI988/Rispens vaccination between line 63 and line 72 birds

There were 34 identified miRNAs (8 known-miRNAs and 26 novel miRNAs) differentially expressed in bursae between the line 63 and 72 birds 26 days post HVT inoculation. Eight of those miRNAs were significantly upregulated (log2 FC: 3.49 to 7.64; 3.66E−6 < p < 5.57E−4; 1.37E−4 < FDR < 1.31E−2) and 26 were downregulated (log2 FC: −11.45 to −1.24; 6.58E−47 < p < 1.89E−3; 5.09E−44 < FDR < 4.06E−2) in line 63 birds in contrast to line 72 in response to HVT inoculation. Thirty identified miRNAs (6 known-miRNAs and 24 novel miRNAs) were differentially expressed between line 63 and 72 birds post CVI988/Rispens inoculation. Twenty-four were significantly upregulated (log2 FC: 1.45 to 11.29; 5.67E−101 < p < 1.97E−3; 4.43E−98 < FDR < 4.80E−2), and the other 6 were downregulated (log2 FC: −6.71 to  −1.69; 4.31E−21 < p < 1.59E−3; 1.12E−18 < FDR < 4.01E−2) in response to CVI988/Rispens inoculation in line 63 birds in contrast to line 72 (Table 3).

Table 3 Differentially expressed microRNAs in bursae between the line 63and the line 72birds 26-days post HVT or CVI988/Rispens inoculation after hatch

ddPCR validation of miRNA expression determined by small RNA_Seq

To validate the miRNA expression determined by small RNA_Seq analysis, ten of the identified miRNAs were selected to subject to absolute quantification of the expression by ddPCR. The selected miRNAs and the designed ddPCR primers are listed in Table 4. The summary statistics including correlation coefficients between the normalized small RNA_Seq reads and the ddPCR absolute quantification counts are given in Table 5. The correlation coefficients ranged from r = 0.63 (gga-miR-140) up to r = 0.99 (both gga-miR-31 and gga-miR-499) with a single p value smaller than 0.0001 (Table 5). Manhattan bivariate plots depicting the relationship between normalized expression of small RNA_Seq data (TPM) and the ddPCR absolute quantification reads for four of the selected miRNAs are given in Figure 3, which visually illustrates the validation. Taking all together, Table 5 and Figure 4 provided experimental evidence that highly positively supports the microRNA expression estimates derived from the small RNA_Seq data of this study.

Table 4 Primers used in ddPCR to validate the expressions of a small subgroup of identified miRNAs, which were determined by small RNA_Seq analysis
Table 5 Summary statistics of ddPCR absolute counts and RNA_Seq reads for validation of miRNA expression data detected by small RNA sequencing
Figure 3
figure 3

Bivariate plots illustrating the relationship between small RNA_Seq data and droplet digital PCR data. The normalized expression of small RNA_Seq data, transcripts per million, (TPM), and droplet digital PCR (ddPCR) data were analyzed by fitting a bivariate model. The bivariate plots for four of the validated miRNAs are given graphically showing the two sets of expression data. The correlation coefficients between the small RNA_Seq and the ddPCR data for these four miRNAs ranged from r = 0.94 (gga-miR-193b) to r = 0.99 (both gga-miR-31 and gga-miR-499), which provided highly positive support to the estimates of microRNA expression derived from the small RNA_Seq data.

Figure 4
figure 4

Manhattan plots illustrating GO term enrichments of target genes. The target genes of differentially expressed microRNAs of the line 63 birds were analyzed by g:Profiler and the enrichment in GO terms (MF: molecular function; BP: biological process; CC: cellular component) and KEGG pathways across Reactome pathways (REAC), WiKi-Pathways (WP), transcription factor (TF), and microRNA target base (MIRNA) were graphically depicted in two plots. The plots depicted the enrichments of target genes for differentially expressed microRNAs of the line 63 birds in response to HVT (top) and CVI988/Rispens (bottom) vaccination, respectively. A clearly visible difference in enrichments between the two treatment groups (HVT and CVI988/Rispens vaccination) is also demonstrated.

Predicted target genes of differentially expressed miRNAs

Varied numbers (12 up to 6153) of target genes were predicted for the differentially expressed miRNAs per contrast for all contrast groups except one, the line 72 HVT treatment group over its counterpart of control group, which resulted in no differentially expressed miRNA. The numbers of differentially expressed miRNAs, predicted target genes, involved gene ontology (GO) terms and KEGG pathways are given in Table 6 by contrast group.

Table 6 Identified numbers of differentially expressed miRNAs, predicted target genes, enriched in number of GO terms and KEGG pathways per contrast group

Target-gene-set enrichment in GO terms and pathways for differentially expressed microRNAs in response to HVT or CVI988/Riepsnes vaccination of line 63 and 72 birds

Target genes of differentially expressed miRNAs in response to HVT or CVI988/Rispens in line 63 birds were highly enriched in a variety of GO terms and pathways. In contrast, HVT failed to induce any differentially expressed miRNA and CVI988/Rispens induced a single differentially expressed miRNA in line 72 birds, which reportedly targets 12 functional genes. Therefore, there was only a very limited number of GO terms and pathways in which those target genes were enriched.

The target genes of HVT-induced differentially expressed miRNAs in the line 63 birds were significantly enriched (Bonferroni corrected p < 0.01) in a total of 237 GO terms, which included 17 molecular function (GO:MF), 172 biological process (GO:BP), and 48 cellular component (GO:CC) terms (Additional file 6; Figure 4), in addition to reactomes (REAC), WiKiPathways (WP), and transcription factors (TF). Most of the MF terms present binding functions, including protein binding, nucleic acid binding, organic cyclic compound binding, heterocyclic compound binding, DNA binding, RNA binding, enzyme binding, ion binding, and sequence-specific DNA binding. The BP terms encompassed cellular process, regulation of cellular and biological processes, regulation of gene expression, cell differentiation, cellular response to stimulus, signaling, signal transduction, and regulation of signal transduction. The CC terms are involved in membrane functionalities including membrane-bounded organelle, intracellular membrane-bounded organelle, membrane-enclosed lumen, membrane part, plasma membrane, intrinsic component of membrane, integral component of membrane, and plasma membrane bounded cell projection.

The target genes of CVI988/Rispens-induced differentially expressed miRNAs in the line 63 birds were significantly over-enriched (Bonferroni corrected p < 0.01) in a total of 1057 GO terms, which included 118 GO:MF, 811 GO:BP, and 128 GO:CC terms (Figure 4). MF terms are heavily presented in binding functions including protein binding, ion binding, organic cyclic and heterocyclic compound binding, enzyme binding, anion binding, DNA binding, transcription regulatory region DNA binding, cation binding, small molecule binding, kinase binding, and ATP binding. MF terms are also broadly involved in biological activities including phosphoric ester hydrolase activity, transcription coregulator activity, channel activity, phosphatase activity, protein tyrosine phosphatase activity, ubiquitin protein ligase activity, enzyme activator activity, molecular transducer activity, and signaling receptor activity. BP terms are presented in cell adhesion, communication, cell cycle, cycle phase transition, cycle process, cell death, development, growth, and cell junction assembly and organization, cell migration and motility, cell–cell adhesion, signaling, and signaling by Wnt. The CC terms are involved in centrosome, chromatin, cytoplasm, dendritic tree, endomembrane system, intracellular organelle, organelle lumen, and intracellular vesicle. KEGG pathways including MAPK signaling pathway, Wnt signaling pathway, Hedgehog signaling pathway, and mTOR signaling pathway (Additional file 7).

Target genes of differentially expressed miRNAs between line 63 and line 72 birds in response to HVT vaccination were also heavily enriched in (221) MF terms, (1528) BP terms, and (250) CC terms. MF terms included variety of binding functions, such as ankyrin binding, ATP binding, ATPase binding, GTP binding, cytokine receptor binding, mRNA binding, myosin binding, and nuclear hormone receptor binding. MF terms also cover pyrophosphatase activity, Ras guanyl-nucleotide exchange factor activity, signaling activity, and signaling receptor activity. The BP terms involved in a variety of systems and functions including angiogenesis, apoptotic process, autophagy, biological adhesion, regulation and process, cell differentiation and division, cell-substrate adhesion and junction assembly, and cellular response to varied compounds and stimulus. The CC terms are involved in functions including cell surface, coated membrane, vesicle, vesical membrane, transcription factor complex, and whole membrane. In addition, the target genes were also enriched in KEGG pathways including adipocytokine signaling pathway, calcium signaling pathway, C-type lectin receptor signaling pathway, ErbB signaling pathway, FoxO signaling pathway, mRNA surveillance pathway, mTOR signaling pathway, and TFG-beta signaling pathway (Additional file 8). Additional GO terms and pathways for the target genes from comparisons between CVI988/Rispens and HVT in lines 63 and 72, between lines 63 and 72 in response to CVI988/Rispens, and line 72 birds in response to CVI988/Rispens are given in rest of the Additional files 9, 10, 11, 12.

Discussion

Marek’s disease has been well under controlled since the 1970s, which, in large, is attributable to the wide use of MD vaccines in poultry flocks [9]. Commonly used commercial MD vaccines include the first-generation anti-virus-induced tumor vaccine HVT, then the second comer SB-1, and the current gold-standard MD vaccine CVI988/Rispens [39, 40]. While most researchers and industry professionals, if not all, fully recognize the great good that MD vaccines have done for the poultry industry, few, if any, claim how the MD vaccines protect chickens against the MDV-induced tumors is thoroughly understood, including immunologists [41]. The reality that this paradox persistently remains bars the advancement of knowledge-based new vaccine design and development.

Genetic mechanism underlying how MD vaccine inserts protection against MD tumor formation is also far from fully understood. Earlier studies demonstrated that MHC plays an important role in modulating MD vaccine protective efficacy in chicken [21,22,23, 42, 43]. Our recent studies using the highly inbred lines of chickens, lines 63 and 72 carrying the same MHC B*2 haplotype, showed that non-MHC genomic contents also play a significant role in enabling chicken’s ability to convey protection against MDV-induced tumor formation in response to MD vaccines. Our experimental data clearly showed the line 63 birds convey a very high protection rate, which strikingly differs from the line 72 birds in response to HVT. The CVI988/Rispens delivers similar protection to the line 63 but a little better protection to line 72 than HVT [25, 28].

Epigenetics has been demonstrated to play an important role in orchestrating key biological processes, which are implemented through a number of factors including DNA methylation, histone modification, and ncRNAs [44, 45]. The epigenetic factors, in turn, are continuously modified throughout life in response to environmental exposures [45, 46]. Epigenetics refers to the study of changes in gene expression that occur without a change in DNA sequence. DNA methylation involves the addition of a methyl group to the carbon-5 of the cytosine pyrimidine ring and typically occurs at CpG sites containing cytosine-guanine nucleotides in a linear sequence. CpG islands, short stretches of DNA with rich CpG sites, are often found at promoters of mammalian genes. DNA methylation at these sites is highly correlated with transcription status of corresponding genes. Histone modification defines discrete chromatin regions with distinct structures associated with distinct transcription states of genes. miRNAs are small-non-coding RNAs. Mature miRNAs are ~23 nucleotides in length and regulate gene expression by preventing the translation of specific mRNAs. Epigenetic mechanisms are multifaceted and complex, which provide an additional layer of transcriptional control and a layer of post-transcriptional control to regulate gene expression, and, consequently, gene function [47,48,49,50,51,52]. A very recent study using small RNA_Seq identified 24 cellular miRNAs with altered expression in response to HVT, MDV, or both inoculations. The report claims cellular miRNAs are of critical players both in protection against and mediating progression of MD [53].

Up to date, there are 2114 unique gga-miRNAs that have been coined with accession numbers (miRbase release 22.1: October 2018). This study identified a total of 693 miRNAs including 198 reported gga-miRNAs and 495 novel miRNAs in the line 63 and 72 chickens. The identified novelMiRs were over twice as many as the identified known gga-miRNAs in bursae of the line 63 and 72 chickens 26 days post MD vaccine inoculation (Additional file 1).

Notable difference in the number of differentially expressed miRNAs were observed in the line 63 birds. Vaccine and challenge trials repeatedly showed that HVT provides equally well or even better protection against MD in response to very virulent plus MDV challenge [25, 28], yet only four differentially expressed miRNAs were identified in line 63 in response to HVT inoculation, but 13 differentially expressed miRNAs were identified in response to CVI988/Rispens in this line. Further, a direct comparison between CVI988/Rispens and HVT inoculation of line 63 birds resulted in a total of 28 differentially expressed miRNAs. All these differentially expressed miRNAs (Table 1) might in different ways at different levels through their target genes as well as the GO terms and pathways insert influence to mediate each vaccine’s protective efficacy against MD.

In contrast, direct comparisons between the lines 63 and 72 birds 26 days post-vaccination resulted in a total of 34 and 30 differentially expressed miRNAs in response to HVT and CVI988/Rispens inoculation, respectively (Table 3). It is interesting to note that the four differentially expressed miRNAs identified between line 63 HVT and line 63 control groups (Table 1) were also among the 34 differentially expressed miRNAs of the line 63 and 72 comparison inoculated with HVT. Not only so, the log2 FC positive and negative signs of the four differentially miRNAs were also consistent, and the magnitude of fold change were similar. Knowing that HVT does not induce much protection to line 72 but great protection for line 63 birds, it is speculated here that the additional 30 differentially expressed miRNAs identified between lines 63 and 72 post HVT inoculation might be also important in modulating HVT protective efficacy in chickens like the line 63 birds.

In summary, a relatively large number of differentially expressed microRNAs were identified in two highly inbred lines of chicken post MD vaccine inoculation. These two genetic lines of chickens differ in the capacity to convey protective efficacy against MDV-induced tumor formation in response to MD vaccination based on previous studies. Therefore, this finding may serve as the first piece of experimental evidence that the epigenetic factor, microRNA, is highly likely involved in modulating vaccine protective efficacy in chicken through their target genes in combination with varied biological processes and pathways.

References

  1. Marek J (1907) Multiple Nervenentzündung (Polyneuritis) bei Hühnern. Dtsch Tierärztl Wochenschr 15:417–421

    Google Scholar 

  2. Gimeno IM, Witter RL, Reed WM (1999) Four distinct neurologic syndromes in Marek’s disease: effect of viral strain and pathotype. Avian Dis 43:721–737

    Article  CAS  PubMed  Google Scholar 

  3. Xie QM, Chang S, Dong KZ, Dunn JR, Song JZ, Zhang HM (2017) Genomic variation between genetic lines of White Leghorns differed in resistance to Marek’s disease. J Clin Epigenet 3:7

    Article  Google Scholar 

  4. Niu S, Jahejo AR, Jia FJ, Li X, Ning GB, Zhang D, Ma HL, Hao WF, Gao WW, Zhao YJ, Gao SM, Li GL, Li JH, Yan F, Gao RK, Bi YH, Han LX, Gao GF, Tian WX (2018) Transcripts of antibacterial peptides in chicken erythrocytes infected with Marek’s disease virus. BMC Vet Res 14:363

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Calnek BW (2001) Pathogenesis of Marek’s disease virus infection. Curr Top Microbiol Immunol 255:25–55

    CAS  PubMed  Google Scholar 

  6. Brochu NM, Guerin MT, Varga C, Lillie BN, Brash ML, Susta L (2019) A two-year prospective study of small poultry flocks in Ontario, Canada, part 2: causes of morbidity and mortality. J Vet Diagn Invest 31:336–342

    Article  PubMed  PubMed Central  Google Scholar 

  7. Gimeno IM (2008) Marek’s disease vaccines: a solution for today but a worry for tomorrow? Vaccine 26(Suppl 3):C31–41

    Article  CAS  PubMed  Google Scholar 

  8. Rozins C, Day T, Greenhalgh S (2019) Managing Marek’s disease in the egg industry. Epidemics 27:52–58

    Article  PubMed  Google Scholar 

  9. Biggs PM, Nair V (2012) The long view: 40 years of Marek’s disease research and Avian Pathology. Avian Pathol 41:3–9

    Article  PubMed  Google Scholar 

  10. Smith J, Sadeyen JR, Paton IR, Hocking PM, Salmon N, Fife M, Nair V, Burt DW, Kaiser P (2011) Systems analysis of immune responses in Marek’s disease virus-infected chickens identifies a gene involved in susceptibility and highlights a possible novel pathogenicity mechanism. J Virol 85:11146–11158

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Schat KA, Calnek BW (1978) Protection against Marek’s disease-derived tumor transplants by the nononcogenic SB-1 strain of Marek’s disease virus. Infect Immun 22:225–232

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Rispens BH, van Volten H, Mastenbroek N, Maas JL, Schat KA (1972) Control of Marek’s disease in the Netherlands. II. Field trials on vaccination with an avirulent strain (CVI 988) of Marek’s disease virus. Avian Dis 16:126–138

    Article  CAS  PubMed  Google Scholar 

  13. Witter RL, Offenbecker L (1978) Duration of vaccinal immunity against Marek’s disease. Avian Dis 22:396–408

    Article  CAS  PubMed  Google Scholar 

  14. Baigent SJ, Smith LP, Nair VK, Currie RJ (2006) Vaccinal control of Marek’s disease: current challenges, and future strategies to maximize protection. Vet Immunol Immunopathol 112:78–86

    Article  CAS  PubMed  Google Scholar 

  15. Wozniakowski G, Niczyporuk JS (2015) Detection of specific UL49 sequences of Marek’s disease virus CVI988/Rispens strain using loop-mediated isothermal amplification. J Virol Methods 221:22–28

    Article  CAS  PubMed  Google Scholar 

  16. Baigent SJ, Nair VK, Le Galludec H (2016) Real-time PCR for differential quantification of CVI988 vaccine virus and virulent strains of Marek’s disease virus. J Virol Methods 233:23–36

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Liu S, Sun W, Huang X, Zhang W, Jia C, Luo J, Shen Y, El-Ashram S, He C (2017) A promising recombinant herpesvirus of turkeys vaccine expressing PmpD-N of Chlamydia psittaci based on elongation factor-1 alpha promoter. Front Vet Sci 4:221

    Article  PubMed  PubMed Central  Google Scholar 

  18. McKimm-Breschkin JL, Faragher JT, Withell J, Forsyth WM (1990) Isolation of very virulent Marek’s disease viruses from vaccinated chickens in Australia. Aust Vet J 67:205–209

    Article  CAS  PubMed  Google Scholar 

  19. Sharma JM, Witter RL (1983) Embryo vaccination against Marek’s disease with serotypes 1, 2 and 3 vaccines administered singly or in combination. Avian Dis 27:453–463

    Article  CAS  PubMed  Google Scholar 

  20. Spencer JL, Robertson A (1975) Influence of vaccination with avirulent herpesvirus on subsequent infection of chickens with virulent Marek’s disease herpesvirus. Am J Vet Res 36:1235–1239

    CAS  PubMed  Google Scholar 

  21. Bacon LD, Hunt HD, Cheng HH (2000) A review of the development of chicken lines to resolve genes determining resistance to diseases. Poult Sci 79:1082–1093

    Article  CAS  PubMed  Google Scholar 

  22. Bacon LD, Witter RL (1994) Serotype specificity of B-haplotype influence on the relative efficacy of Marek’s disease vaccines. Avian Dis 38:65–71

    Article  CAS  PubMed  Google Scholar 

  23. Bacon LD, Witter RL (1994) B haplotype influence on the relative efficacy of Marek’s disease vaccines in commercial chickens. Poult Sci 73:481–487

    Article  CAS  PubMed  Google Scholar 

  24. Witter RL (2001) Protective efficacy of Marek’s disease vaccines. Curr Top Microbiol Immunol 255:57–90

    CAS  PubMed  Google Scholar 

  25. Chang S, Xie Q, Dunn JR, Ernst CW, Song J, Zhang HM (2014) Host genetic resistance to Marek’s disease sustains protective efficacy of herpesvirus of turkey in both experimental and commercial lines of chickens. Vaccine 32:1820–1827

    Article  CAS  PubMed  Google Scholar 

  26. Tomasi TB, Magner WJ, Khan AN (2006) Epigenetic regulation of immune escape genes in cancer. Cancer Immunol Immunother 55:1159–1184

    Article  PubMed  Google Scholar 

  27. Dougan G, Highfield P (1985) Molecular biology, the new genetics and vaccine development. Med Lab Sci 42:393–398

    CAS  PubMed  Google Scholar 

  28. Chang S, Dunn JR, Heidari M, Lee LF, Song J, Ernst CW, Ding Z, Bacon LD, Zhang H (2010) Genetics and vaccine efficacy: host genetic variation affecting Marek’s disease vaccine efficacy in White Leghorn chickens. Poult Sci 89:2083–2091

    Article  CAS  PubMed  Google Scholar 

  29. Heidari M, Zhang HM (2019) PRJNA543524. NCBI. https://www.ncbi.nlm.nih.gov/sra

  30. An J, Lai J, Lehman ML, Nelson CC (2013) miRDeep*: an integrated application tool for miRNA identification from RNA sequencing data. Nucleic Acids Res 41:727–737

    Article  CAS  PubMed  Google Scholar 

  31. Kozomara A, Birgaoanu M, Griffiths-Jones S (2019) miRBase: from microRNA sequences to function. Nucleic Acids Res 47:D155–D162

    Article  CAS  PubMed  Google Scholar 

  32. Kozomara A (2019) miRBase. http://www.mirbase.org/ftp.shtml

  33. Lewis BP, Burge CB, Bartel DP (2005) Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell 120:15–20

    Article  CAS  PubMed  Google Scholar 

  34. Balcells I, Cirera S, Busk PK (2011) Specific and sensitive quantitative RT-PCR of miRNAs with DNA primers. BMC Biotechnol 11:70

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Anders S, Pyl PT, Huber W (2015) HTSeq–a Python framework to work with high-throughput sequencing data. Bioinformatics 31:166–169

    Article  CAS  PubMed  Google Scholar 

  36. g:Profiler (2019) ELIXIR. http://biit.cs.ut.ee/gprofiler/index.cgi

  37. Reimand J, Kull M, Peterson H, Hansen J, Vilo J (2007) g:Profiler–a web-based toolset for functional profiling of gene lists from large-scale experiments. Nucleic Acids Res 35:W193–200

    Article  PubMed  PubMed Central  Google Scholar 

  38. Zhang HM, Heidari M (2019) PRJNA543526. NCBI. https://www.ncbi.nlm.nih.gov/sra/PRJNA543526

  39. Garcia M (2017) Current and future vaccines and vaccination strategies against infectious laryngotracheitis (ILT) respiratory disease of poultry. Vet Microbiol 206:157–162

    Article  CAS  PubMed  Google Scholar 

  40. Reddy SM, Izumiya Y, Lupiani B (2017) Marek’s disease vaccines: current status, and strategies for improvement and development of vector vaccines. Vet Microbiol 206:113–120

    Article  CAS  PubMed  Google Scholar 

  41. Boodhoo N, Gurung A, Sharif S, Behboudi S (2016) Marek’s disease in chickens: a review with focus on immunology. Vet Res 47:119

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Bacon LD, Witter RL (1992) Influence of turkey herpesvirus vaccination on the B-haplotype effect on Marek’s disease resistance in 15.B-congenic chickens. Avian Dis 36:378–385

    Article  CAS  PubMed  Google Scholar 

  43. Bacon LD, Witter RL (1993) Influence of B-haplotype on the relative efficacy of Marek’s disease vaccines of different serotypes. Avian Dis 37:53–59

    Article  CAS  PubMed  Google Scholar 

  44. Aslani S, Jafari N, Javan MR, Karami J, Ahmadi M, Jafarnejad M (2017) Epigenetic modifications and therapy in multiple sclerosis. Neuromol Med 19:11–23

    Article  CAS  Google Scholar 

  45. van Otterdijk SD, Michels KB (2016) Transgenerational epigenetic inheritance in mammals: how good is the evidence? FASEB J 30:2457–2465

    Article  PubMed  CAS  Google Scholar 

  46. Marsit CJ (2015) Influence of environmental exposure on human epigenetic regulation. J Exp Biol 218:71–79

    Article  PubMed  PubMed Central  Google Scholar 

  47. Hamilton JP (2011) Epigenetics: principles and practice. Dig Dis 29:130–135

    Article  PubMed  PubMed Central  Google Scholar 

  48. Ling C, Groop L (2009) Epigenetics: a molecular link between environmental factors and type 2 diabetes. Diabetes 58:2718–2725

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Reddy MA, Natarajan R (2011) Epigenetic mechanisms in diabetic vascular complications. Cardiovasc Res 90:421–429

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. van Weerd JH, Koshiba-Takeuchi K, Kwon C, Takeuchi JK (2011) Epigenetic factors and cardiac development. Cardiovasc Res 91:203–211

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Wulff BE, Nishikura K (2012) Modulation of microRNA expression and function by ADARs. Curr Top Microbiol Immunol 353:91–109

    CAS  PubMed  Google Scholar 

  52. Yan C, Boyd DD (2006) Histone H3 acetylation and H3 K4 methylation define distinct chromatin regions permissive for transgene expression. Mol Cell Biol 26:6357–6371

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Hicks JA, Liu HC (2019) Impact of HVT vaccination on splenic miRNA expression in Marek’s disease virus infections. Genes (Basel) 10:E115

    Article  CAS  Google Scholar 

Download references

Acknowledgements

Dr Jiyuan An is gratefully appreciated and acknowledged for the construction of chicken genome build index files used in miRDeep* analysis of the chicken small RNA_Seq data. The authors also thank Bernice Li and Yanfen Zhai for their excellent technical support on conducting the animal experiment, tissue sampling, ddPCR analysis, and GO term enrichment work; and Melanie Flesberg, on total RNA extraction, qualification, and quantification.

Funding

This work was solely supported by Congressional Prepared Budget for Research ministered by US department of Agriculture and Agricultural Research Service. Lei Zhang was financially supported by the China Scholarship Council (No. 201803250035) sponsoring visiting scholars studying abroad. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture.

Author information

Authors and Affiliations

Authors

Contributions

LZ and CZ conducted bioinformatic analyses of the small RNA sequence dataset, prepared the draft of the manuscript, tabulated and depicted the results. MH revised the manuscript. KD prepared the R customer scripts and revised the manuscript. SC and QX participated in the animal experimental design, helped conduct the animal experiment, and revised the manuscript. HZ conceived and supervised this study. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Huanmin Zhang.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Complete list of unique microRNAs identified in all biosamples of this study.

Unique microRNAs including known microRNAs and novel microRNAs identified in two lines of chickens in response to vaccination.

Additional file 2: miRNAs identified in bursae of Line 6

3chickens 26 days post HVT inoculation. A microRNA profile of line 63 in response to HVT vaccination.

Additional file 3: miRNAs identified in bursae of Line 6

3chickens 26 days post Rispens inoculation. A microRNA profile of line 63 in response to CVI988/Rispens vaccination.

Additional file 4: miRNAs identified in bursae of Line 7

2chickens 26 days post HVT inoculation. A microRNA profile of line 72 in response to HVT vaccination.

Additional file 5: miRNAs identified in bursae of Line 7

2chickens 26 days post Rispens inoculation. A microRNA profile of line 72 in response to CVI988/Rispens vaccination.

Additional file 6: GO terms enriched by target genes of differentially expressed miRNAs in line 6

3bursae 26 days post HVT inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs of the line 63 birds in response to HVT vaccination.

Additional file 7: GO terms enriched by target genes of differentially expressed miRNAs in line 6

3bursae 26 days post Rispens inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs of the line 63 birds in response to CVI988/Rispens vaccination.

Additional file 8: GO terms enriched by target genes of differentially expressed miRNAs in bursae between line 6

3and line 72birds 26 days post HVT inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs between line 63 and line 72 birds 26 days post HVT vaccination.

Additional file 9: GO terms enriched by target genes of differentially expressed miRNAs in bursae between Rispens and HVT vaccinated line 6

3birds 26 days post-inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs between Rispens and HVT inoculated line 63 birds 26 days post-vaccination.

Additional file 10: GO terms enriched by target genes of differentially expressed miRNAs in bursae between Rispens and HVT vaccinated line 7

2birds 26 days post-inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs between Rispens and HVT vaccinated line 72 birds 26 days post-vaccination.

Additional file 11: GO terms enriched by target genes of differentially expressed miRNAs in bursae between line 6

3and line 72birds 26 days post Rispens inoculation. Gene Ontology analysis of a list of target genes of differentially expressed microRNAs between line 63 and line 72 birds 26 days post CVI988/Rispens vaccination.

Additional file 12: GO terms enriched by target genes of a differentially expressed miRNA in bursae of line 7

2birds in response to Rispens vaccination. Gene Ontology analysis of a list of target genes of a differentially expressed microRNA of line 72 birds 26 days post CVI988/Rispens vaccination.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, L., Zhu, C., Heidari, M. et al. Marek’s disease vaccines-induced differential expression of known and novel microRNAs in primary lymphoid organ bursae of White Leghorn. Vet Res 51, 19 (2020). https://doi.org/10.1186/s13567-020-00746-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13567-020-00746-4