Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle.

Bovine papillomavirus (BPV) is associated with bovine papillomatosis, a disease that forms benign warts in epithelial tissues, as well as malignant lesions. Previous studies have detected a co-infection between BPV and other viruses, making it likely that these co-infections could influence disease progression. Therefore, this study aimed to identify and annotate viral genes in cutaneous papillomatous lesions of cattle. Sequences were obtained from the GEO database, and an RNA-seq computational pipeline was used to analyze 3 libraries from bovine papillomatous lesions. In total, 25 viral families were identified, including Poxviridae, Retroviridae, and Herpesviridae. All libraries shared similarities in the viruses and genes found. The viral genes shared similarities with BPV genes, especially for functions as virion entry pathway, malignant progression by apoptosis suppression and immune system control. Therefore, this study presents relevant data extending the current knowledge regarding the viral microbiome in BPV lesions and how other viruses could affect this disease.

Introduction

Bovine papillomavirus (BPV) is a double-stranded circular DNA virus which mainly infects epithelial tissue from the Bovidae subfamily, which comprises the genus Bos, known for domestic cattle. In normal conditions, BPV infection can be asymptomatic and the warts, when they exist, can easily regress to normality. However, depending on the virus type and environmental cofactors, such as the ingestion of bracken fern,[3,4] the lesions can progress into malignant warts.[2,5]

The evolution of these lesions into a cancerous state is an alarming sign to countries that rely on livestock. BPV has been associated with economic losses in several countries, such as Brazil, Italy, Iraq, the United States, Japan, and Germany.[5-9] In Brazil, the virus was found in at least 9 states,[2,10-20] and makes the commercialization of cattle more difficult due to its unpleasant sight, leather deprecation, and some degree of infertility.[13,16] Females will also suffer from mild to severe milk production incapacitation and/or painful milking.

The BPV genome is divided into 2 parts plus the upstream regulatory region (URR). These parts are named early, which contains 5 or 6 open reading frames (ORFs) (E1, E2, E4, E5, E6, and E7), and late which comprises the L1 and L2 capsid proteins. The early genes are responsible for the mediation of all stages of the replication cycle inside the host cell. E1, E2, and E4 show helicase function, transcription and replication of the viral genome, and viral transmission facilitator, respectively.[1,21-24] On the other hand, E5, E6, and E7 act by inducing errors in the cell cycle due to the obstruction of suppression mechanisms and stimulate cell proliferation and survival, as well as keratinocyte differentiation, thus being considered oncogenes.[25-27]

The infection cycle begins with a lesion of the cutaneous or mucous tissues that exposes the basal cell layer. Once inside the host cell, the virus replicates itself as a plasmid, making use of the host’s cellular machinery. The virus forces the host cells to re-enter their cell cycle in order to amplify viral DNA, which culminates in the formation of epithelial warts with a rough appearance and variable size.

Despite the increasing number of studies regarding BPV, there is a lot that has not yet been answered, especially when it comes to infection pathways and the reason why some virus types tend to induce cancerous lesions while others remain as benign warts. The functions of all of the viral proteins are not fully understood, and we have yet to discover the factors that may interfere with the infection, such as protein expression, their interaction with the host and possible co-infections.

In this context, opportunistic microorganisms have been reported to be found in association to BPV, such as Staphylococcus (the main agent of bovine mastitis). In humans with HPV infection, co-infections have been reported with Chlamydia in the cervix intraepithelial tissue, in pharynx and larynx tumors with the polyomavirus BK and Epstein-Barr virus, or even in gastric adenocarcinomas with the Epstein-Barr virus and Helicobater pylori. Additionally, the interaction between HPV, cytomegaloviruses, and Merkel cell polyomavirus has already been documented by Adnan Ali et al and Vazquez-Guillen et al.

Studies published by Yaguiu et al and Carvalho et al have already detected co-infections with different BPV types in the blood stream of a single animal. also found the presence of feline papillomavirus in bovine epithelial warts. In an older study conducted by Goldstein et al, human cytomegalovirus raised the carcinogenic potential of BPV in an in vitro assay. Bovine herpesvirus has also been related to multiple diseases in ruminants, such as vulvovaginitis, encephalitis, pneumonia, and mastitis, with the latter being a point of interest in co-infection studies since BPV can also infect the udder and teats of females. However, the role of other viruses, in co-infection with BPV, in the development of papillomatous lesions is still poorly understood. Therefore, this study aimed to identify and characterize viral genes in cutaneous papillomatous lesions of cattle using a metatranscriptomic approach.

Materials and Methods

RNA-seq data acquisition

RNA-seq data used in this study were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). These RNA sequences were generated in a previous study by Barreto et al, which focused on the identification of differentially expressed bovine genes associated with BPV infection. In short, the sequencing of 3 libraries has been performed, which was constructed from the epithelial tissue of 3 samples of bovine cutaneous papillomatous lesions. All of the data used in this study can be found in the GEO database, under the accession number GSE122853.

Identification of host sequences

The reads were assembled into contigs by using Trinity version 2.9.1+galaxy2. The sequences used were paired-end. The parameters remained as default. Singletons were discarded. The sequences were aligned and mapped to a reference bovine genome (UMD3.1 Bos taurus 8) using Bowtie2 version 2.3.4.3+. RNA sequences aligned to the bovine genome were discarded. All of these tasks were performed with the web-based Galaxy platform for data intensive biomedical research version 20.01 (https://usegalaxy.org).

Identification of viral RNA sequences

The viral sequences were separated from those mapped to the bovine genome using the BAM filter tool (version 0.5.9), and the files were converted to FASTA format using the BAM conversion tool (version 2.4.0.0), under the BAMTools kit, both of which are available in the Galaxy platform. After filtering, they were suitable to undergo a BLASTx search to identify the viral genes that were expressed in the samples. To do this, we used the standalone BLAST+ suite (version 2.7.1), provided by the National Center for Biotechnology Information (NCBI). The installers and source codes are available from the NCBI FTP page (ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.7.1/). The database of choice was SwissProt version 5 (accessed in March 2nd, 2021), which is compatible with NCBI taxonomy filter, thus allowing us to filter the search for viral proteins only through the taxon ID 10239 (Viruses). In addition, we set an e-value cut off of 1e-5 in order to obtain significant hits. In order to compare the taxonomic structure between samples, Fisher’s Exact Test was used.

In order to confirm the taxonomy of each sequence, we have performed a taxonomy binning analysis with MEGAN version 6.21.2. The idea of this analysis was to differentiate widely conserved regions from species-specific genes, thus confirming the viral taxon level for each sequence and enabling further functional annotation of their genes. In this analysis, the naive LCA algorithm with 100% coverage was used.

Functional annotation

The genes identified by the BLASTx searches and confirmed through MEGAN were functionally annotated using Blast2GO version 5.2, part of the OmicsBox version 1.3.11 package (https://www.biobam.com). Gene Ontology (GO) analysis was carried out to retrieve information regarding the molecular functions, biological processes and cellular components of the query sequences. The annotation was submitted to Fisher’s Exact Test with a P-value filter set to >.05 in order to retrieve only the most statistically significant annotations and avoid generic functions. The genes were also mapped to metabolic pathways using Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Protein-protein interactions were analyzed using the STRING Viruses Consortium version 10.5 database (http://viruses.string-db.org) and the software Cytoscape version 3.8.0, with its NetworkAnalyzer app version 4.4.8 was used to merge interactions into a single graph.

Results

Read quantification and transcript mapping

The complete GEO dataset had a total of 52 879 412 trimmed reads and it is made up of 3 files. Each file had a sum of 17 917 849, 17 165 233, and 17 796 330 reads, respectively. After the assembly stage, we obtained a total of 254 833 contigs, divided into 79 967, 89 967, and 84 899 per file. Following the mapping with the bovine genome, we obtained a total of 169 889 (66.66%) contigs mapped and 84 944 (33.33%) contigs unmapped. This means that we have retrieved 27 755 (32.67%), 29 424 (34.63%), and 27 765 (35.68%) unmapped contigs per file.

Analysis of local alignment with BLASTx

The output of the BLASTx algorithm was XML and text files containing a generic name for each sequence and their respective scores against the viral database. The files reached a total of 9146 (10.76%) successful hits, while the amount of no-hits reached 75 798 (89.23%). The identity scores showed a mean value of 36.12%, and 302 (3.30%) hits had an identity score above 70%.

We retrieved a total of 177 viruses overall that belonged to 25 different viral families. The 3 most frequent viral families in terms of different viruses found within them were Retroviridae (90.96%), Poxviridae (38.98%), and Herpesviridae (35.59%), followed by Baculoviridae family in fourth place (17.51%) (Online Resource 1). All other families displayed very low frequency, remaining below 10% frequency.

The binned taxonomic tree generated by MEGAN can be seen in Table 1, where they explicit the exact taxa found for all samples and their respective frequency. This frequency refers to the number of sequences associated with a certain clade. It is possible to notice that not all viruses retrieved after the BLASTx search are displayed because in case of ambiguity when trying to determine the exact taxon for a specific sequence, the algorithm classifies it by a higher clade.

Table 1.

Taxonomic binning arrangement of retrieved taxa for all samples. The number inside parentheses refers to the frequency in absolute values of contigs associated with each taxon. The total number of contigs used in this analysis was 84 944, while the number of no hits summed 77 621.

Superkingdom	Kingdom	Class	Order	Family	Subfamily	Genus	Undefined rank	Species
Viruses (32)	—	—	—	Baculoviridae (15)	—	Alphabaculovirus (32)	—	Autographa californica multiple nucleopolyhedrovirus (171)
								Orgyia pseudotsugata multiple polyhedrovirus (105)
			Herpesvirales (16)	Herpesviridae (20)	Alphaherpesvirinae (50)	Varicellovirus (15)		—
					Betaherpesvirinae (4)	—		Human betaherpesvirus 5 (14)
					Gammaherpesvirinae (19)	—		Alcelaphine gammaherpesvirus (69)
						Rhadinovirus (39)		Human gammaherpesvirus 8 (110)
			Caudovirales (36)	Siphoviridae (25)	—	—		Bacillus virus Spbeta (119)
			—	—	—	—		Sulfolobus islandicus filamentous virus (15)
			—	—	—	—		Microptilis demolitor bracovirus (139)
	Riboviria (21)		—	Retroviridae (1)	Orthoretrovirinae (16)	Alpharetrovirus (116)		Avian leukosis virus (342)
								Fujinami sarcoma virus (104)
							Unclassified alpharetrovirus (29)	Avian retrovirus IC10 (30)
							—	Avian sarcoma virus (97)
								UR2 sarcoma virus (13)
								Y73 sarcoma virus (93)
						Gammaretrovirus (123)		Feline leukemia virus (141)
								Hardy-Zuckerman feline sarcoma virus (96)
								Kirsten murine sarcoma virus (15)
								Moloney murine sarcoma virus (132)
								Murine leukemia virus (617)
							Unclassified gammaretrovirus (1)	Abelson murine leukemia virus (63)
	Bamfordvirae (4)	Megaviricetes (13)	—	Phycodnaviridae (1)	—	Chlorovirus (14)	—	Paramecium bursaria chlorella virus (20)
			—	—	—	—		Acanthamoeba polyphaga mimivirus (2609)
			Pimascovirales (9)	Iridoviridae (21)	—	—		Invertebrate iridescent virus 3 (52)
						Iridovirus (5)		Invertebrate iridescent virus 6 (145)
		—	—	—	—	—		African swine fever virus (125)
			—	Poxviridae (2)	Chordopoxvirinae (72)	Avipoxvirus (3)		Fowlpox virus (520)
						Leporipoxvirus (1)		Myxoma virus (67)
						Orthopoxvirus (34)		Camelpox virus (14)
								Cowpox virus (83)
								Ectromelia virus (48)
								Vaccinia virus (274)
								Variola virus (32)
								Orf virus (16)

Functional analysis of identified genes and pathway mapping

GO term enrichment analysis annotated a total of 2299 molecular functions, 1648 biological processes, and 450 cellular components for all samples (Online Resources 2, 3, and 4).

Figure 1 shows the 15 most frequent functions listed for all files and divided into biological processes, cellular components, and molecular functions. Most of the functions listed are related to any kind of binding function or as a protein/enzyme with catalytic activity and being mostly expressed in the host nucleus. The categories of binding or catalytic activity are the most representative molecular functions. Most of the identified genes act preferably in the host nucleus, plasma membrane, or cytoplasm, which is very typical for viruses. In addition, genes related to the viral structure were also identified. The biological processes were mainly related to the viral infection into the cell and its replication inside the host.

Figure 1.

Quantification, in individual percentages, of the viral gene functions annotated with Blast2GO. Only the 5 most frequent functions for each category are shown.

KEGG pathway mapping has retrieved a total of 33 metabolic maps. Online Resource 5 shows a list with the name, map codes, and overall quantification for the maps retrieved. Of the 19 maps, 15 belonged to the metabolism class, 2 to the organismal systems class, 1 to the genetic information processing class, and 1 to the human disease class. 26.31% of the maps were related to the metabolism of cofactors and vitamins, 26.31% to nucleotide or amino acid metabolism, 15.78% to the metabolism of other biomolecules, 10.52% to biodegradation or biosynthesis, 10.52% to the immune system, 5.26% to translation, and 5.26% to cancer. All of the described pathways could be directly or indirectly related to the development of papillomatous lesions, but special attention must be paid to those of nucleotide metabolism, translation, the immune system, and cancer.

Of the 177 viruses identified, only 17 had a description of their protein-protein interactions (PPI) on the STRING Viruses Consortium. Nevertheless, we have managed to merge together all proteins identified even with different hosts, which generated a graph with 114 nodes and 207 edges with 11 connected components (Figure 2).

Figure 2.

Protein-protein interaction network provided by STRING Viruses Consortium and made into a visual graph by Cytoscape between the viral proteins found (in red) and the host’s (in blue). According to STRING description, most of these viral proteins carry functions related to cellular proliferation, viral replication, suppression of host’s immune response and apoptosis blockade. (A) Interactions among Enterobacteria phage lambda proteins and its host proteins (Escherichia coli). (B) Interactions among Varicella-zoster virus proteins and its host proteins (Homo sapiens). (C) Interactions among Equine Herpesvirus proteins and its host proteins (Equus caballus). (D) Interactions among Escherichia phage T7 proteins and its host proteins (Escherichia coli). (E) Interactions among Human Herpesvirus type 8 proteins and its host proteins (Homo sapiens). (F) Interactions among Vaccinia Virus proteins and its host proteins (Bos taurus).

The predicted functional partners of the input proteins that could be related to the papillomatous infection are the genes C3, PPIF, PPIC, PPIB, C4A, and CASP8 from Vaccinia virus (VACV); DAXX, IRF1, v-cyclin, CDKN1B, CDKN1A, CDK6, and vIRF-1 from Human herpesvirus 8 type P (HHV8P); gB, USP7, US3, UL34, ICP4, and ICP0 from Human herpesvirus 1 (HHV1); SRC, UL18, and US28 from Human cytomegalovirus (HCMVM); TNFRSF10A, Gp5, 1 from Enterobacteria phage T7 (BPT7); YES1, SRC, and PTGS2 from Equine herpesvirus 2 (EHV2); MDV092 from Gallid herpesvirus 2 (GAHVM); 66, CDC25C, and PRKD1 from Varicella-zoster virus (VZVD); and CrmB and T for Variola virus (VAR67). These genes present functions of immune defense inhibition, the override of cellular transcription processes, control over the apoptotic cascade and cell growth by control of the cell cycle.

Discussion

The present study has used next generation sequencing data to acquire, treat, and present data from bovine cutaneous papillomatous lesions for the analysis of the possible role of other viral agents within the BPV pathogenic progression. This made the identification of viral RNA possible, with the further elucidation of their functions and interactions with host genes through metabolic pathways. These results are important to increase the knowledge of virus-host interactions associated with bovine papillomatosis.

After the annotation of expressed genes, a predominance of processes such as DNA regulation, virion attachment to host cells, protein phosphorylation, and immune response suppression was observed. According to Pagliarini et al, these constitute the typical characteristics of oncogenes, which are in close contact with the DNA and induce mutations in it, especially causing alterations in the cell cycle and the regulation of host apoptosis. Some of the BPV proteins are also able to induce such alterations, like E5, E6, and E7, since they must override the cell cycle control in order to replicate and produce typical epithelial warts. The fact that some genes found in this study are situated in either the viral membrane or the host nucleus, and their main molecular functions are DNA/RNA/ATP binding, corroborates this statement.

It is worth mentioning that the annotation process was not uniform, with the absence of consistent information in the different datasets used. Furthermore, some of them were very specific genes, with few or no information about their functions, or they belonged to less studied organisms, which also made it difficult to describe them at a metabolic level. Besides, not all BLASTx hits may result in actual biological signals. After the taxonomic binning reconstruction, it was noticeable that not all retrieved viruses were displayed on the tree, such as the case for the Human immunodeficiency virus (HIV) and BPV, for instance. This means that the sequences regarded as HIV were aggregated in a higher taxon, possibly due to homology to highly conserved regions and therefore should not be considered as a real signal. The same situation occurs with Bicaudaviridae, Papillomaviridae, Flaviviridae, Myoviridae, Autographiviridae, Coronaviridae, Malacoherpesviridae, Lavidaviridae, and Herelleviridae families, along with some species within Retroviridae, Poxviridae, Iridoviridae, Alloherpesviridae, and Herpesviridae families.

Some sequences were identified as poxviruses, such as Myxoma virus (MYXV), Camelpox virus (CMLV), Ectromelia virus (ECTV), Fowlpox virus (FWPV), Orf virus (ORFV), and Cowpox virus (CPXV), whose hosts are cows, rabbits, camels, mice, chickens, humans, caprine, and a variety of mammals, respectively. VACV, along with FWPV, have been identified in humans.[50,51] Also, VACV was also found in milk samples from buffaloes in Brazil, and belongs to the Orthopoxvirus genus along with CMLV, Horsepox virus (HPXV), Monkeypox virus (MPXV), ECTV, Rabbitpox virus (RPXV), and CPXV, known for its considerable importance to both human and veterinary health due to its zoonotic nature. In general, viruses with large DNAs like poxviruses, adenoviruses, asfarviruses, and herpesviruses tend to express pro-survival cellular factors in order to extend the lifespan of the infected cells, basically stopping natural apoptosis. The protein B-cell lymphoma 2 (Bcl-2), which is one of the main apoptosis regulators in humans, was also identified in the form of homologs in MYXV. MYXV has not yet been documented in cattle infections to date, but it is known to cause a rapidly lethal disease called myxomatosis that creates cutaneous fibromas in the site of infection, in a way that may be similar to BPV. In addition, this poxvirus is transmitted through a mosquito bite and, although it is species-specific, it is believed that the pathogen is going through an event of cross-species transmission.[54,55] Similarly, Sheeppox virus (SPPV) is known to affect goats and sheep, causing a variety of skin surface growths such as papules, vesicles, pustules, and/or crusts. Although not much is known about SPPV infection in bovines, it is important to mention that SPPV is a member of the Capripoxvirus genus, which is responsible for some of the most important diseases in ruminants, including Lumpy skin disease virus (LSDV), which is characterized by fever, enlarged lymph nodes, and skin nodules. It is important to point out that this study identified sequences that presents homology with genes from these viruses, but this do not necessarily implicate that they were actually present in the lesion. Also, we cannot exclude the possibility that other bovine viruses, with high similarity to these identified viruses, may be present at the lesions, instead. Thus, further studies are still needed to confirm which virus can really be in these bovine lesions.

Some of the Herpesviridae family viruses found in this study include Human betaherpesvirus 5, also called Human cytomegalovirus (HCMV), Alcelaphine gammaherpesvirus 1 (AIHV-1), and Human gammaherpesvirus 8, also known as Kaposi’s sarcoma associated herpesvirus (KSHV). They are part of the Herpesvirales order, which is known for inflicting a large group of mammals, being mostly present in the mucous epithelium.[58,59] They act by reprogramming the signaling pathways inducing cellular factors (STATs, for example) to support their replication. In this way, they establish a latent infection, where the BPV E7 protein may modulate the host apoptotic process. CMVs, in general, are known to be highly species-specific, but since they are part of the Betaherpesvirinae subfamily, they are very much related to other mammal CMVs, especially the ones that infect primates and rodents. Despite this fact, HCMV has been extensively studied regarding its cell tropism, and it is known that they can infect a wide range of tissues within the host, which is mainly mediated by its UL24 protein. This protein is also shared among other herpesviruses, including Bovine alphaherpesvirus 1 (BoHV-1), and it has been reported in milk samples of clinically healthy cows, which means that they can establish a latent infection. They are also highly present in the respiratory tract and genital organs, mostly leading to infertility, much like BPV. AIHV-1 induces malignant catarrhal fever in bovines with the growth of papules and erosive lesions in many different tissues, mostly in the digestive tract. Lastly, KSHV is one of the human oncoviruses that mainly infect immunosuppressed patients, such as the ones with AIDS or lymphomas. Its envelope glycoprotein, gp350/220, was predicted in other gammaherpesviruses, like BoHV-4 and AIHV-1, implying similarities in structure and function. All these viruses are also known for its opportunistic nature, inducing symptoms only when the host has its immune defenses lowered and are capable of maintaining lifelong latency in a way that is similar to BPV.

Most of the retroviruses retrieved infect a variety of vertebrates, like mice, felines, and specially birds. Although not targeting bovines specifically, these viruses displayed some proteins in common, especially the SRC gene. This gene and its homologs were the ones with the best annotations and were present in 36 KEGG metabolic maps, most of them related to bacterial infections (shiguellosis, tuberculosis, signaling in H. pylori pathway and bacterial invasion into host cells), hormone signaling pathways (relaxin, oxytocin, thyroid hormones, estrogen, and GnRH), viral infections (hepatitis B, human cytomegalovirus, and Kaposi sarcoma associated herpesvirus infection), cancer (viral carcinogenesis, bladder cancer, and proteoglycans in cancer), cellular functions (endocytosis, mitophagy in animals, platelet activation, and regulation of actin cytoskeleton), and cellular factor signaling (EGF-R, ErbB, MAPK, Rap1), among others. In addition, it is also present in the pathway for the PD-L1 expression and PD-1 checkpoint in cancer, where the super-expression of the epidermal growth factor (EGF) stimulates its receptor (EGF-R) and culminates in the expression of the programmed death-ligand 1 (PD-L1) that stimulates PD-1 and could interfere with cell cycle progression.

The majority of these functions are related to the promotion of host invasion by microorganisms, infection state management, immune suppression and cell proliferation and survival. In this context, SRC is a typical oncogene, which has a strong transforming capacity down to a genetic level. According to Hbibi et al, and Roskoski the SRC genes hold the capacity to cause a variety of cancers, including colorectal, ovarian, neck and head, lung and esophageal. In the first study, the authors reported that tumor progression and posterior metastasis is directly linked to the super expression of EGF-R, which interacts with tyrosine-kinase p60-SRC receptor to raise the DNA synthesis in fibroblasts. Thus, it is possible that this super-expression of EGF-R could also increase the synthesis of keratinocytes during a co-infection between BPV and a tumor emergent virus.[19,69] Furthermore Cheng et al, have shown that another bovine virus, the Bovine ephemeral fever virus (BEFV) (ie, a Rhabdovirus), triggers the SRC/JNKAP1 pathway in a similar way, leading to viral entry enhancement and posterior replication by the inhibition of PI3K/Akt/mTORC1 pathway. Additionally, due to its activities of membrane fusion and immune suppression, SRC genes could be the key to transforming benign warts in cancerous lesions.

With the alignments obtained after running BLASTx program, not only did we find viral families related to infections in vertebrates, but a variety of taxa that range from archaea and algae infections or even invertebrates. This could be attributed to one of the main aspects of metatranscriptomic analyses; since the samples are acquired from the environment, they are not pure and can be contaminated before, during and/or after the library assembly.[72-74] There are many explanations for the contamination, varying from the air itself, being a microenvironment rich in microbes; the skin of the animal is also a non-sterile tissue and may present latent viral particles or cells from other individuals due to contact.[72-75] Some viruses might even infect the host, allowing its transcripts to be identified bioinformatically, but they remain inactive or as commensals and refrain from developing symptoms, which is probably the case for the papillomaviruses, even though we have presumably excluded them after the mapping step. After all, some PV types could remain asymptomatic within the host, even in association with other pathogens.[19,76] Another explanation for this would be the possibility of cross-contamination. Kazemian et al assessed the possibility of cross-contamination of HPV38 in human samples, which was not expected to be found in endometrial cancer. The authors identified HPV38 reads in the same batch, with the same viral mutations and the low number of reads being considered evidence for this contamination. In our findings, although we have also identified low read number, we observed different types of BPV in different samples, which might demonstrate that they were not a consequence of cross-contamination. The same could be said to other viruses retrieved after BLAST, for instance HIV. It is very unlikely that HIV is infecting the bovine cells, and its presence may be explained by the sequence similarity and/or the homology that there is among them and other bovine viruses, such as Bovine leukemia virus (BLV). This similarity, mostly in small sequences, could interfere in assigning each read to a taxon. Therefore, a hit might present similarity with a human and a bovine virus at the same time, causing confusion in the taxonomic classification. However, the homology between them could be associated with similar protein functions, regardless of host. Although, we have used BLASTx for phylogenetic assignment, which is a standard method, it is known that highly conserved regions are less informative for taxonomic classification. In the taxonomic binning analysis, the sequences prior assigned as HIV in BLASTx were classified only as Orthoretrovirinae subfamily, which includes bovine retroviruses. Therefore, because of the homology among them, we cannot differentiate the viruses inside this subfamily in order to classify these sequences.

One limitation is that this study is a computational-based identification metatranscriptomic study based on experimental data previously deposited in genomic databases, and thus lacks experimental assays to confirm the presence of these viruses in the samples and discard the possibility of contamination. Contamination is an important issue; therefore, the studies should take into account the quality assurance methods, including reaction controls. Eisenhofer et al highlights the importance of 3 types of negative controls to monitor contaminants in processing the samples. In this context, it is important to point out that the experiments that gave rise to the RNA sequences utilized in this study used blank controls in the RNA extractions and in the amplification controls, but no sampling blank controls were used. Although extreme precautions were taken, contamination can occur in NGS data, mainly due to exposure to a long experimental pipeline. Therefore, it is prudent to take into account that some of the viruses identified might have their origin in contamination. However, it is also possible that these sequences may have originated from limitations of the computationally phylogenetic assignment of small RNA sequences using BLASTx, such as loss of information from non-global alignment, ignorance of population genetic and phylogenetic issues, and the use of sequence-level rather than clade-level confidence metrics.[78,81] To further comprehend the possible role of those viruses in the disease and actually confirm their presence as true biological signals, more studies are needed.

For this reason, the alignment scores and thresholds, as well as the taxonomic binning analysis, should be taken into consideration, to avoid taking into account organisms that were not effectively expressing genes in the animal tissue. For example, it was possible to identify sequences of viruses belonging to the Ascoviridae, Autographiviridae, Baculoviridae, Bicaudaviridae, Herelleviridae, Iridoviridae, Lavidaviridae, Lipothrixviridae, Mimiviridae, Myoviridae, Nodaviridae, Phycodnaviridae, Polydnaviridae, and Siphoviridae families, which are bacteria-, algae-, or even insect-specific viruses. However, these were the families with the lowest identity scores and mainly low frequency. They are part of a separate order of viruses, called nucleocytoplasmic DNA viruses (NCLDVs) or Megavirales, that comprise viruses with large or even giant virions and genomes (sometimes surpassing 500 kb). They share various common traits, like protein similarities and topological features, and are able to infect a variety of eukaryotes, ranging from mammals to invertebrates and protozoans. For instance, Baculoviridae is an insect-specific viral family, but has also been used as a vector for gene transmission; therefore, it can enter both insect and human cancer cells.[84,85] The same can be said for mimiviruses, since they were isolated from the stool samples of patients with pneumonia[84,85] and urine samples from a kidney-transplant recipient.

As for other viruses, despite being present in a large quantity in our samples, it is possible that this happened due to homology to mammalian sequences, since the identity scores of these viruses were very low; it is speculated that around 8% of the human genome is remnant of ancient endogenous viruses. Also, a recent Brazilian study has assessed the presence of bacteriophages in fecal samples of birds and small mammals in the cerrado biome, the same place where the Brazilian livestock is mostly held. Among their findings, the Smacoviridae family was detected and, although not much is known about its biology, it is believed that they infect archaea, which might explain the presence of viruses which infect thermophilic archaea in our results. In addition, the presence of viruses from the Lavidaviridae family might be explained through the simultaneous presence of the Mimiviridae family, since these viruses are the most commonly associated with virophages from Lavidaviridae.

Another way of explaining these viral families that are not commonly found in bovine infections could be their similarity to very conserved regions throughout the taxa. Taking into account the binned trees generate by MEGAN, some viruses (eg, HIV, HCoV, VZVD, BPV, and many more) were classified in a broader taxonomic rank, being aggregated in clades such as kingdom or phylum. This might be a clear signal that these viruses are not actually present in any sample, rather these sequences could be homologous to a certain conserved region among these viruses, as stated by Huson et al. In this context, we can say that we found a virus belonging to some taxonomic rank but we are not able to differentiate them.

It is worth mentioning that we performed an RNA-seq pipeline targeting the sequencing of mRNA from bovine samples, therefore identifying the mRNAs through their poly-A tails. Byrne et al reported the occurrence of polyadenylation in mimivirus transcripts, possibly leading to them being sequenced more than normal, which helps to explain why they appeared as the family with the highest number of transcripts in our results. This case is the same for some herpesviruses and some retroviruses.

The majority of viruses found did not target specifically bovines, rather they infected preferably mice, rats and model animals,[94-96] primates or humans,[97-99] and animals that may have contact with bovines in rural properties (ie, swines, birds, and felines) or with multi-specific viruses.[36,100-104] Some of these hosts are also farm animals and could possibly be in direct or indirect contact with the cows, facilitating viral transmission. In addition, these viruses share the same metabolic routes of infection and disease maintenance with the ones that are specific to bovines, sometimes displaying very similar symptoms. Viruses from the same families or subfamilies tend to share strong similarities to the others, even if they do not target the same host species.

It is important to point out that a large proportion of publicly-available RNA-seq data is mRNA specific, and thus may hinder meta-transcriptome studies. In order to further evaluate meta-transcriptome of viral lesions, a retro-transcription step previous to library preparation and amplification strategies may be required. These approaches, along with nuclease treatment, concentration, and viral nucleic acid purification steps, may be used in other studies to increase the available data on specific viral RNAs.

Conclusions

This study has made the identification of key genes that could be linked to the BPV infection and somehow alter its pathological manifestation possible, in most cases worsening the disease. It was also possible to notice that the genes characterized in this study share strong similarities with the ones from BPV, especially in terms of virion entry pathway, malignant progression by apoptosis suppression and immune system control, and the fact that possible co-infection with BPV and other viruses could interfere with the disease or worsen its scenario. Nevertheless, it is important to point out that further specific molecular detection experiments and functional studies are needed in order to experimentally detect and identify those viruses, since the methodology used analyzed short non-host sequences that could present high sequence similarity and/or homology with other possibly non-expected viruses, and to better understand their real relationship with bovine papillomatosis. However, this study presents relevant data, which extends the current knowledge regarding the viral microbiome in BPV lesions and how other viruses could affect this disease, highlighting the importance of NGS approaches to elucidate complex virus-host interactions.

Click here for additional data file.

Supplemental material, sj-jpg-1-evb-10.1177_11769343221083960 for Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle by Adriana O Fernandes, Gerlane S Barros and Marcus VA Batista in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-xls-2-evb-10.1177_11769343221083960 for Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle by Adriana O Fernandes, Gerlane S Barros and Marcus VA Batista in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-xls-3-evb-10.1177_11769343221083960 for Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle by Adriana O Fernandes, Gerlane S Barros and Marcus VA Batista in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-xls-4-evb-10.1177_11769343221083960 for Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle by Adriana O Fernandes, Gerlane S Barros and Marcus VA Batista in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-xls-5-evb-10.1177_11769343221083960 for Metatranscriptomics Analysis Reveals Diverse Viral RNA in Cutaneous Papillomatous Lesions of Cattle by Adriana O Fernandes, Gerlane S Barros and Marcus VA Batista in Evolutionary Bioinformatics