Molecular characterization of the meq oncogene of Marek’s disease virus in vaccinated Brazilian poultry farms reveals selective pressure on prevalent strains

Abstract

Marek’s disease virus (MDV) has become an increasingly virulent pathogen in the poultry industry despite vaccination efforts to control it. Brazil has experienced a significant rise of Marek’s disease (MD) outbreaks in recent years. Our study aimed to analyze the complete meq gene sequences to understand the molecular epidemiological basis of MD outbreaks in Brazilian vaccinated layer farms. We detected a high incidence rate of visceral MD (67.74%) and multiple circulating MDV strains. The most prevalent and geographically widespread genotype presented several clinical and molecular characteristics of a highly virulent strain and evolving under positive selective pressure. Phylogenetic and phylogeographic analysis revealed a closer relationship with strains from the USA and Japan. This study sheds light on the circulation of MDV strains capable of infecting vaccinated birds. We emphasize the urgency of adopting preventive measures to manage MDV outbreaks threatening the poultry farming industry.

Keywords:

• Mardivirus gallidalpha2
• meq oncoprotein
• selective pressure
• phylodynamics
• phylogeography

1. Introduction

Marek’s disease is widespread in all poultry-producing countries, and is characterized by rapid-onset lymphoid tumors, immunosuppression, and paralysis (Nair et al. 2020). It is caused by the Marek’s disease Virus (MDV), classified as Mardivirus gallidalpha2 within the subfamily Alphaherpesvirinae (ICTV 2022).

Regarding its pathogenicity, MDV is classified into four pathotypes: mild (mMDV), virulent (vMDV), very virulent (vvMDV), and very virulent plus (vv + MDV) strains (Witter et al. 2005). Since its initial identification, an increasing virulence has been reported, which can occur even in fully vaccinated birds (Osterrieder et al. 2006; Shi et al. 2020).

MDV encodes the oncogene meq, which plays a crucial role in the pathogenesis of the disease. The standard isoform of meq consists of 339 amino acids encompassing a proline-glutamine-rich domain (Pro/Gln), the basic region (BR), the leucine zipper (ZIP), and the transactivation domain (TAD) (Stolz and McCormick 2020). Additional alternative isoforms have also been identified, namely, long-meq (L-meq), short-meq (S-meq), and very short-meq (VS-meq) (Chang et al. 2002). Meq is involved in various biological processes including replication, oncogenesis, and immunosuppression in chickens. The presence of meq polymorphisms in the BR, ZIP, and the proline repeats region (PRR) in the transactivation domain has been associated with the virulence level of MDV strains (Qian et al. 1995; Chang et al. 2002; Nair and Kung 2004; Shamblin et al. 2004; Yu et al. 2013; Trimpert et al. 2017; Dunn et al. 2019; Sato et al. 2022; Song et al. 2022).

Since 1969, vaccination has been the primary strategy for preventing MD. These vaccines can elicit an anti-tumor response, thereby reducing mortality and morbidity rates. However, vaccines do not provide complete protection against viral infection, replication, and transmission among chickens (Davidson and Nair 2005). Therefore, despite vaccination efforts to reduce MDV incidence in poultry farms through these programs, sporadic outbreaks of MDV have continued to occur in several countries in recent years (Atkins et al. 2013; Mescolini et al. 2020a; Deng et al. 2021; Ozan et al. 2021; Kannaki et al. 2022).

In Brazil, vaccination against MDV is mandatory in commercial hatcheries (MAPA 2007). In layer farms, vaccination programs typically involve dual schemes including HVT + CVI988/Rispens. However, epidemiological studies of MDV in Brazil are scarce, and molecular characterization based on complete meq sequencing is currently lacking (Chacón et al. 2019; Torres et al. 2019).

The aim of this study was to explore the molecular epidemiological and characteristics of MDV outbreaks in vaccinated layer farms in Brazil. By identifying potential determinants, we aimed to explain the emergence, prevalence, and virulence characteristics of these outbreaks.

2. Materials and methods

2.1. Samples and MDV detection

This study included 116 samples from MD outbreaks, across 40 poultry layer farms in Brazil from 2018 to 2022 in Brazil (Table S1). These samples were collected from the states of São Paulo (27 farms, 75 samples), Minas Gerais (8 farms, 29 samples), Sergipe (2 farms, 8 samples), Pernambuco (1 farm, 2 samples), Espírito Santo (1 farm, 1 sample), and Paraíba (1 farm, 1 sample). The pathological findings observed were categorized into two MD forms based on the main pathological presentation: visceral (mainly neoplastic with visceral lymphoproliferative lesions) and neural (mainly non-neoplastic). These MD forms (phenotypes) were evaluated using Fisher’s exact test (significance level = 0.05) to identify differences between the Brazilian genotypes based on amino acid profiles.

Vaccination programs against MD on these farms included immunization with bivalent CVI988 (MDV1) and FC126 (HVT) vaccine via in ovo. Each sample constituted a pool collected from a specific organ/source gathered from five different birds within the same flock (Table S1).

Nucleic acid extraction was performed with the BioGene Viral kit (Quibasa Ltda., MG, Brazil) according to the manufacturer’s recommendations. MDV serotype and differentiation from the vaccinal CIV988 strain were performed through real-time PCR targeting the pp38 gene (Gimeno et al. 2014). Additionally, other oncogenic viruses, namely avian leukosis virus (ALV) and reticuloendotheliosis virus (REV), were tested by PCR as previously reported to verify the presence and completeness of their genomes amplifying the genes LTR, gag, pol, and env (García et al. 2003; Cao et al. 2013; Chacón et al. 2020, 2022). The complete meq gene of MDV was amplified as previously described (López-Osorio et al. 2017) and sequenced on a 3500xL Genetic Analyzer.

2.2. Sequence and phylogenetic analysis

Phylogenetic analysis was based on the complete meq gene sequences of 30 Brazilian strains obtained in this study and 351 sequences retrieved from GenBank. Only sequences with available data of isolation year and geographical origin were included. Sequences were aligned using MAFFT online service (Katoh et al. 2019). Best-fit substitution models for phylogenetic analysis were estimated with ModelTest-NG v0.1.7 (Darriba et al. 2020). A maximum likelihood (ML) tree was reconstructed with PhyML v3.3.2 (Guindon et al. 2010) by using a GTR + R substitution model and visualized with iTOL v6 (Letunic and Bork 2021).

The amino acid changes analysis in the meq gene was conducted by establishing a consensus sequence from all the 381 strains. Also, thirty sequenced Brazilian strains from this study were compared with the MDV pathotype representative strains.

2.3. Selective pressure analysis

The identification of positive selection in the meq gene was performed through the Codon-based Z-test of selection using the Nei-Gojobori method. To reject the hypothesis that changes in the meq gene were driven by stochastic processes, Tajima’s test of neutrality was performed using MEGA v11.0.10 (Tamura et al. 2021).

The exploration of sites within the meq gene was carried out using models M1a, M2a, M7, M8, and M8a implemented in the codeml program of PAML 4.9j (Yang 1997). The likelihood ratio test (LRT) was used to compare M1a, M7, and M8a models that assume no positive selection (null models; ω < 1) with the M2a and M8 models that assume positive selection (alternative models; ω > 1). Sites with Empirical Bayes (BEB) posterior probability ≥ 0.95 were considered to be under positive selection. This analysis included only strains with known pathotypes and non-redundant Brazilian.

Positive selection in the meq gene was detected through BUSTED (gene-wide diversifying selection) SLAC, FEL, and FUBAR (pervasive diversifying or purifying selection sites), and MEME (episodic diversifying selection sites). All methods were performed on the server Datamonkey 2.0 (Weaver et al. 2018). For this analysis, 381 meq sequences were used. Localization of selective pressure sites along the meq protein was illustrated with the IBS program (Liu et al. 2015).

2.4. Phylodynamic analysis

Temporal signal (clocklikeness) of the data was evaluated with TempEst v.1.5.3. The sampling dates were converted to decimal notation using the R package ‘lubridate,’ and clocklikeness was assessed using the Best-fitting root technique by analyzing the Correlation Coefficient (root to tip divergence vs time). Evolutionary analyses were inferred by using BEAST, v2.6.7 (Suchard et al. 2018). Redundant sequences were discarded, keeping the oldest one for each case. The model GTR + G + I was selected as the nucleotide substitution model under the Relaxed Log Normal Clock model with Coalescent Constant Population using two independent Markov Chain Monte Carlo (MCMC) runs consisted of 50 million generations with a sampling of 5000 generations each. The Maximum Clade Credibility tree (MCC) was generated with a 20% burn-in with TreeAnnotator software (Suchard et al. 2018). The results were analyzed by corroborating effective sampling sizes (ESS) over 200 and drawn by using Tracer, v1.7.2 and Figtree v1.4.4 (Rambaut et al. 2018). We also estimate the mean evolutionary rates for meq domains as follows: Pro/Gln (1–29), BR (30-84), ZIP (85–120), and TAD (121–399).

2.5. Phylogeographic analysis

The geographic distribution of detected Brazilian MDV genotypes were visualized by using the microreact web application (Argimón et al. 2016). To propose possible MDV migration pathways, an ancestral state reconstruction and identification of time to the most recent common ancestor (tMRCA) were performed using Nextstrain v4.2.0 (Hadfield et al. 2018). The pipeline was built using Snakemake v7.14.2 (Köster and Rahmann 2012), augur v18.0 (Huddleston et al. 2021) was used for tracking evolution from nucleotide alignment, iq-tree v1.6.12 (Nguyen et al. 2015), and TreeTime v0.9.4 (Sagulenko et al. 2018), for the phylogeny and timed tree inference respectively. The results were visualized using auspice v2.38.0 (Hadfield et al. 2018).

3. Results

3.1. Gross examination and MDV detection

After necropsies, the findings revealed that visceral MD accounted for 67.74% (21/31) of the farms with clinically defined cases, distributed across four different states in the Southeast and Northeast regions of Brazil. In contrast, neural MD was detected in 32.25% (10/31) of the farms and was limited to the São Paulo state. Nine farms reported no pathological signs. Among the cases of visceral MD, the most prevalent clinical manifestations were visceral lymphomas, depression, and increased mortality. Conversely, in cases of neural MD, paralysis and increased mortality were the predominant signs observed (Table S1). Furthermore, qPCR assays confirmed the presence of field MDV1 in 94.83% of the samples (110/116), while the vaccinal CIV988 strain was detected in 57.76% of the samples (67/116) from 33 of the 40 farms. The positive samples included the brain, bursa, feather, gizzard, heart, intestine, kidney, liver, lung, muscle, proventriculus, spleen, sciatic nerve, skin, thymus, and trachea (Table S1). Regarding other oncogenic viruses, REV was detected in 8.6% of the samples (10/116), infecting various organs such as feathers, thymus, liver, proventriculus, kidney, lung, spleen, and skin. In all cases of co-infection, REV was identified by presenting the complete genome (LTR, gag, pol, and env genes). However, it could not be verified whether it was integrated into the MDV genome or whether it had a synergistic effect on the pathological findings since these samples belonged to two farms that presented the neural form of MD including paralysis, hemorrhage and folliculitis as the main signs. On the other hand, the presence of ALV was not detected in any sample.

3.2. Sequence and phylogenetic analysis

Complete sequencing of the meq gene was performed in samples from 30 different outbreaks (Table S2). The phylogenetic tree based on complete meq sequences showed the presence of seven clusters (Figure 1). Brazilian MDV strains sequenced in this study were included in three different clusters. Strain USP-1879 was included in cluster 2, along with the attenuated strains (CVI988 and 814), low virulent USA and Australian strains, and three strains isolated from Anseriformes, among others. Cluster 3 included twenty-two Brazilian strains alongside USA strains GA (vMDV), RB-1B (vvMDV), Md5 (vvMDV), and W (vv + MDV). This cluster also encompasses Indian and Japanese strains including three isolated from Anseriformes. Other seven Brazilian strains were included in cluster 5 together with a vMDV strain (571), several Italian strains from poultry, and Asian strains isolated from Meleagrinidae and Anseriformes, among others.

Figure 1. Maximum-likelihood phylogenetic tree of 381 complete meq gene sequences, including 30 Brazilian strains of this study. The tree was inferred under a GTR + R substitution model with support values based on 1000 bootstrap replicates indicated on branches. Additional data of meq gene sequences were indicated as follows: phylogenetic cluster (inner ring), isolation country (Middle ring), and number of PRRs (outer bars). Brazilian strains were highlighted in light green, and strains isolated from meleagrinidae, passeriformes, and anseriformes were highlighted in light yellow, light red, and light purple respectively.

Deduced meq protein sequences of Brazilian strains were analyzed, and 16 of the most polymorphic sites distributed along the BR, ZIP and TAD domains were compared against other pathotyped MDV strains (Table 1). The most frequent polymorphisms detected in the Brazilian strains, when compared to the consensus sequence, included E77K, Y80D, A139T, P153Q, and 194delP. Notably, one Brazilian strain (USP-1879) had an L-meq size (398 aa), while the other 29 had the standard meq (338 and 339 aa). Twenty-one of them had 338 aa caused by a proline deletion 194delP and consequently the loss of a PRR. This deletion is shared across mMDV, vMDV, and vvMDV strains.

Table 1. Molecular characteristics and amino acid substitutions^A of meq oncoprotein of Brazilian strains compared with MDV reference strains.

					BR			ZIP		TAD
				amino acid site^B→	71	77	80	115	119	139	153	176	180	194	216	217	277	283	320	326
Strain/sequence	Pathotype	Size (aa)	PRRs N(N°)	Brazilian Genotype	A	E	Y	A	C	A	P	P	T	P	P	P	L	A	I	T
814	att	398	7		S	K	D	*	*	T	*	*	*	–	*	*	*	*	*	*
CVI988	att	398	7		S	*	D	V	*	T	*	*	*	–	*	*	*	*	*	I
USP-1879		398	7	IV	S	*	D	V	*	T	*	*	*	–	*	*	*	*	*	I
CU-2	mMDV	398	7		S	*	D	V	*	T	*	*	*	–	*	*	*	*	*	*
BC-1	vMDV	398	7		S	A	D	*	*	T	*	*	*	–	*	*	*	*	*	*
FT158	vMDV	398	5		S	A	A	*	*	T	*	*	*	–	*	A	*	*	*	*
MPF57	vMDV	398	5		S	A	A	*	*	T	*	*	*	–	*	A	*	*	*	*
02LAR	vMDV	398	5		S	A	D	*	*	T	*	*	*	–	*	A	*	*	*	*
JM/102W	vMDV	399	7		*	*	D	V	*	T	*	*	*	*	S	*	*	*	*	*
04CRE	vMDV	398	5		S	A	D	*	*	T	*	*	*	–	*	A	*	*	*	*
GA	vMDV	339	5		*	K	D	V	*	T	*	*	*	*	*	*	*	*	*	*
571	vMDV	339	4		*	*	*	*	*	T	*	H	*	*	*	*	*	*	*	*
617A	vMDV	339	4		*	*	*	V	R	T	*	*	*	*	*	A	*	*	*	*
Woodlans1	vvMDV	398	5		S	A	A	*	*	T	*	*	*	–	*	A	*	*	*	*
Md5	vvMDV	339	4		*	K	D	V	*	T	*	*	*	*	*	A	*	V	T	*
549	vvMDV	339	2		*	K	D	V	R	T	Q	A	A	*	*	A	*	*	*	*
595	vvMDV	339	2		*	K	D	V	R	T	Q	A	A	*	*	A	*	*	*	*
643P	vvMDV	339	2		*	K	D	V	R	T	Q	A	A	*	*	A	F	*	*	*
RB-1B	vvMDV	339	5		*	K	D	V	*	T	*	*	*	*	*	*	*	*	*	*
W	vv + MDV	339	4		*	K	D	V	*	T	*	*	*	*	*	A	*	V	T	*
New	vv + MDV	339	2		*	K	D	V	R	T	Q	A	*	*	*	A	*	V	T	*
X	vv + MDV	339	2		*	K	D	V	R	T	Q	A	A	*	*	A	*	*	*	*
648A	vv + MDV	339	2		*	K	D	V	R	T	Q	A	A	*	*	A	P	*	*	*
USP-1015		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1023		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1322		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1323		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1324		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1325		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1334		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1379		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1482		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1483		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1487		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1525		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1577		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1876		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1881		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1882		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1920		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1936		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1937		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1948		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-2449		338	3	I	*	K	D	*	*	T	L	*	*	–	*	*	*	*	*	*
USP-1328		339	4	III	*	K	D	*	*	T	L	*	*	*	*	*	*	*	*	*
USP-1284		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-1305		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-1374		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-1458		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-1480		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-1481		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*
USP-2193		339	5	II	*	*	*	*	*	T	*	*	*	*	*	*	*	*	*	*

^A Amino acid substitutions: (-) = gap. (*) = conserved. (colored) polymorphism. Brazilian strains are highlighted in bold.

^B Amino acid position according to the 339 aa‐standard meq isoform.

We classified four genotypes for Brazilian strains based on the amino acids profiles (Table 1). Genotype I included 21 Brazilian strains with 338 aa, 3 PRRs, and the 194delP deletion. Also, they presented the substitution E77K which is prevalent in vvMDV and vv + MDV strains and the novel mutation P153L, leading to the loss of a PRR. This site is also mutated in vvMDV and vv + MDV strains (595, New, and 648 A) as P153Q. Genotype II consisted of 7 strains with 339 aa and 5 PRRs. This genotype was characterized by presenting a meq aa profile similar to the consensus sequence. Genotype III included only one strain with 339 aa and 4 PRRs. The meq aa profile of this genotype closely resembled genotype I, differing only in the presence of the consensus proline 194 instead of the deletion. Genotype IV included only one strain with 398 aa, 7 PRRs, and the 194delP deletion. This strain presented the same profile as the CVI988 vaccine and a Colombian strain, including polymorphisms A71S, T326I, and 194delP.

The results of the Fisher’s exact test between the two most abundant Brazilian genotypes (I and II) showed a significant difference (p < 0.00001) in the two genotypes when compared with respect to the MD form. It was observed that genotype I was associated with the visceral MD form, while genotype II was associated with the neural MD form.

3.3. Selective pressure analysis

The codon-based Z-test displayed positive selection with a dN/dS ratio of 2.07 indicating a higher rate of non-synonymous than synonymous substitutions (p-value < 0.02). Tajima’s test of neutrality resulted in D −2.45, suggesting that the changes in the meq gene were not random, but rather the result of positive selection.

According to codeml, the overall dN/dS for selected sites meq sites was 1.625, indicating evolution under positive selection (Table 2). The null models were rejected (p < 0.0001) and ten sites have evolved under positive selection with BEB > 0.95. Furthermore, BUSTED also found evidence of positive selection in meq (p-value = 0.0014). The pervasive diversifying selection was detected in 14 sites for FUBAR, 3 sites for FEL, and 2 sites for SLAC. Finally, episodic diversifying selection was detected in 9 sites by MEME (Table 3). The distribution of positive selective pressure sites along meq protein is shown in Figure 2 and were located mostly in ZIP and TAD domains.

Figure 2. Distribution of positive selective pressure sites in the meq oncoprotein. The proline-glutamine-rich (pro/gln), the basic region (BR), the leucine zipper (ZIP), and the transactivation domain (TAD) are indicated. Pins with yellow and green diamonds indicate the sites under pervasive or episodic positive selection, respectively. Pins with red circles indicate the sites under pervasive positive selection detected in Brazilian strains.

Table 2. Overall dN/dS rate and likelihood ratio test (LRT) for positive selection for the MDV meq gene estimated by codeml.

Group	Overall dN/dS^A	Model comparison^B^,C	2Δl^D	df^E	p-value
Pathotypes	1.625	M1a vs M2a	39.6174	2	2.00E-09
		M7 vs M8	39.6421	2	2.00E-09
		M8 vs M8a	39.6174	1	< 0.00001

^A dN/dS: Rate of nonsynonymous (dN) to synonymous (dS) substitutions.

^B Null models (Neutral): M1a. M7. M8a.

^C Alternative model (Selection): M2a. M8.

^D Δl: Difference between the likelihood scores.

^E df: Degrees of freedom.

Table 3. Sites under pervasive and episodic diversifying (positive) selection pressure in the meq gene estimated by datamonkey and PAML methods.

Codon position	Datamonkey				PAML
	FUBAR^A Probability α < β	SLAC^B P-[dN/dS > 1]	FEL^B p-value	MEME^B p-value	Codeml^C Pr(ω > 1)
					M2a	M8a
2	0.945	–	–	–	–	–
3	-^D	–	–	0	–	–
77	0.965	–	–	–	0.975	0.987
80	0.903	–	–	–	0.974	0.986
88	0.991	–	–	0.03	–	–
93	0.907	–	–	–	–	–
101	0.925	–	–	–	–	–
115	0.955	–	0.0006	–	0.998	0.999
119	–	–	–	–	0.974	0.986
139	0.975	–	–	–	–	–
153	–	–	0.0004	–	0.979	0.989
176	0.998	0.0392	–	0	0.971	0.985
180	0.973	–	–	0.08	–	–
194	0.944	–	–	0.1	–	–
217	0.998	–	–	0	0.997	0.999
277	0.903	0.0585	–	–	0.968	0.983
283	–	–	–	–	0.971	0.985
285	–	–	–	0.01
320	–	–	–	–	0.973	0.986
326	–	–	0.0002	–	–	–
329	–	–	–	0	–	–
339	0.943	–	–	0.09	–	–

^A FUBAR (p > 0.9).

^B MEME, FEL, SLAC (p < 0.1).

^C codeml p > 0.95.

^D -: no signal detected.

3.4. Phylodynamic analysis

This dataset comprised 96 meq sequences (Figure 3). The estimated mean tMRCA for the total tree was 1963-09-17 [95HPD: 1962-10-30, 1964-01-01]. The Brazilian strains from these epidemiological waves clustered together and their estimated tMRCA was 2016-05-20 [95HPD: 2015-02-13, 2017-08-09]. These strains would have originated from a common ancestor with the RB-1B strain approximately in 1978-07-18 [95HPD: 1975-01-18, 1980-12-31]. Since then, other strains from this cluster have been identified in Anseriformes in Japan (2005) and China (2013), and, progressively since 2010, in chickens in China, India, Iraq, and Italy.

Figure 3. Time-resolved maximum clade credibility tree (MCC) of non-redundant retrieved meq sequences plus Brazilian MDV sequenced strains. Most relevant dates and 95HPDS (tMRCA) are indicated with red lines. Blue bars indicate the 95HPD for each node time. Brazilian MDV strains are shaded in yellow.

The global MDV population dynamics exhibited two expansive waves with approximate duration of 10 years according to the HPD prior to the 1970s and 1980s, until the years 1980s and 1990s, respectively, and followed by a stable effective population size (Ne) until a decline before 2020 (Figure 4). On the other hand, the Brazilian MDV population belonging to the current epidemiological waves shows a notorious increase in this short period evaluated. We estimated the mean evolutionary rates by domains and compared them to explore their specific contributions to the total meq gene. The obtained rates were 0.8 × 10⁻², 1.9 × 10⁻², 4.9 × 10⁻², and 0.4 × 10⁻² for Pro/Gln, BR, ZIP, and TAD domains, respectively (Figure 5).

Figure 4. Bayesian skyline plot of the world vs Brazilian MDV strains. The y-axis represents the log number of the effective (Ne) population and the X-axis represents the period of years. The center black lines represent the mean of the Ne and the shaded area the 95% HPD of the Ne.

Figure 5. Mean evolutionary rates by meq domains. Rates represent the number of substitutions by site by year. The proline-glutamine-rich (pro/gln), the basic region (BR), the leucine zipper (ZIP), and the transactivation domain (TAD) are colored in salmon, green, light-blue and fuchsia, respectively.

3.5. Phylogeographic analysis

The geographical distribution of Brazilian genotypes is represented in Figure 6 according to the place of collection (Brazilian state) of the sequenced strains. Genotype I (21 strains) was broadly distributed across five different Brazilian States (Minas Gerais, Paraíba, Pernambuco, São Paulo, and Sergipe) all presenting exclusively the visceral MD form (Table S2). Conversely, genotype II (7 strains) was limited to São Paulo State, only showing the neural MD form. Genotypes III and IV were found only in Sergipe and Minas Gerais States, respectively.

Figure 6. Genotype distribution of marek’s disease virus strains throughout Brazil. The size of the circles is proportional to the number of available sequences.

The phylogeographic analysis encompassed 381 MDV meq sequences collected from 1962 to 2022, facilitating the visualization of global transmission routes for MDV strains (Figure 7, Figure S1). The outcomes of this analysis unveiled that Brazilian MDV strains potentially originated from three distinct routes. Specifically, the ancestral origin of Brazilian genotypes I and III can be traced back to the United States, which subsequently spread to other countries such as Hungary, Turkey, India, Japan, and eventually Brazil (Figure S2). Remarkably, within this cluster, the GA strain emerges as the oldest member, with its origin dating back to 1964. On the other hand, Brazilian genotype II would be originated from an ancestor in Japan which was dispersed to Italy, USA, and Brazil, and includes a Japanese relative strain (Tokachi-m2) isolated from Anas platyrhynchos (Figure S3).

Figure 7. Global phylogeographic distribution of marek’s disease virus strains. Colored circles are proportional to the number of available sequences. Connecting lines are colored in accordance with the country of origin for each route of transmission.

4. Discussion

MDV is an avian pathogen responsible for detrimental outbreaks in the poultry industry. Since the virus’s early pathological descriptions, its severity and acuteness have been constantly increasing (Osterrieder et al. 2006). The worldwide use of live commercial vaccines has exerted evolutive pressure on the field strains, resulting in insufficient protection from current vaccination programs (Nair et al. 2020). According to a retrospective analysis of cases in Brazil, it is estimated that MDV has circulated in a similar way to the rest of the world, being endemic during the 1960s until the introduction of the HVT vaccine at the end of the 1970s. The increasing cases of vaccination failure until the end of the 1980s gave rise to the introduction of the CVI988/Rispens vaccine in the 1990s (Katayama Ito NM 1999). Historically, Marek’s disease (MD) has been limited or effectively ­controlled, with minimal reports based on histopathology since 1999 (de Sousa 2010), and PCR between 2014 and 2016 (Blume et al. 2016; Torres et al. 2019). Torres et al. (2019) reported the widespread occurrence of pathogenic MDV with increased virulence in the field. Furthermore, since 2018, the increase in Marek’s disease (MD) cases referred to the Laboratory of Avian Diseases at the School of Veterinary Medicine—University of São Paulo in vaccinated poultry farms across Brazil has motivated this study. Recognizing this knowledge gap, we aimed to deepen the understanding of MDV dynamics in Brazil and elucidate potential factors contributing to the observed increase in outbreaks surge despite vaccination efforts.

The presence of field-type MDV was confirmed in all the farms analyzed (94.83% of the samples). Other oncogenic viruses were evaluated, resulting in the detection of REV in 8.6% of the samples, corresponding to two farms that presented the neural form of MD. Previous studies have shown the presence of REV in commercial birds in Brazil (Chacón et al. 2019, 2020, 2022). Furthermore, the co-infection of MDV and REV is a well-documented event and may cause a synergistic effect (Cui et al. 2010; Sun et al. 2017; Chacón et al. 2019). However, the co-infection status in the present study did not appear to worsen the clinical picture.

To gain deeper insight into the clinical presentation of MD, we categorized our studied strains into visceral or neural forms based on the main pathological findings, attempting to correlate genotypes with clinical information. This categorization can provide valuable insights into the relationship between MDV genotypes and clinical manifestations, contributing to a better understanding of disease dynamics in Brazil.

Numerous studies have employed phylogenetic inferences of meq sequences to classify MDV strains into clusters (López-Osorio et al. 2017; Mescolini et al. 2020a; Deng et al. 2021; Ozan et al. 2021; Kannaki et al. 2022). In this study, we incorporated a comprehensive set of complete meq sequences along with available epidemiological information to broaden the coverage. Our strains were included in three clusters (C2, C3, and C5). Cluster 2 is the most polymorphic and higher PRR-containing cluster (2 to 9), where several strains present 59 aa insertion in the C-terminal region with a suppressive effect on the meq expression and oncogenicity (Chang et al. 2002; Shamblin et al. 2004) including the less virulent pathotyped strains including vaccinal and other strains of reduced virulence (Renz et al. 2012; López-Osorio et al. 2017; Mescolini et al. 2019). Also, several Brazilian strains were included in cluster 3 together with highly virulent MDV strains from the USA (Shamblin et al. 2004; Witter et al. 2005). The clustering of Brazilian MDV strains within this cluster suggests a potential origin from virulent strains exhibiting the visceral MD form with multiple visceral lymphomatosis. Interestingly, all clusters with Brazilian strains included also Japanese strains isolated from migratory Anseriformes (Murata et al. 2012).

Extensive research has been published on MDV virulence classification based on the number of PPPs repeats of meq gene (Qian et al. 1995; Chang et al. 2002; Nair and Kung 2004; Shamblin et al. 2004; Renz et al. 2012; Mescolini et al. 2020b, 2022; Davidson et al. 2024). Apparently, a lower number of PPP repeats is associated with greater virulence (Chang et al. 2002; Shamblin et al. 2004). Furthermore, strains with reduced virulence present additional PPPs within a 59–60 amino acid insertion in the transactivation region, were associated with a repressive effect (Qian et al. 1995; Chang et al. 2002; Nair and Kung 2004; Shamblin et al. 2004; Renz et al. 2012; Yu et al. 2013). In the case of the Brazilian MDV genotypes studied (I, II, III, and IV), they presented 3, 4, 5, and 7 PPPs, respectively, suggesting a variety of virulence levels in circulating strains.

Moreover, studies have reported associations between meq amino acid substitutions and the virulence and oncogenicity of MDV (Shamblin et al. 2004; Conradie et al. 2020; Song et al. 2022). These studies suggest that novel polymorphisms P153L and 194delP, observed in genotypes I and IV of MDV Brazilian strains, may be involved in increasing virulence. However, the association molecular characteristics such as number of PPs and substitutions in the Meq protein needs further validation through in vivo experiments. Additionally, it is noteworthy that these newly identified polymorphisms are not present in MDV Brazilian strains isolated in 2014 (Blume et al. 2016; Torres et al. 2019).

Moreover, it has been observed that MDV strains in Brazil exhibited other amino acid changes as E77K and C119R, consistently reported as recurrent substitutions in highly virulent vv + MDV strains (Shamblin et al. 2004; Zhang et al. 2016; Conradie et al. 2020; Song et al. 2022). In contrast, the polymorphism A71S, which is absent in Brazilian genotypes I, II, and III, has been described in MDV strains with low virulence (Shamblin et al. 2004; Zhang et al. 2016; Conradie et al. 2020; Song et al. 2022). These aa changes serve as potential important classifier markers for MDV strains with different levels of virulence within the Brazilian context.

The molecular characterization of the meq gene is considered an accurate approximation and rapid method for determining MDV virulence compared to the state-of-the-art method of MDV pathotyping (Witter et al. 2005). This approach allows for the assessment of specific molecular characteristics associated with virulence, providing a valuable alternative for assessing the pathogenic potential of MDV strains more efficiently (Shamblin et al. 2004; Conradie et al. 2020; Davidson 2020). Although it is true that the present study did not include the evaluation of virulence through the gold-standard pathotyping, the molecular characteristics found in the different genotypes detected suggest the presence of various levels of virulence among the circulating strains. Nevertheless, we were able to show a clear association between the most frequent genotypes (I and II) and their visceral and neural MD forms, respectively.

Positive selection analysis of the meq gene revealed overall signals of positive selection evolution including episodic and pervasive selection. These results are consistent with the study of Padhi and Parcells (Padhi and Parcells 2016). Interestingly, detected sites under positive selection are involved in processes of cell oncogenesis (88 A) (Brown et al. 2009), the level of pathogenicity (77E, 80Y, 115 A, 153, 194, and 217 P) (Shamblin et al. 2004; Zhang et al. 2016; Conradie et al. 2020; Song et al. 2022), and even forming part of epitopes (77E and 80Y) (Lee et al. 2003; Sun et al. 2021). In particular, six of the 22 sites under positive selective pressure were detected in sequences of Brazilian strains, mainly in genotype I strains (which presents five of these sites). This suggests that the Brazilian strains of this genotype are evolving under selective pressure.

The substitution rate of meq was reported to be higher than that of typical dsDNA viruses and even RNA viruses (Padhi and Parcells 2016; Trimpert et al. 2017). Our results confirm this estimation. Moreover, we compared the substitution rate per domain and observed a closer similarity among all domains, with slightly faster rates for ZIP and BR domains. Selective pressure could potentially influence these rates (Padhi and Parcells 2016). In addition, the ZIP and BR domains play a crucial role as part of epitopes in the function and evolution of meq (Lee et al. 2003; Sun et al. 2021). Therefore, the combination of high substitution rates and processes of selective pressure could accelerate the variation of the meq gene enabling the emergence of new high-virulent strains (Trimpert et al. 2017).

The evolutionary analysis of the meq gene revealed that the Brazilian MDV strains investigated in this study originated from a common ancestor in 2016, with the emergence of genotype I becoming predominant since 2018. The end of the viral diversity expansive waves shown in the Bayesian skyline plot of meq synchronized with the commercial use of the HVT vaccine (until the end of the 1970s) and the CVI988 vaccine (at the beginning of the 1990s) (Reddy et al. 2017). These episodes could imply that vaccination temporarily reduces viral diversity but, instead, results in selective pressure and subsequent vaccine failure (Padhi and Parcells 2016).

Our phylogeographic analysis revealed that all examined genotypes were restricted to specific regions, except for genotype I, which exhibited a wider distribution across various Brazilian coastal states. This finding raises concerns since genotype I not only represents the most prevalent but also the most pathogenic genotype. The widespread presence of this highly virulent genotype throughout multiple regions underscores the significance of addressing and monitoring its impact on poultry health and the potential for further transmission.

Regarding the potential origin of MDV strains, our findings suggest potential hypothesis and scenarios. The most prevalent Brazilian strains are related with strains from Japan and the USA. A previous phylogeographic study reported MDV migration routes from the USA to China and then to Japan (Deng et al. 2021). Our study expands this view to a global level, including many other routes. Additionally, our proposed MDV migration is supported by the close phylogenetic relationship between Brazilian genotypes and other sequences in the analysis. Genotype II exhibited a phylogenetic association with a Japanese strain, while genotypes I and III were close to American and Japanese strains. The proximity of oncogenic viruses between the USA and Brazil was also previously documented, suggesting a possible etiological relationship (Chacón et al. 2022). The obtained potential routes associating the Brazilian MDV strains with those from Japan and the USA are probably related to commercial exchange, including the trade in poultry meat and derivatives, given Japan’s high interaction with Brazil (ABPA 2023), or the importation of chicken lineages or vaccines for avian use in Brazil from the USA (Ferreira 2011). However, it is important to highlight as a limitation that our analysis is based exclusively on the meq gene. The use of complete genomes could clarify and optimize these analyses.

In conclusion, we detected the presence of various circulating MDV strains and genotypes in Brazil. These strains caused an increase in MDV outbreaks in vaccinated poultry farms. The present study combined extensive molecular and evolutionary analysis with a simple phenotypic classification to demonstrate the relevance of a multifaceted approach. We highlighted the biological significance and impact on the poultry industry by focusing on the molecular and clinical findings related to the Brazilian genotype I. Notably, this genotype, which displayed the most severe clinical condition (visceral MD form), has emerged as the most prevalent and geographically widespread, exhibiting signs of evolution under positive selective pressure. Furthermore, we found phylogenetic and phylogeographic linkage between strains from the USA and Japan. Additional studies as well as the inclusion of more Brazilian sequences from different geographic and temporal spaces are necessary to confirm the potential associations between the molecular characteristics described in the Brazilian strains of the present study and their clinical implications in MD pathology.

Ethics statement

This study was approved by the Ethics Commission on Animal Use of the School of Veterinary Medicine, University of São Paulo (FMVZUSP), under CEUAVET protocol no. 1727010620.

Acknowledgments

The authors are grateful to the staff of the farms for supplying the biological material and epidemiological information.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The obtained sequences were submitted to GenBank under the accession numbers ON228225 to ON228254.

Additional information

Funding

This work was supported in part by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Grant [2017/05296-6]. A.J.P.F. was supported by the Conselho Nacional de Pesquisa e Desenvolvimento Tecnológico (CNPq) under Grant 301084/2019-0.