Phylogenetic and mutational analysis of bovine leukemia virus (BLV) tax gene in specialized dairy production systems in Antioquia, Colombia

Abstract

The bovine leukemia virus (BLV) is a pathogen of high importance for the dairy industry. Currently, twelve genotypes have been described worldwide with different pathogenicity and virulence, so it is critical to evaluate the circulating genotypes in each country/region to associate this information with risk situations. The aim of this work was to perform a phylogenetic and mutational analysis of the BLV tax gene in cows that belong to specialized dairies in the Department of Antioquia, Colombia. A conventional PCR for the tax gene was performed on 86 bovine samples. Sanger sequencing was carried out on 22 PCR products with a size of 959 bp. The sequences obtained were aligned and analyzed using the Maximum Likelihood and Bayesian phylogenetic approaches. A predictor was used to analyze the possible impact of amino acid substitution on the Tax structure and function. Although all sequences were found to belong to genotype 1, four of the 22 sequences were grouped into a different subclade G1A. Fifty percent of the samples showed punctual mutations in their amino acids. Mutation S104L was identified as "possibly harmful," while the V146A change found in all subclade G1A samples was identified as "possibly benign." Although further studies are necessary to determine whether there is an effect of these mutations on the development of the disease, this study presents part of the evolution of the virus and the changes at the amino acid level that are occurring in cattle from specialized dairy farms in Antioquia.

Introduction

Bovine leukemia virus (BLV) is the etiological agent of enzootic bovine leukosis (EBL), a disease of economic importance in bovine farming due to the negative impact it has on the productive and reproductive parameters of infected cattle (Bartlett et al. 2020). This virus is classified within the Retroviridae family, Deltaretrovirus genus, and is genetically related to the human T cell lymphotropic virus type 1, 2, and 3 (HTLV-1, HTLV-2, and HTLV-3 (Abdala et al. 2019). BLV has two copies of its genome (ssRNA) and is characterized because, in addition to the genes that code for structural proteins, it has a region that codes for accessory Tax and Rex proteins, which participate in the transcriptional and post-transcriptional regulation processes of the virus, and are essential in the infection and oncogenicity processes (Pluta et al. 2020). In South America, BLV is widely distributed and has been found in all countries where studies have been conducted (Selim et al. 2021). In Colombia, genotypes 1, 2, 3, and 6 have been identified as circulating in the national territory (Benavides et al. 2017, Corredor-Figueroa et al. 2020, Úsuga-Monroy et al. 2023). The availability of complete BLV sequences from different genotypes worldwide can define genotype-dependent pathogenesis and the association between genetic variability in each genotype and its infectivity (Polat et al. 2017). The tax gene, which codes for the Tax protein, has essential functions during viral replication and pathogenesis progression (Pluta et al. 2021). It is also one of the most polymorphic gene regions of the viral genome compared to the other genes (Polat et al. 2016). The Tax protein is a transcriptional activator of the viral genome and is essential in the infection process. Mutations in the sequence of this gene are associated with increased white blood cell count in infected cattle during the persistent lymphocytosis stage (Zyrianova et al. 2020), increased proviral load, and improved virus transmissibility (Pluta et al. 2020). Despite the importance of this protein for BLV, there are few studies on the genetic variations that this region has (Pluta et al. 2020, Buehring et al. 2019). The above added to recent studies that have established the zoonotic potential of BLV through the identification of different genes of the BLV provirus in samples of both benign and malignant human breast tissue and blood samples (Buehring et al. 2019, Le et al. 2020), demonstrates the need to delve deeper into the functions and mutations of the tax gene in BLV. Accordingly, the aim of this study was to perform a phylogenetic and mutational analysis of the BLV tax gene in naturally infected cows that belong to specialized dairies in the Antioquia Department in Colombia.

Material and methods

Ethical aspects

This research is part of the macro-project "Bovine leukosis in diaries of Antioquia: Evaluation of the zoonotic potential and the effect on reproductive performance," endorsed by the Institutional Committee for the Care and Use of Animals (CICUA) of Universidad Nacional de Colombia [020-2020].

Samples

The samples selected for this study were 86 bovines positive for BLV infection that were identified using the env molecular marker obtained from the work carried out by Castillo-Rey et al. (2023). In that study, the sample size was 575 milking cows from specialized mechanical and traditional dairies from three regions of Antioquia (North, East, and Valle de Aburrá) in Colombia. Figure 1 shows the geographical location of the samples analyzed in the current study.

Figure. 1. Georeferencing of the samples of the Department of Antioquia, Colombia. The map represents the geographical location of the analyzed samples from three regions of Antioquia. North: the lilac dot represents the samples from San Pedro (2 positive samples out of 12); the purple dot represents the samples from Entrerrios (1 positive sample out of 6); the pink dot represents the samples from Santa Rosa (4 positive samples out of 19), and the blue dot represents the samples from Belmira (2 positive samples out of 8). Valle de Aburrá: the orange dot represents the samples from Medellín (2 positive samples out of 16); the brown dot represents the samples from Bello (3 positive samples out of 11); and the green dot represents the samples from Girardota (6 positive samples out of 10). East: the purple dot represents the samples from the el Santuario area (2 positive samples out of 4).

PCR for the tax gene

The PCR for the tax gene was performed in a final volume of 25 μL with 150 ng of DNA, 0.2 mM of dNTPs, 1X buffer, and 1.5 μL of each oligonucleotide (10 mM): taxFW and taxRV (Úsuga-Monroy et al. 2023). The reaction conditions were: 5 min at 94°C, 40 cycles of 30 s at 94°C, 30 s at 56°C, 3 min at 72°C, and 5 min at 72°C. At the end of the reaction, a 959 bp fragment of the BLV tax gene was obtained and verified on a 2% agarose gel.

Sequence analysis

Sanger sequencing was performed for each of the samples in both directions of the genome at a commercial company (Macrogen Inc., Korea). According to the quality score of the sequences obtained, the contigs were generated, edited, and assembled in the SeqMan Pro (DNAStar Lasergene™) and MegAligment programs (Kumar et al. 2020). The BLAST tool was used to create the datasets for the phylogenetic analyses, seeking to find the homology between the obtained sequences and BLV reference sequences.

Phylogenetic analysis

For the phylogenetic analysis of the BLV tax region, the sequences from this study were aligned with 39 partial sequences of the tax gene and deposited in GenBank for the region between nucleotides 7299 and 8139. Sequences from BLV genotypes 5, 7, 8, and 11 were not included because, to date, no complete genome sequences/tax fragments are available for these genotypes. Sequences were aligned using the MUSCLE procedure in the MEGA X V 10.2.6 software (Kumar et al. 2020), and a maximum likelihood analysis was carried out using the best substitution model according to the Akaike Information Criterion (AIC). The bootstrap values were determined with 1000 repetitions, and only values higher than 70% were considered significant. A second phylogenetic tree was constructed through Bayesian methods using the MrBayes V 3.2.2 software (Ronquist et al. 2012). The evolution model utilized was selected based on the sampling parameters of the sample, and trees were sampled every 1000 generations. The consensus tree was obtained once the initial generations were discarded after 25%. The final trees obtained through both methods were edited with the FigTree V1.4.4 software (Rambaut, 2012).

Amino acid analysis of the Tax regulatory protein

The alignments of the partial amino acid sequences of the Tax protein from the sequences of the current study were predicted using the MEGA X V 10.2.6 software (Kumar et al. 2020). Genotype 1 sequences used in creating the phylogenetic trees were taken as reference sequences for the alignment. The amino acid substitutions from the viral sequences circulating in the Antioquia dairy were analyzed to determine if these amino acid changes were uniformly distributed in the three regions North, East, and Valley of Aburrá, or if, on the contrary, the substitutions were associated with a region, municipality, specific race. A maximum likelihood analysis was performed, and the best substitution model was chosen according to the Akaike information criterion (AIC). Bootstrap values were determined with 1000 repetitions, and only values higher than 70% were considered significant. An analysis was performed to predict the possible impact of an amino acid substitution on the structure and function of a protein using the PolyPhen-2 tool (http://genetics.bwh.harvard.edu/pph2/).

Results

This study only includes 22 sequences obtained from the tax gene due to the low quality of some of the PCR products obtained. As shown in Figure 2, the phylogenetic tree of the tax gene obtained by the maximum likelihood method presented a topology similar to the phylogenetic tree obtained by the Bayesian method in Figure 3. All the tax gene sequences analyzed in both trees were grouped in the BLV genotype 1 (G1) clade.

In both phylogenetic trees, two subclades of G1A and G1B are highlighted (Figures 2 and 3). The first highlighted in green has four sequences (OP414480, OP414490, OP414491, and OP414479) recorded in this study, along with three sequences obtained from the north and Valle de Aburrá areas of the Department of Antioquia (OP328415, OP328414, and OP328401) (Úsuga-Monroy et al. 2023). Therefore, it could be considered a new subtype of genotype 1 (G1A), considering branch supports of 87% obtained with the maximum likelihood and 100% with the Bayesian methods. The other subclade highlighted in pink locates three study samples (OP414439, OP414487, and OP414486) with branch supports of 85% obtained with the maximum likelihood and 86.92% with the Bayesian methods.

From the multiple alignment of the nucleotides of the 22 sequences obtained in the study, pairwise genetic distances (p) were determined, fluctuating between 0.000 and 0.007 with an average of 0.003. These distances are an estimator of the mutation rate of the population and a standard measure of variation. There are five sequences between which the highest difference was found in the nucleotide sequence, including OP414472, OP414479, OP414480, OP414484, and OP414486. However, these are grouped in common clades, and only OP414480 is within a subclade. Therefore, the type of mutations resulting from this genetic distance must be analyzed.

An alignment of the partial amino acid sequences of the 22 samples of this study was carried out together with 17 sequences from pure Holstein cattle from 2017 from Antioquia (Úsuga-Monroy et al. 2023), one from Uruguay and one from Paraguay belonging to G1, to identify possible significant variations in the amino acid profile of the samples analyzed in the study. The analysis included a total of 246 amino acids located in the region from amino acids 18 to 298 of the BLV Tax protein. Some samples obtained in this study were observed to have two mutations of interest in the Tax protein of the bovine leukemia virus. At position 104, a change from Serine to Leucine (S104L) was found, which was reported in three sequences, and at position 146, a change from Valine to Alanine (V146A) was registered in four of the study sequences and three of the sequences belonging to the Department of Antioquia (Úsuga-Monroy et al. 2023). The results of the amino acid profile were corroborated by creating a phylogenetic tree with the amino acid sequences obtained in the study and the reference sequences from GenBank for genotype 1 (Figure 4).

According to the analysis of the probability that a substitution is harmful (PolyPhen-2), the S104L change was found to be associated with a possibly harmful mutation with a score of 0.663 (sensitivity: 0.86; specificity: 0.91). In contrast, the V146A mutation is benign with a score of 0.001 (sensitivity: 0.99; specificity: 0.15) (Figure 5).

Figure. 2. Phylogenetic tree of the bovine leukemia virus (BLV) tax gene by the Maximum Likelihood method. The evolutionary history of the tax gene was reconstructed using the GTR+G+I model. Bootstrap values (1000 replicates) were plotted near the nodes. The tree is represented to scale where branch lengths are measured in the number of substitutions per site (Bar: 0.005 substitutions per site). The analysis was based on a combination of 39 partial sequences obtained from the GenBank database and 22 partial sequences isolated in the current study from three regions (north, east, and Valle de Aburrá) of the Department of Antioquia, Colombia, marked with (●). Genotypes (G) are indicated with numbers next to the sequence name, and clades are highlighted with different colors on the branches (blue G1, green G2, pink G3, sky blue G6, purple G9, orange G10, and black G4. Only one sequence is included). The data set used to build the tree consists of 843 positions. All evolutionary analyses were carried out using the MEGA X V.10.2.6 software.

Figure. 3. Phylogenetic tree of the bovine leukemia virus (BLV) tax gene obtained by the Bayesian method. The GTR+G substitution model was used. Bootstrap values (1000 replicates) are shown near the nodes. The tree is drawn to scale, with branch lengths measured as the number of substitutions per site (Bar: 0.002 substitutions per site). The analysis was based on a combination of 39 partial sequences obtained from the GenBank database and 22 partial sequences isolated in the current study from three regions (north, east, and Valle de Aburrá) of the Department of Antioquia, Colombia, marked with (●). Genotypes (G) are indicated with numbers next to the sequence name, and clades are highlighted with different colors on the branches (blue G1, green G2, pink G3, sky blue G6, purple G9, orange G10, and black G4). The fragments used to obtain the tree had a total of 843 positions in the final data set. Evolutionary analyses were performed in MrBayes V3.2.2., and the final tree was edited in FigTree V1.4.4.

Figure. 4. Phylogenetic tree of the evolution of the BLV Tax protein for genotype 1 using the Maximum Likelihood method, including the Jones-Taylor-Thornton (JTT) model with a discrete Gamma distribution of rate differences between sites. The evolutionary history of the Tax protein was reconstructed using the JTT model. Bootstrap values (1000 replicates) were plotted near the nodes. The tree is represented to scale, where branch lengths are measured in the number of substitutions per site (Bar: 4.0 substitutions per site). The analysis was based on a combination of 24 partial sequences obtained from the GenBank database and 22 partial sequences isolated in the current study from three regions (North, East, and Valle de Aburrá) of the Department of Antioquia, Colombia, marked with (●). A total of 269 amino acids were involved in the final data set. All evolutionary analyses were carried out using the MEGA X V.10.2.6 software. The clade with the sequences with the S104L mutation is highlighted in green, and the sequences with the V146A mutation are in pink. All evolutionary analyses were carried out using the MEGA X V.10.2.6 software, and the final tree was edited in FigTree V1.4.4.

Figure. 5. Mutant protein sequences of the BLV tax gene and the Tax protein structure diagram using the PolyPhen-2 mutation analysis. A. Mutant S104L protein sequence, B. Mutant V146A protein sequence, and C. Tax protein structure diagram. a: amino-terminal zinc finger, b: leucine-rich activation domain, c: multifunctional domain, P: phosphorylation site, S1: S104L amino acid substitution, S2: V146A amino acid substitution.

Discussion

The BLV tax oncogene encodes the Tax protein, which has a structure of 309 amino acids, and its functions include interfering with DNA repair mechanisms, preventing apoptosis, and inhibiting tumor suppressor genes (Gao et al. 2020). This genome region shows high genetic variability, so it is not commonly used to diagnose infection (Dao et al. 2018). Previous studies of the molecular prevalence of BLV using different genes have shown that the diagnosis by PCR of this region of the genome is less precise than that obtained using the env gene, finding a smaller number of positive animals with tax (Corredor-Figueroa et al. 2020, Canova et al. 2021). The above is consistent with the results obtained in this study, where only 45% of the animals positive for BLV due to the env gene were also positive due to the tax gene. Therefore, it is not widely used in diagnosis. However, some problems have been reported in the characterization of some genotypes with the amplification of the 444 bp segment of the env gene (Buehring et al. 2014), which is the basis for the classification of the 11 currently known genotypes of BLV and is one of the reasons why discordances can be found when performing phylogenetic studies with other regions of the genome such as the tax gene (Zyrianova et al. 2020, Dao et al. 2018).

In addition to being an important molecular marker for phylogenetic studies of the virus, mutations in the tax gene are relevant to understanding the development of the disease since it has been found to have relevance in it (Pluta et al. 2020, Zyrianova et al. 2020). The 22 tax gene sequences analyzed in the current study belong to genotype 1 (Figures 2 and 3). The phylogenetic tree obtained by the Maximum Likelihood method (Figure 2) coincides in the groupings obtained by the Bayesian method (Figure 3), registering in both branches supports values mostly higher than 70%. Furthermore, the 22 sequences obtained in the study from cattle of different breeds and racial crosses (Holstein, Jersey, Viking Red, Ayrshire, and Jerhol) highly correlate with each other, grouping and forming subclades within genotype 1 sequences. Further, there is also a closeness to the sequences obtained in the Department of Antioquia in a study carried out only with cattle of the Holstein breed (Úsuga-Monroy et al. 2023). The samples are from the three study regions (North, East, and Valle de Aburrá), so the genotype G1 in the department is common. However, a more extensive analysis should be performed, including sequences from different geographic regions and different types of bovine production systems, to draw conclusions about the origin of this grouping and whether it could indeed be classified as a subgroup of the genotype.

An essential part of a phylogenetic study is the analysis of variations in the amino acid profile of the protein because it has been found that mutations in its amino acid sequence can translate into changes in the secondary and tertiary structure of the protein that can further lead to changes in virus pathogenicity and disease progression (Tajima and Aida 2000, Tajima and Aida 2002, Twizere et al. 2000, Inoue et al. 2013, Mori et al. 2019). Despite the importance of this protein for the pathogenesis of BLV, studies on the known genetic variations are very few (Pluta et al. 2020, Zyrianova et al. 2020, Tomiyasu et al. 2021), and this is the first report of amino acid variation of the Tax protein of samples from the specialized dairy of Colombia.

The first S104L mutation is located close to one of the two phosphorylation sites located at serine position 106 (Figure 5C). Although these phosphorylation sites are not crucial for the function of the protein in transactivation processes in vivo (Twizere et al. 2000), an amino acid alteration near this position can affect the structure of the protein, favoring or not the phosphorylation process, which may have effects on infection development. The PolyPhen-2 mutation analysis shows that the S104L change is associated with a possibly harmful mutation (score of 0.663) (Figure 5A). In this case, the change occurs due to a larger amino acid of different polarity, a possible interference point in the protein phosphorylation process. When analyzing the origin of the sequences with this mutation, they all came from the same region (Valle de Aburrá). Likewise, they all correspond to the F1 Jersey x Holstein racial component (Jerhol), so there may be a link between the sequences obtained, the breed, and the origin. Furthermore, these sequences are grouped in the subclade highlighted in green in Figures 2 and 3, being genetically close and separating themselves from the rest of the sequences in this same subgroup (G1A).

The V146A mutation of the Tax protein in the current study was identified in three sequences (OP414439, OP414487, and OP414486). These sequences are grouped in subgroup G1B by both phylogenetic reconstruction methods. Furthermore, this same amino acid substitution was found in three samples from the previous study carried out on pure Holstein cows from the Department of Antioquia (Úsuga-Monroy et al. 2023). In this work, there is a substitution of valine, an essential amino acid, with alanine, one that is smaller in size, both hydrophobic. This substitution of an essential amino acid for a non-essential one can drastically alter the secondary and tertiary structure of the Tax protein. However, the PolyPhen-2 mutation analysis shows that the V146L change is associated with a possibly benign mutation (score of 0.001) (Figure 5B). The regions of origin of the samples that presented this mutation are the north and Valle de Aburrá areas. Nonetheless, it cannot be ruled out that it could occur throughout the territory because the sample size for this analysis does not allow us to conclude otherwise. Likewise, the data do not allow us to establish whether there is any association between the racial component and the results because the samples studied were not classified.

A study carried out by Tajima and Aida (2020) revealed a correlation between the substitution of valine at position 146 with isoleucine in the BLV Tax protein and an increase in the transactivation activity of the protein. This increase in transactivation activity leads to greater transcription and production of viral particles. Because the substitutions found in this study were located in the same position (V146A), it is relevant to evaluate the possible effects of changing the amino acid from valine to alanine in this position. However, although both modifications imply a change in the amino acid in position 146, the introduction of alanine instead of isoleucine may have different consequences on the function and activity of the Tax protein.

In the phylogenetic tree constructed with the amino acid sequences (Figure 4), two clades in the group are formed by the previously mentioned mutations, S104L (green subclade) and V146A (pink subclade). Apart from this, the sequence from the United States (EF428061) is aligned in a clade apart from the others indicated in the figure in blue. This is due to having two unique mutations to this genotype in its amino acid sequence. This phylogenetic tree allows the identification of sequences related to having similar amino acid substitutions in greater detail. More studies are necessary to follow up on the mutations found (S104L and V146A) to validate if they can be correlated with a change in the tumor capacity of the virus or the development of the disease.

The molecular diagnosis of bovines infected with BLV using the tax marker presented fewer positive animals than that carried out using the env gene. However, this gene exhibits mutations that make the study of its genetic evolution desirable. The tax gene is useful in identifying the different genotypes of the bovine leukemia virus, and it is also important to monitor mutations in this region due to their influence on the development of the disease and virulence. However, it is recommended to perform phylogenetic analysis using the env gene sequences to corroborate that the samples undoubtedly correspond to the genotype and that the region used in this study serves as a molecular marker for phylogenetic analyses of BLV. According to the results of this study, in the specialized dairies of the Antioquia Department, there are currently two possible circulating subvariants of BLV genotype 1 (G1A and G1B), which involve changes in the amino acid profile of the BLV tax protein, finding the amino acid change S104L related to the Valle de Aburrá area of the Department of Antioquia and the Jerhol racial component, in addition to V146A that circulates in the north and Valle de Aburrá areas, without association with any racial component. Therefore, carrying out a more in-depth study about the possible effects of these amino acid changes on viral pathogenicity and disease development in cattle is essential.

Acknowledgments

The authors thank the specialized dairy producers, who allowed the collection of samples and the use of data from their herds to carry out this research. Finally, the authors thank the Technical Assistance Department of Colanta for its essential support in executing this project.

Funding

This research was funded by "Convocatoria conjunta de proyectos de I+D+i en el marco de la agenda regional de I + D -> i del grupo G8+1″ [Joint call for R&D+I projects within the framework of the regional R&D agenda -> i of the G8+1 group] and the Government of Antioquia. The financers had no role in the study design, data collection and analysis, publication decisions, or manuscript preparation.

Ethical approval

This research is part of the macro project "Bovine leukemia virus in dairies of Antioquia: Evaluation of the zoonotic potential and the effect on reproductive performance," endorsed by the Institutional Committee for the Care and Use of Animals (CICUA) of Universidad Nacional de Colombia [020-2020]. No permits or licenses were required for the development of the project.

Conflict of interest

The authors declare that they have no conflicts of interest.

Author Contributions

Conceptualization, C.U.M.; methodology, D.C.R., and C.U.M.; formal analysis, D.C.R., and C.U.M.; investigation, D.C.R.; writing original draft preparation, D.C.R.; writing, review and editing, D.C.R., C.U.M., and A.L.H.; visualization, C.U.M.; supervision, A.L.H.; project administration, A.L.H.; funding acquisition, A.L.H., and C.U.M. All authors have read and agreed to the published version of the manuscript.

Data availability

All data analyzed during the current study are available from the corresponding author upon request.