Complete genome and phylogenetic analysis of bovine papillomavirus type 15 in Southern Xinjiang dairy cow.

BACKGROUND: Bovine papilloma is a neoplastic disease caused by bovine papillomaviruses (BPVs), which were recently divided into 5 genera and at least 24 genotypes.
OBJECTIVES: The complete genome sequence of BPV type 15 (BPV Aks-02), a novel putative BPV type from skin samples from infected cows in Southern Xinjiang China, was determined by collecting warty lesions, followed by DNA extraction and amplicon sequencing.
METHODS: DNA was analyzed initially by polymerase chain reaction (PCR) using the degenerate primers FAP59 and FAP64. The complete genome sequences of the BPV Aks-02 were amplified by PCR using the amplification primers and sequencing primers. Sequence analysis and phylogenetic analysis were performed using bio-informatic software.
RESULTS: The nucleotide sequence of the L1 open reading frame (ORF) of BPV Aks-02 was 75% identity to the L1 ORF of BPV-9 reference strain from GenBank. The complete genome consisted of 7,189 base pairs (G + C content of 42.50%) that encoded 5 early (E8, E7, E1, E2, and E4) and 2 late (L1 and L2) genes. The E7 protein contained a consensus CX₂CX29CX₂C zinc-binding domain and a LxCxE motif. Among the different members of this group, the percentages of the complete genome and ORFs (including 5 early and 2 late ORFs) sequence identity of BPV Aks-02 were closer to the genus Xipapillomavirus 1 of the Xipapillomavirus genus. Phylogenetic analysis and sequence similarities based on the L1 ORF of BPV Aks-02 revealed the same cluster.
CONCLUSIONS: The results suggest that BPV type (BPV Aks-02) clustered with members of the Xipapillomavirus genus as BPV 15 and were closely related to Xipapillomavirus 1.

INTRODUCTION

Papillomaviruses (PVs) are a diverse crescent group of viruses whose genomes comprise small non-enveloped and circular double-stranded DNA viruses [12]. PVs have been reported to cause infections in a large variety of amniote species [2]. In cattle, bovine papilloma, also known as a wart, is the most common skin tumor caused by bovine papillomavirus (BPV) [3]. The BPV types are found by broad-spectrum detection in different places of the animals' body, even in healthy skin [45]. Several studies have reported that the viral types of BPVs have limited relationships with the clinical status, type of herd, and age of animals [67]. This suggests that traditional clinical diagnostic techniques may not be effective methods to determine the BPV types.

Recently, culture-independent molecular techniques are used to detect the PV without a virus culture system [8]. Based on the sequence similarity of the highly conserved open reading frame [ORF] L1 of PVs, the BPV types have been characterized and divided into 5 genera: Deltapapillomavirus, Epsilonpapillomavirus, Xipapillomavirus, Dyoxipapillomavirus, and Dyokappapapillomavirus [91011]. Thus far, 24 genotypes have been reported [111213]. BPVs 1, 2, 13, and 14 are classified in Deltapapillomavirus genus; BPVs 5 and 8 are grouped in Epsilonpapillomavirus genus; the Xipapillomavirus genus comprise BPVs 3, 4, 6, 9, 10, 11, 12, 15, 17, 20, 23, and 24; BPV 7 is classed as a member of the Dyoxipapillomavirus genus; the Dyokappapapillomavirus genus includes BPVs 16 and 18. Other new BPV 19 and BPV 21 belong to an undefined genus [111213141516].

BPVs have been studied widely in China. Recently, several PVs have been detected and reported in Xinjiang, China [171819202122]. In 2014, the authors examined 11 warty lesion samples from a dairy farm in Southern Xinjiang, China. Although the L1 gene sequence (GenBank accession No. KM455050.1) was submitted to the GenBank database, the complete genome sequence of the novel BPV (named BPV 15) has not been characterized. Therefore, this study conducted a complete genome and phylogenetic analysis of BPV type 15 in Southern Xinjiang dairy cows.

MATERIALS AND METHODS

Sample collection

In 2014, 11 warty lesion samples were obtained from a dairy farm located in the region of Aksu in Southern Xinjiang, China. All these samples were administered to the laboratory and examined by L1 gene and phylogenetic analysis. This paper reports the novel BPV 15, which was isolated from the warty lesion sample collected from an infected cow. The sample was called Aks-02. The lesion was dark grey, 0.5–1 cm in diameter, and located on the neck region below the jaw. These samples were suspended in a 50% glycerol phosphate buffer and stored at −20°C.

DNA extraction

The frozen warty lesion specimens were homogenized. The total genomic DNA was extracted from warty lesions specimens using a TIANamp Genomic DNA Kit (TIANGEN BIOTECH, China) according to the manufacture's protocol. The extracted DNAs were dissolved in 50 μL TE buffer and stored at −20°C until used.

Polymerase chain reaction (PCR) assay

Papillomaviral DNA was initially analyzed by PCR using the degenerate primers, FAP59 (forward: 5′-TAACWGTIGGICAYCCWTATT-3′) and FAP64 (reverse: 5′-CCWATATCWVHCATITCICCATC -3′) [23]. The expected length was approximately 480 base pair (bp).

To characterize the complete genome sequences of the BPV Aks-02 further, amplification primers and sequencing primers were designed by DNAStar version 5.0 software according to available genomes in GenBank under the following accession number: AY300819 and HQ612180. The primers set for amplification of the complete BAPV3 genome was BAPV-3_F (forward: 5′-CAGTGACACCTATTCCAAGAGGTT-3′ and BAPV-3_R2 (reverse: 5′- GCATGGACCCTAAACAAGTGCAAC-3′). Supplementary Table 1 provides details of the sequencing primers. Both amplification and sequencing primers were synthesized by Sangon Biotech (China). The amplification of viral DNA by PCR was carried out based on the manufacturer's recommendations of the LA PCR Kit (TaKaRa Biomedical Technology Beijing, China). Briefly, in a total volume of 20 μL: 2 × LA Buffer 2.0 μL, dNTP (2.5 mM) 0.5 μL; Primer BAPV-3_F (10 mM) 0.4 μL; Primer BAPV-3_ R2 (10 mM) 0.4 μL; template DNA 1.0 μL; LATaq enzyme (5 U/μL 0.4 μL; ddH2O 15.3 μL, were mixed homogeneously and short spun. The expected length was approximately 7,200 bp.

The PCR was performed in a thermocycler (TECHNE, TC-1500, UK) using the following cycling profile: 10 min at 94°C, followed by 32 cycles of 30 sec at 95°C, 45 sec at 56°C and 80 sec at 72°C; a final extension step of 10 min at 72°C was performed. A 5 μL aliquot of the PCR products was loaded on 1.5% agarose gel in Tris-acetate EDTA buffer at constant voltage (90V) for approximately 45 min and visualized under ultraviolet light.

Sequencing and sequence analysis

All products were purified using AxyPrep DNA Gel Extraction kits (Axygen Biotechnology, China) according to the manufacturer's instruction. PCR amplification of the complete genomes was sequenced bi-directionally using the amplification prime and multiply sequencing primers utilizing the Applied Biosystems 3730XL DNA Analyzer (Applied Biosystems, USA). Data of sequencing was assembled using DNAStar version 5.0 software and submitted to Blastn search. The characteristics of the complete sequence of the PV, including predictions of putative ORF, molecular weight, motifs, and regulatory sequences, were predicted using DNAStar version 5.0 software.

Phylogenetic analysis

The complete genome sequences of the L1 gene of the 23 related BPV genotypes from GenBank were imported and aligned using Clustal W in DNAStar version 5.0 software. Phylogenetic trees were constructed from the alignment of L1 sequences using the Neighbor-joining method in MEGA 6.0 software.

Gene sequence accession number

The complete sequence of the novel BPV 15 was deposited in the GenBank database with the accession number KM983393.1. Supplementary Table 2 lists the GenBank accession numbers of the different genotype sequences used in sequence analysis and phylogenetic analysis.

RESULTS

Clinical and pathological findings

In 2014, 11 warty lesion samples were taken from nine 12-month-old and two 18-month-old cows in a dairy farm, including 1,268 cows in southern Xinjiang. Eleven cows had a normal temperature, respiration, and appetite. All these warty lesions, which were mostly distributed in the face and neck, were approximately 0.8 cm in diameter and 1 cm in vertical height, grey and hyperkeratotic epidermis. In terms of appearance, there was little difference between them according to clinical and pathological findings.

PCR sequencing

DNA extraction from papilloma samples was amplified successfully by the PCR using the degenerate primers (FAP59/FAP64). The expected length of the PCR product was approximately 480 bp. The negative control of PCR with Millipore water was the non-amplified product.

The BPV types identified and obtained GenBank accession number in this study are as follows: BPV 2 (KM455049.1, KP663622, KP271447.1) in 81.8% (9/11) of the lesions, BPV 14 (KP336743.1) in 0.9% (1/11), and BPV 15 (KM455050) in 0.9% (1/11). These sequences were analyzed with the NCBI database, and BLAST showed that a putative new BPV type was discovered during this investigation.

PCR, sequencing and sequence alignments

The present identity between the ORF L1 nucleotide sequences of the BPV Aks-02 and BPV 3, -4, -6, -9, -10, -11, and -24 ranged from 71.8% to 76%, respectively (Supplementary Fig. 1). The results showed that BPV Aks-02 and these related BPV types were of the same species, according to the difference between 71% and 89% based on the L1 nucleotide sequence [1112425].

Based on the results of sequence alignments analysis of the L1 gene, the complete genome sequences of the identified BPVAks-02 were amplified successfully and sequenced with 17 sequencing primers (Supplementary Table 1), resulting in amplicons, approximately 7,189 bp size (Supplementary Fig. 2). According to E1, E2, E7, L1, and L2 nucleotide sequence alignment between BPV Aks-02 and related BPV, those of BPV Aks-02 shared more than 65.3% similarity with 5 ORFs of Xipapillomavirus 1. The similarities between the complete genomes of BPV Aks-02 and those of Xipapillomavirus 1 were a maximum of 69.5%–72.5% identity within 5 different genera (Table 1). Overall, these results suggest that BPV Aks-02 should be classified in Xipapillomavirus 1 of the Xipapillomavirus genus.

Table 1

Percentage nucleotide similarity (%) between BPV Aks-02 and several related BPVs

Genera	Deltapapillomavirus	Xipapillomavirus			Epsilonpapillomavirus	Dyoxipapillomavirus	Dyokappapapillomavirus
		Xipapillomavirus 1	Xipapillomavirus 2	Xipapillomavirus 3
Complete genomes	49.7–50.5	69.5–72.5	63.0	62.9	49.5–49.9	52.1	49.8–51.6
E1	53.3–53.6	74.1–78.3	63.5	66.3	52.1–52.6	54.0	53.4–54.4
E2	46.6–47.9	66.9–71.3	58.3	62.0	47.1–47.4	49.5	50.0–50.2
E7	43.5–44.1	65.3–75.2	55.7	55.2	37.2–41.7	46.2	42.4–47.0
L1	56.6–58.5	70.7–75.1	69.1	68.3	56.0–56.2	58.8	56.2–58.4
L2	44.0–45.2	69.8–73.9	64.3	61.1	29.6–31.3	47.8	42.1–46.6

BPV, bovine papillomavirus.

The phylogenetic tree of the L1 ORF sequence was constructed with 32 BPVs using the Neighbor-joining method in MEGA 6.0 software. PVs are classified based on the conserved L1 nucleotide sequence. As evident from the phylogenetic tree, the BPV Aks-02 was related more closely to BPV 9 (AB331650) of the species Xipapillomavirus 1, classified in the Xipapillomavirus genus (Fig. 1). Therefore, phylogenetic analysis indicated that the BPV Aks-02 appeared in the same cluster as BPV 9 (AB331650) within the genus Xipapillomavirus 1.

Fig. 1

Phylogenetic tree based on the papillomavirus complete L1 gene nucleotide sequences.

The phylogenetic tree was inferred using the maximum likelihood method based on the complete L1 gene nucleotide sequence alignment of related BPVs and BPV Aks-02 (the BPV Aks-02 indicated with a black diamond). The tree was divided into the previously determined genera Xipapillomavirus, Epsilonpapillomavirus, Deltapapillomavirus, Dyokappapapillomavirus, Dyoxipapillomavirus, and unassigned papillomavirus genus (BPV 19, BPV 21). The scale bar indicates the number of nucleotide substitutions per site.

BPV, bovine papillomavirus.

Complete BPV 15 genome structure

The complete genome structure of the BPV Aks-02 was predicted to contain 7,189 bp with a 42.5% GC content and to consist of 5 early (E8, E7, E1, E2, E4) and 2 late (L1 and L2) ORFs. The long control region (LCR) was located between the early and late ORFs of the circular virus. The E4 ORF was embedded within E2. E7 and E1, and E1 and E2, partially overlapped. The BPV Aks-02 was a lack of characteristic E6 ORF [1]. Fig. 2 and Table 2 show the predicted ORFs and characteristics of the genome structures.

Fig. 2

Complete genome structure analysis of the bovine papillomavirus Aks-02. The boxes represent the early (E) and late (L) open reading frames. LCR encompassing nucleotides 1–473: ▲, E1 binding sites; ●, Polyadenylation signals; ▼, Nuclear factor-1; ◆, TATA box of the viral promoters; △, E2 binding sites.

LCR, long control region.

Table 2

Main features of the proposed bovine papillomavirus genome

ORF	Position	Length (nt)	Length (aa)	Molecular mass (Da)	Isoelectric point
LCR	1–473	473	-	-	-
E8	474–602	129	42	4,858.97	5.411
E7	808–1,110	303	100	11,140.66	3.938
E1	1,100–2,929	1,830	609	69,497.13	6.556
E2	2,871–4,085	1,215	404	45,700.32	9.723
E4	3,391–3,819	429	142	16,302.39	5.938
L2	4,099–5,673	1,506	501	57,318.08	7.786
L1	5,684–7,189	1,575	524	57,206.94	4.759

LCR, long control region.

The LCR

LCR is the non-coding region (NCR) or up-stream regulatory region (URR) that was located at nt 1–473 between the stop-codon of L1 and the start-codon of E8 in the BPV Aks-02. LCR of the BPV Aks-02 contained 1 E1BS (AACAAT at nt 28–31), 2 canonical E2BS (ACCGN4CGGT at nt 322–333, nt 410–421), and 1 non-canonical E2BS (ACC-N6-GGT at nt 395–406), respectively. E2BS is essential for activating or repressing the transcription of the viral genome [26]. One modified E1BS (TAACAA at nt 366–371) was presented in the LCR. It contained only 1 polyadenylation (polyA) site (AATAAA at nt 81–86) and 2 TATA boxes (TATAAAA at nt 261–267, nt 425–431) in the LCR. Nuclear factor 1 (NF-1) binding sites (TTGGCA at nt 113–118) were also located in the LCR.

BPV 15 early region

The early region of the BPV Aks-02 includes 5 ORFs, such as E8 (128 bp), E7 (302 bp), E1 (1,829 bp), E2 (1,214 bp), and E4 (428 bp). The BPV Aks-02 genome included a putative E8 gene in the early region [127]. An important element of the E8 protein was presented to be the conserved C-terminal amino acids LXGWD and repeated sequences TGTCAACTGT (nt 568–577). The BPV Aks-02 E8 ORF shared 65.1% identity of the nucleotide sequence with BPV3 E8 ORF (AF486184.1) [11].

The BPV Aks-02 E7 ORF encoded a 100 aa protein and contained conserved retinoblastoma tumor-suppressor protein-binding domain (pRbBD: LXCXE). Furthermore, 1 consensus zinc-binding domain (ZnBD: CX2CX29CX2C) was located in the C-terminal region. The CXXC motifs were separated by 29 aa residues.

The E1 ORF coded approximately 609 aa as the largest protein of BVP 15. The conserved ATP-binding site of ATP-dependent helicase (GX4GKS) is located in the C-terminal region of E1. Two modified E1BS (TAACAA) were found at 1,357–1,362 and nt 1,540–1,545 in the E1 ORF. There was a lack of a leucine-zipper domain (LX6LX6LX6LX6L) in the E2 protein of the BPV Aks-02 [28]. The Putative E4 ORF showed a start codon and nested within the E2 ORF, which showed 71.8% similarity to BPV 3 (AF486184.1).

BPV 15 late region

The late region of the BPV Aks-02 contained the major capsid protein L1 and minor capsid protein L2. Two E1BS (AACAAT at nt 5,439–5,444, nt 6,917–6,922) and 1 modified E1BS (TAACAA at nt 5,438–5,443) were found in the late region. For the late viral mRNA transcription, polyadenylation signals PolyA (AATAAA) were identified at nt 7,040–7,045 of the L1. In addition, polyadenylation signals PolyA (ATTAAA instead of AATAAA) was found at nt 4,093–4,098 between the stop-codon of the E2 and the start-codon of L2 ORF. One canonical E2BS (ACC-N6-GGT at nt 7,126–7,137) in the L1 gene and 1 non-canonical E2BS (ACC-N7-GGT at nt 4,308–4,320) in the L2 gene were found. NF-1 binding sites (TTGGCA at nt 5,278–5,283) were also found in the L2 ORF.

DISCUSSION

In this study, the detection of papillomaviral DNA from warty lesion samples by PCR revealed different BPV types in southern Xinjiang, including BPV 2, BPV 14, and BPV 15. The present study confirmed a putative new BPV type, which was previously unreported and referred to herein as BPV Aks-02, from a skin wart excised from one of these cows.

BPVs have been divided into different genotypes based on the conserved region within the L1 ORF of different genera of PVs. The different genera are defined when the identity in the L1 nucleotide sequence is less than 60%, whereas species within a genus share 60%–70% identity [29]. Compared to other BPV genera, the Xipapillomavirus genus included 3 PV species: Xipapillomavirus 1 (BPV 3, 4, 6, 9, 10, 11, 15, 23, and 24), Xipapillomavirus 2 (BPV 12), Xipapillomavirus 3 (BPV 20) [1114151627]. On the other hand, the other 2 new BPV types, namely, BPV 19 and BPV 21, are still an unclassified genus, which requires further research [111415]. This paper described the identification of the complete genome sequence of BPV type 15, a novel putative BPV type from skin samples of cows in Southern Xinjiang, China. The BPV Aks-02 can be classified as the genus Xipapillomavirus 1 of the Xipapillomavirus genus based on sequence comparisons and phylogenetic analyses.

Different genomic structures of the PV genomes may have different evolutionary histories and complicating classification. Genetically diverse PV genotypes show remarkable conserved genomic organization in a wide range of host species. The core ORFs E1, E2, L1, and L2 are present in all characterized PV genomes [230]. The LCR is located between the early and late ORFs of the circular virus. The LCR, which is known as the NCR or URR, includes the origin of replication (ori) and multiple transcription binding sites [31].

In this study, the complete genome structure of the BPV Aks-02 contained 5 early (E8, E7, E1, E2, E4) and 2 late (L1 and L2) ORFs. Nevertheless, the percentages of the complete genome and ORFs sequence identity among different members of this group were closer to the genus Xipapillomavirus 1 of Xipapillomavirus genus, as shown in Table 1.

The early region mainly encodes the non-structural viral proteins for virus replication, genome maintenance, and promotion of cell growth [21431]. The early region of the BPV Aks-02 includes 5 ORFs, such as E7, E1, E2, E4, and a putative E8 gene in the early region [127].

The E1 and E2 of the early region in all PV are necessary for genome transcription and replication [14]. In addition, the E1 ORF contains a cyclin interaction RXL motif (the native KRRLL recruitment motif), which is essential for initiating PV replication. This could be a potential target for developing therapeutic agents or vaccine development [1429].

Many sites were found in the early region of the complete genome. Most PVs possess some of the transcription factor binding sites in the LCR, such as E2BS, E1BS, and NF-1 [323334]. E2BS are essential for activating or repressing the transcription of the viral genome [26]. In this study, LCR of the BPV Aks-02 contained 1 E1BS, E2BS, polyadenylation (polyA) sites, and TATA boxes. NF-1 binding sites (TTGGCA) were also located in the LCR of BPV Aks-02. The pRbBD found in the E7 protein of most PVs members plays a vital role in interfering with the host cell cycle and the immortalization and transformation of the host cell [111435]. The BPV Aks-02 E7 ORF encoded a 100 aa protein and contained a conserved retinoblastoma tumor-suppressor protein-binding domain (pRbBD: LXCXE). The BPV contained 1 canonical E2BS in the L1 gene and 1 non-canonical E2BS in the L2 gene. NF-1 binding sites were also found in the L2 ORF. The L1 protein showed a high proportion of positively charged residues (K or R) in the C-terminal end [36].

In conclusion, these findings confirm the genetic diversity of BPVs in South Xinjiang, China. Based on the complete L1 gene nucleotide sequences identified, the complete genome sequences, and gene nucleotide sequences of main ORFs, it is suggested that BPV Aks-02 can be classified in the Xipapillomavirus 1 of Xipapillomavirus genus. According to Nucleotide BLAST of the complete genome of BPV 15, there is little information on the BPV 15 gene in the GenBank. Currently, BPV 15 is rarely reported and has much less epidemiological data available globally. New data from this study may contribute to an understanding of the classification, genomic evolution, and epidemiology of BPVs.