Phylogenetic analysis of small ruminant lentiviruses originating from naturally infected sheep and goats from Poland based on the long terminal repeat sequences

Maedi-visna virus (MVV) and caprine arthritis encephalitis virus (CAEV) belong to the Lentivirus genus of the Retroviridae family. First isolated respectively from sheep and goats, MVV and CAEV were considered to be strictly host specific. Today, these viruses are no longer considered species-specific pathogens, since they cross the species barrier between goats and sheep very efficiently and infect both. Therefore, MVV and CAEV are considered as a single group referred to as small ruminant lentiviruses (SRLVs).

Small ruminant lentiviruses are transmitted vertically via infected milk and colostrum uptake and horizontally through respiratory secretion when infected animals are in close contact with uninfected ones (10). These lentiviruses induce a chronic multisystemic disease with inflammatory lesions in the mammary gland, lungs, joints and brain. Symptoms are observed only in one third of the infected animals and are typically pneumonia, arthritis (enlargement of the carpal joints), mastitis or encephalitis (15). The retroviral SRLVs are characterised by a single-stranded RNA genome, which, once retrotranscribed by the viral reverse transcriptase (RT) enzyme, integrates into the host genome in the provirus form. The provirus can then undergo a productive replicative cycle or remain hidden or dormant inside the cell in a state called latency. The genome of SRLVs is comprised of three structural genes encoding the group-specific antigens (gag), the polymerase (pol) and envelope (env) genes, and accessory genes including vpr-like (formerly tat), rev and vif, which have regulatory functions. The pol gene encodes the reverse transcriptase (RT), protease (PR) and integrase (IN) enzymes involved in replication and DNA integration. The gag gene encodes internal structural proteins, which are the matrix (MA), capsid (CA) and nucleocapsid (NC) proteins. The env gene encodes surface and transmembrane glycoproteins. The proviral DNA of SRLVs is flanked by non-coding sequences known as long terminal repeats (LTRs). The LTR is subdivided into three regions, U3, R and U5, and contains principal regulatory sequences such as promoter and enhancer elements, including host transcription factor binding sites required for the replication of the virus and viral gene expression (5).

Based on the gag and pol nucleotide sequences, the current SRLV phylogeny consists of five groups (A–E) which are further divided into multiple subtypes. Group A SRLVs comprise the prototypic strain K1514 (A1 subtype), which was the first lentivirus strain characterised. This group is the most heterogeneous and holds at least twenty-seven subtypes (A1–A27), while group B is subdivided into five subtypes (B1–B5) containing the prototypic Cork CAEV strain belonging to subtype B1. Groups A and B are widespread, while groups C, D and E are restricted to certain geographical areas. Group C was detected only in Norway, group D was detected in Switzerland and Spain, and group E, which comprises two subtypes (E1 and E2), was isolated only in Italy (5, 10). The env gene and LTRs show a high level of variability, while the pol and gag genes are relatively conserved, and it is these which are mainly used for phylogenetic analyses of SRLVs. However, the phylogenetic analysis of variable regions is more informative than that of conservative regions (24). Moreover, the highest level of phylogenetic information for SRLVs seems to be accumulated in LTR sequences rather than in the other genomic regions of SRLVs (29). In Poland, SRLV infections are widespread, with herd-level seroprevalence of 33.3% and 71.9% in sheep and goats, respectively (8, 9). Previous studies on the genetic diversity of SRLVs from Poland provided evidence for the presence of two groups (A and B) and 12 subtypes within SRLVs. Polish SRLVs isolated from sheep and goats belonged to subtypes B1, B2, A1, A5, A12, A13, A16–A18, A23, A24 and A27 (17). The phylogenetic analysis of Polish strains has so far been performed based on gag gene or gag and env gene sequences; however, phylogenetic analysis of the LTR has never been performed. Thus, the aim of this study was to extend genetic/phylogenetic analysis of the previously identified SRLV Polish strains by adding LTR sequences. This study provides additional insight into the genetic diversity of SRLV field strains in Poland, their phylogenetic relationships, and their position in the recently established SRLV classification.

A total of 112 samples analysed in this study were collected from six different voivodeships of Poland over the 13-year period of 2008–2021. Previous studies made 61 LTR sequences available (16, 17, 18) while this study obtained 51 new LTR sequences. The samples originated from 51 sheep and 61 goats from 19 flocks of a single animal species (sheep or goats) or mixed species (goats and sheep). Anticoagulated blood was used as a source of peripheral blood leukocytes (PBL) which were isolated by the standard protocol (16). The phylogenetic affiliation of all samples was previously determined on the basis of the gag and/or env sequences. The sequences were assigned to SRLV subtypes A5, A12, A13, A16, A17, A18, A23, A24, A27, B1 and B2 (Table 1).

The genomic DNA was extracted from the original biological samples (leukocyte pellets) using a NucleoSpin Blood Quick Pure Kit (Macherey-Nagel, Düren, Germany), according to the manufacturer’s recommendations. The quality and quantity of DNA were evaluated in a nanophotometer (Implen, Munich, Germany). All methods were performed in accordance with the relevant guidelines and regulations. Specifically, blood collection was approved (no. 37/2016) by the Local Ethics Committee on Animal Testing at the University of Life Sciences in Lublin, Poland. The LTR U3-R region was amplified using a nested PCR protocol already described (25). The primer pair of LTREFW (5ʹ-ACTGTCAGGRCAGAGAACARATGCC-3ʹ) and LTRERV (5ʹ-CTCTCTTACCTTACTTCAGG-3ʹ) was used in the first round of PCR, and the LTRIFW (5ʹ-AAGTCATGTAKCAGCTGATGCTT-3ʹ) and LTRIRV (5ʹ-TTGCACGGAATTAGTAACG-3ʹ) pair was used in the second round. The PCR products were analysed by electrophoresis on a 1.5% agarose gel. Then, the PCR products were purified using NucleoSpin Gel and PCR Clean-up (Marcherey-Nagel) and directly sequenced on a 3730xl DNA Analyzer using the BigDye Terminator v3.1 Cycle Sequencing Kit (both products of Applied Biosystems, Foster City, CA, USA). The obtained SRLV sequences were edited and analysed using Geneious Pro 5.3 software (Biomatters, Auckland, New Zealand). All new sequences reported in this study were submitted to GenBank under accession numbers ON637539– ON637623. The consensus sequences were aligned to each other and with previously published SRLV sequences of genotypes A–E. Multiple sequence alignment was performed using multiple sequence comparison by log-expectation (MUSCLE).

Model testing was performed to select the best-fit evolutionary model based on the Bayesian information criterion (BIC) and Akaike information criterion (AIC). As best accorded with the results, the Kimura 2-parameter model with gamma distribution (+G) and five rate categories was used to construct phylogenetic trees using the neighbour-joining method (NJ), the maximum-likelihood method (ML), and the unweighted pair group method with arithmetic mean (UPGMA). These methods contrast a statistical approach (in ML) with distance-based methods (in UPGMA and NJ). The statistical confidence of the topologies was evaluated by nonparametric bootstrap analysis with 1,000 iterations. Multiple alignment, model testing, tree building and pairwise genetic distances were calculated with MEGA 7 software (11). The RDP4 application was used to perform a recombination analysis (12). The software used seven primary methods: RDP, GENECONV, BootScan, MaxChi, Chimaera, SiScan and 3Seq. Putative recombinant events were considered significant when P≤0.01 was observed for the same event using five or more methods.

The LTR fragment from 15 sheep and 7 goat samples could not be amplified. The phylogenetic analyses of the LTR region of 54 goat and 36 sheep sequences from Polish SRLV strains and 39 reference strains produced three phylogenetic trees with topologies that were almost identical when using the three different methods ML, NJ and UPGMA. The unrooted phylogenetic trees of LTR sequences are shown in Figs 1–3.

Fig. 1

Fig. 2

Fig. 3

Three main clusters were observed in the phylogenetic trees with clear separation between group A, group B and group E according to the reference sequences. Group A represented MVV-like isolates, while group B represented CAEV-like isolates. Group E comprised two subtypes (E1 and E2) isolated in Italy. Group B consisted of subtype B1, B2 and B3 strains. The LTR sequence of the Norwegian strain 1GA of group C clustered with strains belonging to group B. This strain was the most closely related to the Volterra and Fonni strains representing subtype B3; however, the mean genetic distance between these strains was 30%. Furthermore, the It009 strain, a member of subtype A20, was most closely related to subtype B1 sequences (mean genetic distance 18.4%). The highly divergent group A consisted of many branches. Reference sequences representing subtypes A1, A2/A3, A4, A8 and A19 formed separate clusters. The Polish strains were clearly separated from the strains belonging to groups B and E. Our analysis showed that all Polish sequences were affiliated to group A and grouped in at least 10 clusters. The sequences belonged to the known subtypes A1, A5, A12, A13, A16, A17, A18, A23, A24 and A27. In particular, the LTR sequence of the #3g90472 strain was closely related to the K1514 and LV-1 strains (mean genetic distance 1.4%) and clustered in subtype A1. The LTR sequences of the Polish #13g4742, #13g5994, #13g5826, #13g7592, #13g5819, #13g6038, #13g5962, #13g3038, #17g8008, #17g9692, #15g3540, #15g0788, #14s0334, #13g5870, #17g1318, #17g8039, #15g3535, #15g0599, #17g8046, #15s20 and #18s5023 strains clustered within subtype A5 showing a mean genetic distance of 3.3%. Affiliation of this cluster was supported by high bootstrap values of 89% by ML, 100% by NJ and 97% by UPMGA. The sequences of the Polish #15g8699, #15g9509, #15g0580, #5sTryk6, #5sTryk2, #15g0599(2), #15g9510, #15g3533, #14g7096, #14g8891, #14g7219, #14g6808, #14g7102 and #14g7134 strains belonged to subtype A12; however, the existence of this subtype was not confirmed by high bootstrap values. The intra-subtype similarity of sequences belonging to this subtype was 5.0%. The LTR sequences of the #6922, #9179, #8063, #6981, #9155, #1406, #1304 and #1911 strains originating from sheep from flocks no. 10 and 12 clustered together (mean genetic distance 1.6%) and formed subtype A13, supported by high bootstrap values of 99% by NJ and 99% by UPMGA. Strain #7g2993 was closely related to strain #8g1202 (nucleotide distance 9.6%) and together with this strain formed a separate cluster named A16, supported by high bootstrap values of 62% by ML, 79% by NJ and 89% by UPMGA. The sequences originating from goats from flock 9 (#9431, #8344, #0042, #5675, #6909, #5621, #5616, #1561, #3085, #5654, #5686, #1580, #8172, #1485 and #8370) clustered in subtype A17 (mean nucleotide distance 3.3%). Affiliation of this subtype was supported by high bootstrap values of 89% by ML, 99% by NJ and 99% by UPMGA. Sequences of this region of strains #1, #4, #16, #14, #3225, #3249 and #3188 originating from sheep from flocks 15 and 16 formed a separate cluster, subtype A24, supported by high bootstrap values of 69% by NJ and 80% by UPMGA. The sequences of the #goat2, #goat3, #goat4, #goat5, #goat6, and #goat7 strains clustered in subtype A27, supported by high bootstrap values of 91% by ML, 99% by NJ and 100% by UPMGA. The LTR sequences of strains #17s2590, #17s4315, #17s1622 and #17s3275 clustered in subtype A23, while those of strains #18s9855, #18g4464, #10s7041, #11s7020 and #11s0090 clustered in subtype A18.

The affiliation of the A23 subtype was supported by high bootstrap values of 98% by NJ and 83% by UPGMA while affiliation of the A18 subtype was supported by high bootstrap values of 71% by NJ and ML. The LTR sequences of strains #3, #16, #40, #13 and #33 originating from a sheep from flock 15 formed a unique cluster within group A, which was supported by high bootstrap values of 59% by ML, 77% by NJ and 95% by UPMGA.

Most of the analysed strains belonged to the same subtype by the indication of the gag, env and LTR genomic regions (Table 1). However, the affiliations of some strains were different. In particular, strain #0016 belonged to subtype A13 on the basis of its gag and env sequences, while on the basis its LTR this strain formed a separate cluster within group A, different from other subtypes. The LTR sequence of this strain was closely related to that of a strain belonging to subtype A12 (mean genetic distance of 10.5%). Strains #14s0334, #15g3540, #15g0788 and #15s20 had gag/env sequences which classified the strains as B2 but had LTR sequences which assigned them to subtype A5. Strain #15g0580, which had previously been classified as B2, affiliated to subtype A12 when heeding its LTR sequence. Furthermore, the gag/env sequences of strains #15s1, #15s4 and #15s6 were classified as A12 while the LTR sequences of these strains clustered in subtype A24. Strains #15s3, #15s13, #15s16, #15s33 and #15s40, of which the gag/env sequences indicated them to be in A12, revealed LTR sequences placing them in a new cluster within group A. The LTR sequence of strain #17s3275 confirmed its affiliation to subtype A23 predicated on its env sequence, although the gag sequence of this strain belonged to subtype B2. The gag sequence of strain #18s5023 matched subtype A24, the env sequence was not assigned to any known subtype, while the LTR sequence of this strain corresponded to subtype A5. Discrepancies in affiliation also occurred in strain #17s3691 (B2/B/A), #17s4018 (A23/B/A), #15g0599 (B2/A/A5/A12), #18s9855(B2/B/A18) and #18g4464 (A/A18). No clear affiliation was also observed in the case of strains isolated from animals co-infected with different subtypes. When the gag sequences were analysed, subtypes A12 and B2 were found in goat 3535 and sheep 14 from flock 15, while the LTR sequences isolated from these animals clustered in subtypes A5 and A24. Furthermore, goat 1202 from flock 8 co-infected with strains A12 and B2 (env/gag sequences) was apparently infected with a subtype A16 strain when the LTR sequence was consulted.

The WebLogo corresponding to generated alignments of LTR sequences representing SRLV subtypes A1, A2/A3, A4, A5, A8, A12, A13, A16, A17, A18, A19, A23, A24, A27, B1, B2, B3, C, E1 and E2 and the inferred U3/R and R/U5 boundaries are shown in Fig. 4. The LTR region analysed in this study contained two AP-1 sites, one AML(vis), one AP-4 motif and one TATA box in the U3 region and the AATAAA motif in the R region (Figs 4 and 5). This representation of the alignment revealed that the U5 region is generally more conserved than the R and U3 regions. The AP-4 site, the TATA box and the polyadenylation site (AATAAA) were the most conserved elements between all strains belonging to group A, B, C and E SRLVs. Our results revealed that strains of subtypes A17, A20 A27 and B3 had a unique T to A substitution in the fifth position of the TATA box. The CAEV-like sequences did not have an AML sequence close to the TATA box present in any MVV-like sequences. The AP-1 sites revealed rather group-specific conservation. Our results also showed that the LTR sequence reflected subtype-specific patterns of sequence diversity (Fig. 5). The GAGAAGCTTTG and TAAGAGCTTTG insertions next to the TATA box were only found in sequences belonging to subtypes B1 and B2, respectively, and CTTGCTACT and TCAGACGCT insertions were found exclusively in sequences belonging to groups C and E, respectively. Sequences from subtype A1 had a unique TCGAAG GAAAGA insertion in the R region, while the sequences from subtypes A17, B3 and C had respective unique CCGAAGGAAAG, TCGAAGGAAAGA and TCGAAGGAAAGAG insertions. The sequence from subtype A4 had a unique GCTTTGCC insertion, while the sequences from subtype A8 had a unique CTGGTCGC insertion close to the polyA site. All subtype A20, B1, and B2 sequences had a unique GTACCGAGACCT insertion located downstream of the polyA site. Two insertions located upstream of polyA, GATTGCC and TGCCGAGTG, were identified only in sequences belonging to subtypes B1 and B2.

Fig. 4

Fig. 5

Fig. 6

Some indications of recombination events were observed; however, these possible recombination events were confirmed only by two out of seven methods used in this study. On the basis of LTR alignment, two putative recombination events were detected. The MaxChi and SiScan methods detected a recombination event in almost all subtype A12 strains. In this recombination event, the beginning and end breakpoints were located at 88 and 264 nucleotides in alignments and the major and minor parents were the Turkish #Mk3 strain and an unknown strain, respectively (Fig. 6a). The MaxChi and 3Seq methods detected a recombination event in all subtype A5 strains between positions 258 and 114 in alignment with #It007 (subtype A8) as the major parent and the unknown strain as the minor parent (Fig. 6b).

The first attempt to establish the phylogeny of SRLVs was published in 1998 by Zanoni (29). At that time, SRLVs were phylogenetically classified into at least six clades, with no clear separation by host species or geographical origin (24, 29). The phylogenetic organisation of SRLVs proposed by Shah et al. (27) based on sequences of 1.2 kilobase (kb) pol and 1.8 kb gag-pol fragments sorted strains into four main groups, termed A–D. However, detection of these fragments with classifiable sequences often could not be achieved, most probably because of the high genetic heterogeneity of SRLV strains. Therefore, many other strains have been added to the phylogenetic tree using different fragments and regions. This led to the new classification that currently affiliates SRLVs into five groups (A–E) and at least 34 subtypes (A1–A27, B1– B5 and E1–E2). Previous studies using gag and env gene sequences revealed that Polish SRLVs isolated from sheep and goats belonged to subtypes B1, B2, A1, A5, A12, A13, A16–A18, A23, A24 and A27 (17). The LTR sequence analysis in this research elaborates the phylogeny of known SRLVs from this country and augments the genetic information for subtyping.

The phylogenetic trees constructed in this study showed the distribution of sheep- and goat-derived sequences in three out of five clades in the SRLV phylogeny (groups A, B and E). Group A was the most heterogeneous and consisted of many subtypes. Group B consisted of subtype B1, B2 and B3 strains, while group E comprised two subtypes isolated in Italy (E1 and E2). The phylogram of the LTR did not include group D because LTR sequences of strains belonging to this group were not available. Our results revealed that the LTR sequence of the Norwegian 1GA strain classified into group C by Shah et al. (27) clustered with strains in group B. This strain was the most closely related to the subtype B3 Volterra and Fonni strains. In this aspect, our results are in line with those of Mendiola et al. (13), who showed that the LTR sequences of the strains from genotypes C and B3 had a common root. Furthermore, in the phylogeny based on the gag or gag-pol sequence, the Norwegian 1GA strain also belonged to group B (14). Analysing the phylogenetic trees of SRLVs, it can be observed that the affiliation of some strains is inconstant and depends on the analysed fragment. For example, on the basis of the pol fragment, the Belgian VB.5.6 and LB.2.3 strains belong to the B5 subtype, while on the basis of the gag-pol region these strains belong to subtype B1 (14). In the gag/pol phylogenetic tree, the It009.2017 strain affililates to subtype A20, while in the pol tree it finds a place in subtype A1. Moving to the LTR sequence of this strain, it was most closely related to subtype B1 sequences, which may suggest that it is in fact a recombinant strain. A similar situation was observed in the case of B4 strains, which transpired to be A/B recombinant strains (26). In addition, sequences purported to belong to genotype D turned out to be genotype A strains, exhibiting divergence in the pol gene (21). Our results revealed that the Polish caprine and ovine LTR sequences clustered within group A and grouped into at least 10 clusters. Most of the Polish strains (78%) analysed in this study belonged to the same subtype when the gag, env and LTR genomic regions were scrutinised. Reina et al. (23) and Zanoni

(29) also plotted phylogenetic trees on the basis of different regions (gag, pol, env and LTR) which had topologies compatible with one another. Phylogenetic analysis of the LTR fragment confirmed the high heterogeneity of Polish strains which clustered in subtypes A1, A5, A12, A13, A16–A18, A23, A24 and A27. The identical affiliations of most Polish strains emerging from analysis of three different regions confirmed the existence of these subtypes. The topologies of LTR trees were almost identical even when derived from three different methods (NJ, ML and UPGMA), which also increases the reliability of the obtained results. The agreement of the gag, env, and LTR phylogenies suggests that these genomic fragments have co-evolved. The high SRLV heterogeneity in Poland is probably the result of the absence of any SRLV control programmes, and the lack of the data such a programme might have provided may cause molecular and serological misdiagnosis of infections with some divergent strains possessing a high number of mutations.

Discrepancies in affiliation were nevertheless observed in 24 (21%) Polish strains, most of which (22 strains) came from mixed flocks where more than one SRLV genotype circulated. This confirms that mixed flocks are an excellent environment for the emergence of new SRLV variants. The different affiliation of LTR versus gag/env sequences of some Polish strains may suggest that these strains resulted from cross-species transmission and later evolved to adapt to a new host. Erhouma et al. (2) showed that after transmission from goats to wild ibexes via natural contact, important genetic differences were present mainly in the LTR region of proviral SRLV sequences. Circulation of different subtypes in one flock provides an opportunity for co-infection and recombination. On the basis of the LTR sequence, co-infection with strains belonging to subtypes A5 and A12 was evidenced in goat 0599 from flock 15. Double infections are a prerequisite for recombination between heterologous or homologous SRLV strains and an important source of genetic variability. In this article, some evidence of recombination was seen in strains from subtypes A5 and A12; however, possible recombination events were confirmed only by two out of seven methods.

Our results revealed that the LTR sequence of SRLV genotypes A, B, C and E have several homologous regions, such as those encompassing the TATA box, the polyA-signal, the AP-4 site, and some stretches of the U5 region. The observed conservation of these sites confirmed previous findings suggesting that maintenance of these sites is essential to the transcriptional activity of the virus (3). Our results also showed that variations in the LTR sequence were subtype specific. Considerable sequence variation was also observed between different HIV-1 subtypes, and several subtype-specific sequence motifs were identified in the HIV-1 LTR region (22). Furthermore, it was revealed that the subtype-specific LTRs dictated a different replication rate and expression of HIV-1, which resulted in biological differences between subtypes (7). It is known that minor changes or rearrangements within the transcription factor binding sites in the LTR of retroviruses can have a significant impact on cell tropism and pathogenicity and affect viral fitness (7, 28). We suppose that the sequence variation between SRLV subtypes may also correspond to the different transcriptional activity of these strains, which may affect their phenotypes. The relationship between SRLV subtype, biological properties and pathogenicity is still unclear; however, in some cases a certain relationship can be observed. The presence of non-pathogenic SRLVs was reported in goats infected with viruses belonging to the A4 subtype, and it was evidenced that particular LTR mutations may explain their attenuated phenotype (1). The highly divergent SRLV E1 genotype was characterised in Italy as a low pathogenic caprine lentivirus, since it did not cause clinical symptoms in goats (6). Sequences of SRLVs isolated from Spanish sheep in an arthritis outbreak were assigned to subtype B2 (20), while proviral sequences isolated from Spanish sheep with neurological signs were assigned to genotype A2/A3 (4). Our results also revealed that strains representing subtypes A17, A20, A27 and B3 had a unique T to A substitution in the fifth position of the TATA box. The significance of this mutation is still unknown. An identical mutation (TATAAAA) was observed in subtypes E and J of HIV-1 strains (28). On the one hand, in vitro and in vivo transcription analyses indicated that mutation in the TATA box dramatically decreases HIV gene expression (28). On the other hand, Jeeninga et al. (7) revealed that mutation in the TATA box did not affect the activity of the subtype E LTR of HIV-1, which may indicate that the TATAAAA motif functions as an alternative TATA box. Additionally, it was also revealed that the sequences flanking the TATA box could contain multiple recognition sequences that may impart function to the TATA element (19).

In summary, the results of this study extend the current knowledge on the genetic diversity of SRLV field strains in Poland. Our results confirmed the existence of the A1, A5, A12, A13, A16–A18, A23, A24 and A27 subtypes and revealed that the phylogeny of a single genomic fragment does not necessarily reflect the phylogeny of the whole genome sequence because it reflects only a portion of it. Furthermore, we confirmed that mixed flocks favoured the emergence of new SRLV variants. Genetic analysis of LTR sequences revealed variability across SRLV subtypes. Several subtype-specific markers were identified; however, further studies should be conducted to assess whether these variations determine different transcription, viral fitness, and pathogenicity levels among distinct SRLV subtypes.