Epigenetic Factor MicroRNAs Likely Mediate Vaccine Protection Efficacy against Lymphomas in Response to Tumor Virus Infection in Chickens through Target Gene Involved Signaling Pathways

Simple Summary Epigenetic factor microRNAs may modulate vaccine protective efficacy well, supported by experimental data showing that the expression of some microRNAs drastically differs between well-protected and poorly protected genetic lines of chickens. Marek’s disease vaccination and virus challenge trials were conducted in two genetically divergent inbred lines of chickens, and hundreds of microRNAs were identified. However, a small number of miRNAs were observed with significant differential expression per line per vaccination treatment group. One of the vaccines was HVT, which reportedly protects the two lines of birds with a significant difference. The target genes of the differentially expressed miRNAs in response to HVT were reportedly involved with many Gene Ontology terms and pathways, which suggests there is a high complexity of the genetic/epigenetic mechanism modulating vaccine protective efficacy. Abstract Epigenetic factors, including microRNAs (miRNAs), play an important role in affecting gene expression and, therefore, are involved in various biological processes including immunity protection against tumors. Marek’s disease (MD) is a highly contagious disease of chickens caused by the MD virus (MDV). MD has been primarily controlled by vaccinations. MD vaccine efficacy might, in part, be dependent on modulations of a complex set of factors including host epigenetic factors. This study was designed to identify differentially expressed miRNAs in the primary lymphoid organ, bursae of Fabricius, in response to MD vaccination followed by MDV challenge in two genetically divergent inbred lines of White Leghorns. Small RNA sequencing and bioinformatic analyses of the small RNA sequence reads identified hundreds of miRNAs among all the treatment groups. A small portion of the identified miRNAs was differentially expressed within each of the four treatment groups, which were HVT or CVI988/Rispens vaccinated line 63-resistant birds and line 72-susceptible birds. A direct comparison between the resistant line 63 and susceptible line 72 groups vaccinated with HVT followed by MDV challenge identified five differentially expressed miRNAs. Gene Ontology analysis of the target genes of those five miRNAs revealed that those target genes, in addition to various GO terms, are involved in multiple signaling pathways including MAPK, TGF-β, ErbB, and EGFR1 signaling pathways. The general functions of those pathways reportedly play important roles in oncogenesis, anti-cancer immunity, cancer cell migration, and metastatic progression. Therefore, it is highly likely that those miRNAs may, in part, influence vaccine protection through the pathways.


Introduction
Epigenetics is within the stable of heritable phenotypic changes in organisms induced through modification of gene expression and intervention of gene translation other than alteration of the structure of DNA, the genetic code itself [1,2].Over the last two decades, clear and definite evidence from numerous studies showed that epigenetics is deeply involved in modulating host immunity and tumorigenesis processes upon infections [3][4][5][6].The epigenetic factor microRNA (miRNA) is reportedly involved in various bioprocesses, including oncogene functions and tumorigenesis suppression [7].Host cellular miRNAs, like gga-miR-26a, gga-miR-181a, and gga-miR-130a, for instance, play a role in lymphoma cell proliferation suppression [8,9].
Marek's disease (MD) is a contagious disease of domestic chickens caused by an avian α-herpesvirus [10,11], commonly known as MD virus (MDV).MD has been controlled primarily by the use of MD vaccines since 1970 [10], yet sporadic MD outbreaks take place worldwide [12,13].Therefore, the commercial companies are clearly aware that MD remains a threat to the poultry industry, as it continues to impose an annual cost of over USD 2 billion to the industry [14].Commonly used commercial MD vaccines include Herpesvirus of turkeys (HVT), SB-1, and CVI988/Rispens, in addition to others that have been under active tests and development [15].MD vaccine efficacy is dependent upon multiple factors including host genetics [16].Our previous studies showed the protective efficacy of a MD vaccine can differ drastically from one genetic line of birds to the next [17].Advancement in understanding the underlying genetic and epigenetic factors that modulate vaccine protective efficacy against tumor incidence would greatly improve the strategy for design and development of new and more potent vaccines, and therefore, would empower better control of infectious diseases, like Marek's disease, in the chicken.
Epigenetic factors include histone modification, DNA methylation, and non-coding RNAs [49].Examination of histone modifications and differential chromatin marks in chickens in response to MDV challenge revealed significant differences between inbred lines of chickens in their resistance to MD at both cytolytic and latency phases [44,50,51].In a whole genome histone modification study, Luo et al. found that trimethylations of histone H3K4me3 and H3K27me3 marks were positively and negatively correlated with expression of protein coding genes, respectively, both in MDV-challenged and MDVunchallenged chickens [50].These reports strongly suggest that histone modifications are highly likely to be involved in MD and potentially associated with tumorigenesis [52] in chicken.DNA methylation affects DNA transcription, consequently passing on impacts on health and disease status in addition to other phenotypic characteristics [38,53].In a separate study, Luo et al. found that the promoter DNA methylation level of the CD4 gene was downregulated in response to MDV infection, which was negatively correlated with CD4 gene expression in an MD-susceptible line of chickens [42].The promoter region of CD30, a key gene likely associated with tumorigenesis in MD, was hypomethylated in response to MDV infection and in MD lymphoma [46,54].DNA methylation level alteration in response to MDV infection was also detected between a genetic line of chickens highly susceptible to MD and another line relatively resistant to MD [43].
MicroRNA (miRNA) is an important class of short non-coding RNAs.Mature miRNAs are approximately 21-25 nucleotides in length.miRNA negatively regulates gene expression by binding to the 3 ′ -UTR or 5 ′ -UTR of mRNAs to inhibit translation or by initiating mRNA degradation to prevent translation of target genes [55].Dysregulated expression of microRNAs has been observed in numerous disease states, particularly in human cancers, neurologic disorders, autoimmune diseases, metabolic diseases, cardiovascular diseases, and stress-induced emotional and suicidal behavior disorders [56,57].It is reported that gga-miR-21 was observed to have significantly upregulated expression in response to MDV infection, and miR-26a expression was consistently downregulated in MD tumors.gga-miR-21 and miR-26a are reportedly known to facilitate an MDV oncogene, Meq, in lymphomagenesis and to suppress MD tumor formation, respectively [9,58].
Two highly inbred lines of chickens, known as line 6 3 and line 7 2 , were developed at the USDA, Agricultural Research Service facility, Avian Disease and Oncology Laboratory (ADOL, East Lansing, MI, USA).These lines have subsequently been maintained at both the ADOL and the US National Poultry Research Center (Athens, Georgia) facilities.The 6 3 and 7 2 lines share a common major histocompatibility complex haplotype (B*2) but differ in susceptibility to avian leukosis virus infection and in resistance to MD [19,59].Our recent studies revealed that MDV infection induced differential expression of nineteen and nine miRNAs in the two genetic lines of birds, respectively [60], and the MD vaccine HVT and CVI988/Rispens induced differential expression of four and thirteen exclusive microRNAs in the 6 3 line birds, and 0 and one miRNA in the 7 2 line birds, respectively [61].To advance the fundamental understanding of microRNAs and microRNA expression in association with MDV infection post-vaccination, this study was designed to profile miRNAs and to identify miRNAs with dysregulated expression in the primary lymphoid organ, bursa of Fabricius, in the two genetically divergent inbred lines of chickens in response to a very virulent plus MDV (vv+MDV) challenge post-vaccination.Predicted target genes of the differentially expressed miRNAs were subjected to Gene Ontology (GO) terms analysis to elucidate the likely biological significance of the differentially expressed miRNAs, potentially in association with MD vaccine protective efficacy.

Experimental Animals and a Vaccination-Challenge Trial
Chicks on the day of hatching were sampled from two highly inbred lines of White Leghorns, known as line 6 3 and line 7 2 , maintained in the Avian Disease and Oncology Laboratory at East Lansing, Michigan, U.S.A.Both lines are B*2 haplotype homozygous, but the lines of chickens drastically differ in genetic resistance to MD. Line 6 3 is resistant and line 7 2 is highly susceptible [59].The two lines of chickens also differ in response to MD vaccines [17].
The sampled chicks from line 6 3 and line 7 2 were randomly divided into two treatment groups per line.One group was inoculated with HVT at a dose of 2000 plaqueforming units (PFU) each on the day of hatching, and the other group was inoculated with CVI988/Rispens at the same dosage on the same day.Both vaccinated groups were challenged with a very virulent plus strain of Marek's disease virus (vv+MDV), known as 648A, at a dose of 500 PFU each on the fifth day post-hatch, intraperitoneally (IP).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 simultaneously conducted joint project.The corresponding control group datasets (SRA accession numbers SAMN11674924-SAMN11674929 under the PRJNA543524 BioProject (URL https: //www.ncbi.nlm.nih.gov/sra/PRJNA543524;accessed 2 January 2024) were used in this study only for computation of the differential expression of miRNAs in response to HVT and CVI988/Rispens vaccination.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 the Institute for Laboratory Animal Research (2011) were closely followed throughout the experiment.

Extraction of Total RNA Samples
Three chickens from each treatment group were randomly euthanized at 26 days postvaccine 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 Illumina Real Time Analysis (RTA) v2.7.7, and the output of RTA was demultiplexed and then converted to FastQ-format data with Illumina Bcl2fastq v2.19.1 (Illumina, San Diego, CA, USA).The small RNA sequence datasets analyzed and reported in this study were deposited to SRA under NCBI (URL https://www.ncbi.nlm.nih.gov/sra;accessed 1 February 2024).The sequence datasets (accession numbers SAMN11674924 to SAMN11674929) of the unvaccinated lines 6 3 and 7 2 control groups used in the comparison analyses were generated simultaneously for a joint project following exactly the same protocols as described above, and were also deposited to SRA under NCBI.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).

Data Analyses of Small RNA_Seq Reads
The data files of small RNA_Seq reads of all samples that passed QC were subjected to analysis one at a time with the miRDeep* software v3.8 [62], with 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 latest "gga.gff3"file at the miRbase download website (URL http://www.mirbase.org/download/;accessed 21 February 2024), which was constructed in accordance with galGal 5.0 assembly.Target genes of differentially expressed miRNAs were predicted using the built-in target gene prediction function in miRDeep*.

Validation of the Small RNA Sequence Reads Data by Droplet Digital PCR
A small number of random samples from the identified miRNAs were subjected to Droplet Digital PCR (ddPCR) analysis to validate the levels of small RNA sequence reads datasets.Primers were customarily designed for the analysis for each of the sampled miRNAs using the mature miRNA sequences following the procedures described by Balcells et al., 2011 [63].The cDNA samples used in the ddPCR validation were reverse transcribed from the individual total RNA samples using the iScript™ RT Supermix Kit (Cat No. 170-8841) following the manufacturer's instruction (Bio-Rad (Tokyo, Japan)).A ddPCR reaction of a 25 µL 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 primer (200 nM; synthesized by Eurofins Genomics, Huntsville, AL, USA), and 9.5 µL of nuclease-free water.Of which, 20 µL was loaded into one of eight 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.A volume of 40 µL of the generated droplet emulsion per sample was 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).The primers designed for ddPCR validation of the selected miRNAs are given in Table 1.

Identification of Differentially Expressed miRNAs and GO Terms Enrichment
The number of reads per microRNA for each biological sample were counted using HTSeq [64].In each of the pairwise comparisons (between the MD-vaccine-inoculated group and the control group within each chicken line, between the two MD-vaccinated groups within each chicken line, and between the two chicken lines within each vaccinated group), differentially expressed microRNAs were identified by use of a custom R script encompassing the DESeq R package (2.1.0).Filter criteria of FDR < 0.05 and FC > 2 were enforced.For some of the differentially expressed miRNAs that ended up with a zerostatistic estimate for a normalized average of transcripts per million (TPM) in a contrast, an arbitrary small value of 5 was assigned to substitute for the zero to compute a numeric fold change, and then Log 2 fold change values were computed for easier comprehension of the estimates.To better understand the functional involvements of the identified microR-NAs differentially expressed in response to the vaccination, predicted target gene lists of differentially expressed microRNAs for each of the contrasts were subjected to GO terms and pathway analysis using the g:Proflier (URL http://biit.cs.ut.ee/gprofiler/index.cgi; accessed 2 February 2024) online tools with the following options: Organism: Gallus gallus; Statistical domain scope: All known genes; Significance threshold: Bonferroni correction; User threshold: 0.01 [65].

Small RNA Sequence Reads, Reads Quality, and Deposition to NCBI
Small RNA sequencing generated an average of 35.7 million pass-filter (PF) reads (the number of clusters that passed Illumina's "Chastity filter") per biological sample.The PF reads ranged from 23 to 50.8 million with Phred Quality Scores above 30 for the 12 bursal samples of the line 6 3 and 7 2 chickens vaccinated with either HVT or CVI988/Rispens and followed by vv+MDV (648A) challenge.The numbers of PF reads of each biological sample for both chicken lines and treatment groups are listed in Table 2.The raw sequence datasets of this study are available at the Sequence Read Archive (SRA) under National Centre for Biotechnology Information (PRJNA544273-SRA-NCBI (nih.gov)).

Validation of the Small RNA Sequence Reads Data by ddPCR
A randomly selected small set of miRNAs was subjected to ddPCR analysis using custom-designed primers (Table 1).The absolute quantification of the four randomly selected miRNAs in the same sets of total RNA samples by ddPCR was analyzed against the corresponding normalized PF reads (TPM) of the small RNA sequence reads data with a bivariate model.The correlation coefficients between the two sets of expression data for the examined miRNAs ranged from r = 0.87 to r = 0.98, with a statistic of p < 0.001 (Figure 1).This provided highly positive support in validation of the small RNA sequence reads data generated in this study and the estimates of the miRNA expression derived from the small RNA sequence datasets.

MicroRNA Profiles of the Chicken Line by Vaccine Treatment Groups
A total of 630 unique miRNAs were identified in bursal samples of the two lines of chickens among all the treatment groups.Close to one-third of the identified miRNAs were known chicken miRNAs and the rest were novel miRNAs (Table S1).A total of 481 and 524 known and novel miRNAs, respectively, were identified in samples of the line 63 chickens vaccinated with HVT and CVI988/Rispens followed by MDV challenge.A total

MicroRNA Profiles of the Chicken Line by Vaccine Treatment Groups
A total of 630 unique miRNAs were identified in bursal samples of the two lines of chickens among all the treatment groups.Close to one-third of the identified miRNAs were known chicken miRNAs and the rest were novel miRNAs (Table S1).A total of 481 and 524 known and novel miRNAs, respectively, were identified in samples of the line 6 3 chickens vaccinated with HVT and CVI988/Rispens followed by MDV challenge.A total of 457 and 510 known and novel miRNAs, respectively, were identified in samples of line 7 2 chickens subjected to the same vaccination and MDV challenge.The total numbers of known and novel miRNAs identified in the birds of the lines 6 3 and 7 2 in response to HVT or CVI988/Rispens vaccination followed by MDV challenge are summarized in Table 3 for each line by treatment group.The detailed information, including miRNA name (miR_ID), hairpin loci, mature loci, precursor sequence (Sequence), mature sequence (Mature.miR),and normalized number of reads (TPM) is given in Tables S2-S5.A heatmap was constructed, which offers a partial scope of visual illustration of expressional differences and similarities of the identified miRNAs observed among the four line-by-vaccine treatment groups: line 6 3 vaccinated with HVT or CVI988/Rispens followed by MDV challenge, and line 7 2 vaccinated with HVT or CVI988/Rispens followed by MDV challenge (Figure 2).As shown, this subset of identified miRNAs most differed in expression between the line 6 3 and line 7 2 chickens subjected to the CVI988/Rispens vaccination followed by MDV challenge.The HVT-vaccinated and MDV-challenged line 6 3 and line 7 2 groups did differ from each other, but both differed more from the CVI988/Rispens-vaccinated line 7 2 group than the CVI988/Rispens-vaccinated line 6 3 group followed by MDV challenge.A Venn diagram depicts the number of miRNAs identified uniquely within each of the four line-by-vaccine treatment groups, or in common between two treatment groups, or among three or all four treatment groups (Figure 3).The identified numbers of unique miRNAs ranged from 10 to 28 among the four treatment groups, with a dominantly large number of 346 miRNAs in common among all four treatment groups.The HVT-vaccinated MDV-challenged line 6 3 and line 7 2 groups were observed with 5 miRNAs in common, the smallest number of identified miRNAs between any two groups of the line-by-vaccine treatment groups; the largest number was 34 miRNAs in common, observed both between the line 6 3 HVT and CVI988/Rispens-vaccinated MDV-challenged groups and between the CVI988/Rispens-vaccinated MDV-challenged line 6 3 and line 7 2 groups.

Differentially Expressed miRNAs between Treatment Groups
Both HVT and CVI988/Rispens vaccination followed by the vv+MDV (648A) challenge induced differentially expressed miRNAs in contrast to control groups of the line 63 and line 72 birds.The control groups of birds were hatched on the same day as the vaccine and MDV treatment groups of birds, and were housed and reared under the same conditions; however, neither were vaccinated nor MDV challenged.The majority of the differentially expressed miRNAs were novel, along with a few of the known chicken miRNAs (Table 4).

Differentially Expressed miRNAs between Treatment Groups
Both HVT and CVI988/Rispens vaccination followed by the vv+MDV (648A) challenge induced differentially expressed miRNAs in contrast to control groups of the line 6 3 and line 7 2 birds.The control groups of birds were hatched on the same day as the vaccine and MDV treatment groups of birds, and were housed and reared under the same conditions; however, neither were vaccinated nor MDV challenged.The majority of the differentially expressed miRNAs were novel, along with a few of the known chicken miRNAs (Table 4).HVT vaccination followed by the vv+MDV challenge induced seven and six differentially expressed miRNAs in bursae of the line 6 3 and 7 2 chickens, respectively, 21 days post-MDV inoculation.Five out of the seven differentially expressed miRNAs in line 6 3 birds in response to HVT vaccination and MDV challenge were significantly upregulated (p < 6.5 × 10 −5 ), with a range of Log 2 fold change (FC) from 3.04 to 12.76.The other two miRNAs were significantly downregulated (p < 3.1 × 10 −5 ), with Log 2 FC of −6.96 and −8.07, respectively.All the differentially expressed miRNAs in the line 7 2 birds were significantly upregulated (p < 2.5 × 10 −4 ) except one, a novel miRNA (novelMiR_692) with a Log 2 FC of −11.61.The upregulated miRNAs in the line 7 2 birds were observed with Log 2 FC ranging from 5.87 to 12.47 (Table 4).

Differentially Expressed miRNAs between CVI988/Rispens and HVT Vaccination Groups Followed by vv+MDV Challenge within Each Line of Chickens
By comparison between the treatment groups of CVI988/Rispens and HVT vaccinations followed by MDV challenge, only a few miRNAs were identified within each of the lines of birds that were differentially expressed.The novelMiR_508_1, novelMiR_508_2, and novelMiR_508_3 (identical in their mature miRNA sequence but different at the hairpin loci and mature loci; see Table S1 for details) were significantly upregulated in expression (p < 1.2 × 10 −4 ) in response to CVI988/Rispens vaccination and MDV challenge with a Log 2 FC of 6.32, in contrast to HVT vaccination followed by MDV challenge in the line 6 3 birds.By the same comparison, three miRNAs (gga-mir-30a*, gga-mir-205b, and novel-MiR_215) were significantly downregulated in expression (p < 2.5 × 10 −6 ) in line 7 2 birds.The Log 2 FC ranged from −6.15 to −11.10 (Table 5).Five and twenty-two differentially expressed miRNAs were identified between line 6 3 and line 7 2 birds that were vaccinated with HVT and CVI988/Rispens, respectively, followed by MDV challenge.
Ten of the twenty-two differentially expressed miRNAs were observed with significantly higher expression (p < 1.1 × 10 −3 ), and the other 12 miRNAs had significantly lower expression (p < 1.7 × 10 −3 ), in line 6 3 birds than the line 7 2 birds vaccinated with CVI988/Rispens followed by MDV challenge.The 10 miRNAs having significantly higher expression in line 6 3 birds were observed with Log 2 FCs ranging from 1.16 to 14.06.The 12 miRNAs having significantly lower expression in line 6 3 birds were observed with Log 2 FCs ranging from −1.11 to −8.90 (Table 6).

Gene Ontology (GO) Terms Enrichment of Predicted Target Genes of the Differentially Expressed miRNAs
A total of 2035 and 2134 target genes were predicted for the differentially expressed miRNAs of line 6 3 and line 7 2 birds, respectively, in response to the HVT vaccination followed by MDV challenge treatment in contrast to each line's control group.These target genes were highly enriched in a total of 1729 and 1739 GO terms and pathways, respectively.The target genes were highly enriched in biological process (BP) and molecular function (MF) terms in common between the two HVT-MDV treatment groups of the two chicken lines (6 3 and 7 2 ).The BP terms included the following: (a) cell and cellular processes like cell adhesion, communication, development, differentiation, migration, cell morphogenesis, motility, cell death, cell-cell and Wnt signaling, cell surface receptor signaling pathway, cellular component processes including homeostasis, assembly, biogenesis, morphogenesis, cellular localization and metabolic process, cellular response to stimulus and stress; (b) epithelial processes, like cell differentiation, tube morphogenesis, and development; (c) immune system development; (d) intracellular processes, like protein transport and signal transduction; and (e) negative and positive regulation involving cellular processes, signal transduction, signaling, transcription, apoptotic process, cell communication, gene expression, intracellular signal transduction, and protein phosphorylation.The MF terms included binding processes like ATP binding, DNA bindings, enzyme binding, identical protein binding, kinase and lipid bindings, nucleic acid and nucleotide bindings, protein and protein domain-specific bindings, signaling receptor binding, and small molecule binding.There was a total of 182 GO terms enriched exclusively with target genes of the line 6 3 HVT-MDV treatment group.Those GO terms included a variety of channel activity terms like ion and calcium channel activity, gated channel activity, and voltage-gated channel activity, and transmembrane transport activity like active transmembrane transporter activity and monovalent inorganic cation transmembrane transporter activity under the MF category.More GO terms under the BP category were observed, which included angiogenesis, regulation of response to external stimulus, Wnt signaling pathway, negative regulation of cell motility and migration, cell-cell adhesion, immune system process, and transmem-brane receptor protein serine/threonine kinase signaling pathway (Tables S6 and S7).The GO terms for each of the functional categories, which were enriched with target genes of differentially expressed miRNAs of the line 6 3 and line 7 2 HVT-vaccinated treatment group and MDV-challenged treatment group of birds, are visually depicted in the top and middle Manhattan-like plots in Figure 4, respectively.
tivity, and voltage-gated channel activity, and transmembrane transport activity like active transmembrane transporter activity and monovalent inorganic cation transmembrane transporter activity under the MF category.More GO terms under the BP category were observed, which included angiogenesis, regulation of response to external stimulus, Wnt signaling pathway, negative regulation of cell motility and migration, cell-cell adhesion, immune system process, and transmembrane receptor protein serine/threonine kinase signaling pathway (Tables S6 and S7).The GO terms for each of the functional categories, which were enriched with target genes of differentially expressed miRNAs of the line 63 and line 72 HVT-vaccinated treatment group and MDV-challenged treatment group of birds, are visually depicted in the top and middle Manhattan-like plots in Figure 4, respectively.A direct comparison between the line 6 3 and line 7 2 HVT-MDV treatment groups resulted in identification of five differentially expressed miRNAs (Table 6), which predictively target a sum of 774 genes.Those 774 target genes were enriched collectively in 907 GO terms, of which 32 out of a total of 45 MF GO terms are highly involved with binding activities ranging from ion and small molecule binding to enzyme and protein binding activities.Close to half of the BP GO terms are involved with regulation processes including regulation of cellular process, transduction by RNA polymerase II, cell communication, signaling, signal transduction, cell mortality, and cell differentiation.Some of the target genes were also enriched in key KEGG pathways, which included MAPK signaling pathway, TGF-beta signaling pathway, ErbB signaling pathway, and EGFR1 signaling pathway.The target genes of these five differentially expressed miRNAs were especially and highly more enriched in the Transfac (TF) GO terms than any other category of GO terms (Figure 4, bottom plot), which included the critical GO terms of the proto-oncogene c-Jun and c-Fos motifs, as well as the cancer-associated Smad3 motifs (Table S8).
The GO terms enrichment results for target genes of differentially expressed miRNAs for the rest of the treatment groups and comparisons are detailed in Tables S9-S13, which include treatment groups of line 6 3 and line 7 2 birds subjected to CVI988/Rispens vaccination followed by MDV challenge (Table S9 and Table S10, respectively), the comparisons between the line 6 3 CVI988/Rispens-MDV and HVT-MDV groups (Table S11), the line 7 2 CVI988/Rispens-MDV and HVT-MDV groups (Table S12), and between the line 6 3 and line 7 2 groups subjected to CVI988/Rispens-MDV groups (Table S13).

Discussion
As one of the avian tumor virus-induced diseases, Marek's disease has been well under control since the 1970s in most parts of the world, which is primarily attributable to the wide use of MD vaccines in poultry flocks wherever applicable [11].The commonly used commercial MD vaccines include the gold-standard MD vaccine, CVI988/Rispens; the very first anti-tumor vaccine, HVT; and the MDV-2 vaccine, SB-1 [66,67].While most, if not all, researchers and industry professionals fully recognize the great good that MD vaccines have done for the poultry industry, few, if any, claim to thoroughly understand the mechanism of how MD vaccines protect against MDV-induced tumorigenesis and tumor formation in chickens [68].The reality that this paradox persistently remains bars the advancement of knowledge-based new vaccine design and development, as well as better control of MD in poultry.
Protective efficacy of different vaccines varies [69][70][71] and protective efficacy of a single MD vaccine also differs between genetic lines of chickens due to host genetics [16,17].The highly inbred line 6 3 and line 7 2 chickens used in this study, sharing a common major histocompatibility complex haplotype (B*2), strikingly differ in their ability to convey protective efficacy by up to 40% in response to CVI988/Rispens and up to 82% in response to HVT vaccination against very virulent plus MDV challenge [17,72].The observed phenotypic differences in the ability to convey protective efficacy in response to vaccination against MDV-induced tumor formation may result from a rather complex mechanistic system potentially involving host genomic and epigenomic variation and vaccine by chicken line interactions [61,73].This point is well indicated, in turn, by the fact that the MD vaccine's immunologic mechanism of protection remains poorly understood [74].The understanding of genomic and epigenomic mechanisms underlying vaccine protective efficacy is even more limited.
Mounting evidence from research of various biological fields has demonstrated that epigenetic factors like DNA methylations, histone modifications, and non-coding short RNAs, particularly microRNAs, play a pivotal role in differentiation of immune cell subsets and immune functions of those immune cells, which result in alteration of immune responses in reaction to foreign antigens.It is also reported that microRNAs regulate a wide array of cellular processes including, but not limited to, cell proliferation, differentiation, and apoptosis.Dysregulated expression of miRNAs is reportedly associated with malignant transformation and tumor formation [75,76].To investigate the potential roles that miRNAs play in modulating immune response and, subsequently, tumor incidence post-MD vaccination and MDV challenge, two highly inbred lines of White Leghorns, the lines 6 3 and 7 2 , were profiled for miRNAs and compared to identify differentially expressed miRNAs in response to HVT or CVI988/Rispens vaccination followed by a vv+MDV challenge.Over six hundred miRNAs were identified among a set of four chicken line by vaccine treatment groups (line 6 3 HVT+MDV, line 6 3 CVI988/Rispens+MDV, line 7 2 HVT+MDV, and line 7 2 CVI988/Rispens+MDV).Although the total numbers of miRNAs identified within each of the line 6 3 and line 7 2 groups subjected to either HVT+MDV or CVI988/Risepens+MDV treatment were similar (Table 3), the numbers of differentially expressed miRNAs in response to CVI988/Rispens+MDV (in contrast to each line's non-vaccinated non-challenged control group) differed (14 and 9) between the line 6 3 and line 7 2 groups (Table 4).It is also noticeable that both HVT+MDV and CVI988/Rispens induced four miRNAs (gga-mir-19b*, gga-mir-425-5p, novelMiR_530, and novelMiR_547) in common with significantly upregulated expression and one miRNA (novelMiR_91) with significantly downregulated expression in common within the line 6 3 birds (Table 4).Since HVT induced equally good or even better protection of the line 6 3 birds than the CVI988/Rispens vaccination [72,73], all or some of these five miRNAs dysregulated in expression between the line 6 3 and line 7 2 birds in response to the same HVT and MDV treatment would be highly probably associated with the previously observed highly protective efficacy against MD in the line 6 3 birds [17,72].Furthermore, of the five differentially expressed miRNAs, the two miRNAs novelMiR_530 and novelMiR_91 were identified in the HVT+MDV treatment group of the line 6 3 birds but not in the same treatment group of the line 7 2 birds.Therefore, these two miRNAs might play an even more critical role in modulating the protective efficacy of HVT vaccination in birds like those of line 6 3 [17,72,73].
Emerging evidence also shows miRNAs play key roles in fine-tuning critical biological processes that affect host immune homeostasis, and consequently infectious diseases and cancer development [77][78][79][80].It is also shown that tumor viruses are capable of, on one hand, modulating the expression of cellular miRNAs to facilitate their own infection processes by invading the host immune system, and, on the other, to dysregulate oncogenes and tumor suppressor genes [81].Dysregulated expression of miRNAs in response to virus infection is reportedly associated with cancerous disease development [82][83][84] and host antiviral immunity [85].It has been shown that the line 6 3 birds are capable of conveying over 80% protection against vv+MDV-induced MD in response to HVT vaccination, while the line 7 2 birds are virtually incapable of conveying any protection, or convey a poorly minimal protection, against the same vv+MDV-induced MD in response to HVT [72].Seven and six differentially expressed miRNAs were identified within the line 6 3 and 7 2 birds, respectively, in response to HVT vaccination and vv+MDV challenge in this study (Table 4), which were primarily or totally different from the sets of differentially expressed miRNAs identified in the two lines of birds in response only to HVT vaccination or only to MDV challenge alone (Table S14).These results may have shown that the induced dysregulation of miRNA expression in response to HVT and MDV treatment is not a simple sum of the HVT treatment plus the MDV treatment.
Direct comparison between the HVT+MDV treatment groups of the line 6 3 and line 7 2 birds identified five differentially expressed miRNAs.The same comparison between the CVI988/Rispens+MDV treatment groups of line 6 3 and line 7 2 birds identified 22 differentially expressed miRNAs (Table 6).Only two of the differentially expressed miRNAs (gga-miR-1684a and novelMiR-1062) were in common between the two comparison sets of HVT+MDV and CVI988/Rispens+MDV treatment groups.
Gene Ontology analysis of the target genes of the five differentially expressed miR-NAs between the line 6 3 HVT+MDV treatment group and the line 7 2 HVT+MDV group indicated that those miRNAs, through their target genes, are likely involved with multiple signaling pathways, including the MAPK signaling pathway, which elicits many responses in cells evoked by factors including environmental changes and plays a major role in onco-genesis [86]; the TGF-β signaling pathway, which regulates development, homeostasis, and tissue repairments, and reportedly plays a major role in cancer suppression [87]; the ErbB signaling pathway, which promotes autophosphorylation and subsequent downstream signaling cascades through binding with numerous signal transducers and was confirmed as having an important role involving cancers [88][89][90]; and the EGFR1 signaling pathway, which reportedly plays a role in immune defense against pathogen infection, diverse cellular processes, including cell apoptosis, proliferation, differentiation, migration, and the cell cycle, and immune regulation in both vertebrates and invertebrates [91].Further research is warranted to clarify the details by providing additional experimental evidence as proof of the mechanism used by the vaccine to achieve the various levels of protection.The findings of this study suggest that these five microRNAs are likely to be instrumental for the host epigenetics to interpolate between vaccine protection and virus-induced disease incidence collectively through varied GO terms and pathways including signaling pathways, with which the microRNAs' target genes are involved.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/vetsci11040139/s1,Table S1: Identified unique miRNAs in bursae of two genetically divergent lines of chickens in response to Marek's disease vaccination and virus challenge; Table S2: Line 6 3 miRNA profile in response to HVT vaccination and MDV challenge.Table S3: Line 6 3 miRNA profile in response to Rispens vaccination and MDV challenge.Table S4: Line 7 2 miRNA profile in response to HVT vaccination and MDV challenge.Table S5: Line 7 2 miRNA profile in response to Rispens vaccination and MDV challenge.Table S6: GO terms enrichment of target genes predicted for line 6 3 differentially expressed miRNAs in response to HVT and MDV inoculation.Table S7: GO terms enrichment of target genes predicted for line 7 2 differentially expressed miRNAs in response to HVT and MDV inoculation.Table S8: GO terms enrichment of predicted target genes of differentially expressed miRNAs between lines 6 3 and 7 2 in response to HVT and MDV inoculation.Table S9: GO terms enrichment of target genes predicted for line 6 3 differentially expressed miRNAs in response to Rispens and MDV inoculation.Table S10: GO terms enrichment of target genes predicted for line 7 2 differentially expressed miRNAs in response to Rispens and MDV inoculation.Table S11: GO terms enrichment of predicted target genes of line 6 3 differentially expressed miRNAs between Rispens and HVT vaccination followed by MDV challenge.Table S12: GO terms enrichment of predicted target genes of line 7 2 differentially expressed miRNAs between Rispens and HVT vaccination followed by MDV challenge.Table S13: GO terms enrichment of predicted target genes of differentially expressed miRNAs between lines 6 3 and 7 2 in response to Rispens and MDV inoculation.Table S14: Differentially expressed miRNAs in response to HVT+MDV, HVT, or MDV in line 6 3 birds.

21 Figure 1 .
Figure1.A bivariate plot demonstrating the moderately high association between the droplet digital PCR-generated absolute quantifications of miRNA expressions and the normalized small RNA sequence reads data (TMP) of a set of randomly picked miRNAs out of the identified miRNAs of this study.The result suggests that the small RNA sequence reads datasets generated and analyzed in this study were statistically valid.

Figure 1 .
Figure1.A bivariate plot demonstrating the moderately high association between the droplet digital PCR-generated absolute quantifications of miRNA expressions and the normalized small RNA sequence reads data (TMP) of a set of randomly picked miRNAs out of the identified miRNAs of this study.The result suggests that the small RNA sequence reads datasets generated and analyzed in this study were statistically valid.

Figure 2 .
Figure 2. A heatmap graphically demonstrating the identified miRNA expression levels among th four treatment groups.A subsample of identified miRNAs from each of the four line-by-vaccine and MDV treatment groups-line 63 vaccinated with HVT or CVI988/Rispens (Risp), and line 72 vac cinated with HVT or Risp (all followed by Marek's disease virus challenge)-are plotted to illustrat the expressional differences (normalized counts) and similarities.The expressional levels of the iden tified miRNAs of the line 63 Risp-vaccinated and MDV-challenged treatment groups were most distan

Figure 2 .
Figure 2.A heatmap graphically demonstrating the identified miRNA expression levels among the four treatment groups.A subsample of identified miRNAs from each of the four line-by-vaccine and MDV treatment groups-line 6 3 vaccinated with HVT or CVI988/Rispens (Risp), and line 7 2 vaccinated with HVT or Risp (all followed by Marek's disease virus challenge)-are plotted to illustrate the expressional differences (normalized counts) and similarities.The expressional levels of the identified miRNAs of the line 6 3 Risp-vaccinated and MDV-challenged treatment groups were most distant from the Risp-vaccinated and MDV-challenged line 7 2 group; the line 6 3 and line 7 2 HVT-vaccinated MDV-challenged treatment groups did differ from each other, but both were more distant from the Risp-vaccinated MDV-challenged line 7 2 than the line 6 3 Risp-vaccinated MDV-challenged treatment groups.

Figure 3 .
Figure 3.A Venn diagram depicting the numbers of identified miRNAs within each, between two, and among three or all four line-by-vaccine treatment groups.As shown, the numbers of identified miRNAs unique within each of the treatment groups ranged from 10 to 28.More than half (346) of the identified miRNAs (630) were observed in common in all four line-by-vaccine treatment groups.

Figure 3 .
Figure 3.A Venn diagram depicting the numbers of identified miRNAs within each, between two, and among three or all four line-by-vaccine treatment groups.As shown, the numbers of identified miRNAs unique within each of the treatment groups ranged from 10 to 28.More than half (346) of the identified miRNAs (630) were observed in common in all four line-by-vaccine treatment groups.

Figure 4 .
Figure 4. Manhattan-like plots depicting GO terms (MF: molecular function, BP: biological process, CC: cellular component) and pathways (KEGG: Kyoto Encyclopedia of Genes and Genomes, REAC: reactome, WP: WiKiPathways, TF: Transfac, MIRNA: miRTarBase, HP: human protein ontology) of target genes of the differentially expressed miRNAs induced by a combination of HVT Marek's disease vaccine and Marek's disease virus (MDV) inoculations: (top) line 6 3 HVT plus MDV treatment group versus line 6 3 control group; (middle) line 7 2 HVT plus MDV treatment group versus line 7 2 control group (all the control groups of birds were neither vaccinated nor MDV challenged); (bottom) line 6 3 HVT plus MDV treatment group versus line 7 2 of the same treatment group.The numbers of GO terms and pathways differed little between the top and the middle plots.The bottom plot, however, showed there were significant numbers of GO terms and pathways involved with the target genes of differentially expressed miRNAs between the line 6 3 and line 7 2 birds subjected to the same vaccination and MDV challenge.

Table 1 .
Customarily designed primers used in ddPCR validation of the RNA sequencing reads data of the randomly picked miRNAs.

Table 2 .
Sequenced samples and small RNA sequence reads of the samples.

Table 3 .
Number of miRNAs identified in two highly inbred lines of White Leghorn layers in response to Marek's disease vaccination and MDV challenge.

Table 4 .
Differentially expressed miRNAs in bursae of Fabricius of line 6 3 and 7 2 chickens in response to HVT or CVI988/Rispens vaccination followed by a vv+MDV 1 challenge.

Table 6 .
Differentially expressed miRNAs between line 6 3 and line 7 2 chickens challenged with Marek's disease virus post-vaccination.