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

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].

Open Access
*Correspondence: Huanmin.Zhang@usda.gov † Lei Zhang and Chen Zhu contributed equally to this work 1 Avian Disease and Oncology Laboratory, USDA-ARS, East Lansing, MI 48823, USA Full list of author information is available at the end of the article 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 6 3 and 7 2 , maintained at USDA, Agricultural Research Service, Avian Disease and Oncology Laboratory (East Lansing, Michigan, USA), share a common B*2 haplotype, but line 6 3 is relatively resistant and line 7 2 is highly susceptible to MD [3]. Data from vaccination and challenge trials conducted in recent years showed that line 6 3 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 7 2 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.

Experimental animals
Day-old chicks were randomly sampled from two highly inbred lines of White Leghorns, known as the line 6 3 and line 7 2 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 6 3 is resistant and line 7 2 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 6 3 and line 7 2 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 6 3 and 7 2 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 "known-MiR.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 miR-NAs 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 miR-NAs [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 6 3 HVT, line 6 3 CVI988/Rispens, line 7 2 HVT, and line 7 2 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 log 2 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].

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 6 3 and 7 2 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 6 3 birds, and 490 and 630 in bursae of the line 7 2 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 6 3 and 7 2 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.

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

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 miR-NAs 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.

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 7 2 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.

Target-gene-set enrichment in GO terms and pathways for differentially expressed microRNAs in response to HVT or CVI988/Riepsnes vaccination of line 6 3 and 7 2 birds
Target genes of differentially expressed miRNAs in response to HVT or CVI988/Rispens in line 6 3 birds Table 2  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 7 2 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 6 3 birds were significantly enriched (Bonferroni corrected p < 0.01) in a total of 237 GO terms, which included 17 molecular function 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 6 3 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 6 3 and line 7 2 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 guanylnucleotide 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 6 3 and 7 2 , between lines 6 3 and 7 2 in response to CVI988/Rispens, and line 7 2 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-23, 42, 43]. Our recent   studies using the highly inbred lines of chickens, lines 6 3 and 7 2 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 6 3 birds convey a very high protection rate, which strikingly differs from the line 7 2 birds in response to HVT. The CVI988/Rispens delivers similar protection to the line 6 3 but a little better protection to line 7 2 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 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. 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 6 3 and 7 2 chickens. The identified novelMiRs were over twice as many as the identified known gga-miRNAs in bursae of the line 6 3 and 7 2 chickens 26 days post MD vaccine inoculation (Additional file 1).
Notable difference in the number of differentially expressed miRNAs were observed in the line 6 3 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 6 3 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 6 3 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 6 3 and 7 2 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 6 3 HVT and line 6 3 control groups (Table 1) were also among the 34 differentially expressed miRNAs of the line 6 3 and 7 2 comparison inoculated with HVT. Not only so, the log 2 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 7 2 but great protection for line 6 3 birds, it is speculated here that the additional 30 differentially expressed miRNAs identified between lines 6 3 and 7 2 post HVT inoculation might be also important in modulating HVT protective efficacy in chickens like the line 6 3 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.