The tissue virome of black-spotted frogs reveals a diversity of uncharacterized viruses

Abstract

Amphibians are an essential class in the maintenance of global ecosystem equilibrium, but they face serious extinction risks driven by climate change and infectious diseases. Unfortunately, the virus diversity harbored by these creatures has been rarely investigated. By profiling the virus flora residing in different tissues of 100 farmed black-spotted frogs (Rana nigromaculata) using a combination of DNA and RNA viromic methods, we captured 28 high-quality viral sequences covering at least 11 viral families. Most of these sequences were remarkably divergent, adding at least 10 new species and 4 new genera within the families Orthomyxoviridae, Adenoviridae, Nodaviridae, Phenuiviridae, and Picornaviridae. We recovered five orthomyxovirus segments, with three distantly neighboring two Chinese fish-related viruses. The recombination event of frog virus 3 occurred among the frog and turtle strains. The relative abundance and molecular detection revealed different tissue tropisms of these viruses, with the orthomyxovirus and adenoviruses being enteric and probably also neurotropic, but the new astrovirus and picornavirus being hepatophilic. These results expand the spectrum of viruses harbored by anurans, highlighting the necessity to continuously monitor these viruses and to investigate the virus diversity in a broader area with more diverse amphibian species.

Introduction

Amphibians are a class of tetrapod vertebrates that are characterized by their ability to inhabit both aquatic and terrestrial environments. They are key components of various ecosystems throughout the world by serving as both the predator of invertebrates and the prey of other vertebrates (Zug and Duellman 2023). The creatures are globally distributed with exception of the polar areas, particularly concentrated in the neotropical countries (Wake and Koo 2018, Zug and Duellman 2023). They are comprised of three orders: Anura (frogs and toads), Caudata (newts and salamanders), and Gymnophiona (caecilians), each with a distinctive body type (Zug and Duellman 2023). Frogs and toads are the most abundant and diverse amphibians, accounting for ∼88% (n = 7678) of the total 8678 species (Frost 2024). Some anurans secrete substances on their skin that are of medical use, such as antibiotics, anesthetics, and painkillers (Zug and Duellman 2023). They are also known for their ecological and economical values by controlling pests that destroy crops or transmit diseases (Wake and Koo 2018). Especially, some frog species have been exploited as food, attracting a large scale of farming and consumption (Zug and Duellman 2023).

Although amphibians deeply participate in the energy circulation of the global ecosystem and hence are indispensable to maintaining the equilibrium of ecological spheres, they have confronted multifaceted threats that push many species to the edge of extinction (Luedtke et al. 2023). Indeed, amphibians are the most threatened vertebrates, with 40.7% of species having their populations in decline (Luedtke et al. 2023). Human activity and climate change are among the factors responsible for the status deterioration (Wake and Koo 2018). For example, deforestation and urbanization destroy their natural habitats, the wide use of pesticides poisons them via food chains, and the environmental pollution interferes with their reproduction and development (Wake and Koo 2018). From 2008 to 2016, the increasing air temperature and the simultaneous peak drought resulted in a 20% reduction of the California newt to mean body condition (Bucciarelli et al. 2020). Another detrimental factor is pathogens. The widely distributed chytrid fungus Batrachochytrium dendrobatidis has a devastating impact on amphibians worldwide (Wake and Koo 2018). The pathogenic Elizabethkingia miricola infection caused epidemic meningitis-like disease outbreaks in Chinese frogs (Hu et al. 2017). An introduced ranavirus led to the collapse of the amphibian communities in northern Spain (Price et al. 2014). These situations highlight that multiple conservation actions should be urgently planned and implemented to reverse the current trends (Luedtke et al. 2023).

Undoubtedly, investigation of the virus diversity harbored by amphibians is one of conservation measures, which helps recognize viruses with potential risks and take prevention measures in advance. Many viruses are associated with diseases in amphibians (Johnson and Wellehan 2005, Latney and Klaphake 2020). For example, ranid herpesvirus 1 is a causative agent of renal adenocarcinomas in leopard frogs (Rana pipiens) (Johnson and Wellehan 2005). The ranavirus infection can result in abnormal behaviors in frogs, such as lethargy, anorexia, buoyancy problems, etc. (Chinchar et al. 2014). Currently, the knowledge about amphibian viruses is unparallel to the diversity of amphibians. Thanks to the advance of high-throughput sequencing, powerful viromic techniques were occasionally applied to the amphibian virus investigation, which dramatically expanded their diversity (Reuter et al. 2018, Russo et al. 2018, Harding et al. 2022, Parry et al. 2023, Chen et al. 2023a), e.g. influenza-like viruses in Asiatic toads (Bufo gargarizans) (Shi et al. 2016), hepadnaviruses in Tibetan frogs (Nanorana parkeri) (Dill Jennifer et al. 2016), hepatitis E-like viruses in agile frogs (R. dalmatina) (Reuter et al. 2018), etc., indicating that current amphibian virus diversity only represents the tip of the iceberg as regards the virus flora present in these animals, with a vast majority of amphibian viruses awaiting being uncovered.

China has abundant amphibian resources, covering 665 species, with 567 species being anurans (Che et al. 2022). Some anuran species, e.g. R. nigromaculata, R. rugulosa, R. spinosa, etc., have been extensively farmed for food. In recent years, frog farming has been rapidly developed, and production has reached 210 084 tons in 2022 (Ministry of Agriculture and Rural Affairs Fisheries Administration of P. R. China, National Fisheries Technology Extension Center and China Society of Fisheries 2023). The frogs R. nigromaculata, a.k.a., black-spotted frogs (BSFs), are an indigenous species widely distributed in China, particularly concentrated in south China (Che et al. 2022). They favor residing in paddy fields, consuming a large number of diverse pests. They are also a favorite food for humans, and large-scale aquaculture has brought considerable profits to farmers (Xiong et al. 2022). However, due to high farming density, improper management, and other factors, infectious diseases bring a huge economic burden to farmers as well as pose substantial threats to the wild population (Zhong et al. 2018). The development of vaccines is one possible countermeasure; yet, it is currently hindered by the poorly described virus diversity harbored by these animals. In this study, we aimed to profile the complete virus flora among BSFs collected from five farms in Hunan and Jiangxi provinces using a combination of DNA-specific multiple displacement amplification (MDA) and RNA-specific meta-transcriptomic (MTT) viromic techniques. The results reveal multiple novel viruses from at least 11 families, with most vertebrate viruses divergent enough to be candidate for new genera, which greatly refreshes our knowledge of the virus diversity harbored by the species.

Results

Overview of frog whole virome

We sampled 100 BSFs from five farms in Changde, Changsha, Xiangtan, and Xinyu in 2022. Their intestine, lung, liver, brain, and spleen tissues were collected and mixed for viromic analysis according to farms, which resulted in 25 MDA and 25 MTT libraries for Illumina sequencing. The query of contigs against our eukaryotic viral reference database (EVRD) version 2.0 (Chen et al. 2022) resulted in 28 viral contigs, with 2406–104 462 base pairs (bp) in length. These viral contigs were annotated to at least 11 viral families (Fig. 1a), i.e. Solemoviridae (n = 1), Phenuiviridae (n = 1), Partitiviridae (n = 2), Orthomyxoviridae (n = 1), Nodaviridae (n = 3), Endornaviridae (n = 5), Dicistroviridae (n = 2), Picornaviridae (n = 2), Astroviridae (n = 3), Iridoviridae (n = 1), Adenoviridae (n = 3), and other unclassified viruses (n = 4) (Supplementary Table S1). Mapping reads against these viral contigs showed that the intestine samples had the most abundant and diverse viruses, followed by livers, lungs, brains, and spleens (Fig. 1a). Some viruses were prevalent across locations and tissue types (Fig. 1a); e.g. a turtle picornavirus 2 showed high abundance in all liver libraries and was also present in the intestines, lungs, and spleens of almost all sites, but was absent in brains. Interestingly, a plant-related partitivirus was prevalent in the brains, intestines, lungs, and livers of almost all sites (Fig. 1a).

We then compared the within-library diversity across sites and tissue types using Shannon and Simpson indices. At the site level, BSFs from Xiangtan had higher virus richness than from Xinyu (Wilcoxon test, P < .05) but the others did not show significant differences (Fig. 1b), suggesting geolocation is a minor factor impacting the BSF viromes. As shown in the heatmap, BSF viromes were significantly different in virus richness at the tissue type level, with the intestines having the richest virus species composition, while the spleens came in last (Fig. 1b). The analysis of between-library diversity showed that these libraries were largely aggregated according to the tissue types, and libraries of these solid tissues, i.e. livers, lungs, spleens, and brains, showed more similar virome composition within the same tissue types as compared to that of intestine libraries (Fig. 1c). These should be attributed to the fact that some viruses showed specific tropisms for certain tissue types and that the intestine viromes were more susceptible to external environmental factors (Chen et al., 2023b).

Characteristics and phylogeny of a novel orthomyxovirus

Orthomyxoviruses are a group of single-stranded negative-sense RNA (-ssRNA) viruses, with 6–8 segments as their genome (ICTV 2023). They gain frequent public health attention due to viruses within genera Alphainfluenzavirus, Betainfluenzavirus, Thogotovirus, and Quaranjavirus capable of causing diseases in humans (Shaw and Palese 2013). Our MTT viromics revealed a contig from the intestine library of Changde BSFs (JKCDC22) showing ≤33.9% amino acid (aa) identities over the entire coding region with the segments polymerase basic 1 (PB1) of Influenza C viruses (CIVs), indicating the discovery of a novel orthomyxovirus. To recover the remaining segments, we initiated an iterative PSI-BLAST search of all orthomyxovirus protein sequences against contigs, which confined orthomyxovirus-like candidates to 2689 contigs. We then examined their sequencing depth, length in the range of 800–3000 nucleotides (nt), presence or not of inverted complementary extreme ends, and presence or not in the same libraries, which was further followed by manual verification using HHblits search against UniRef30 version 2023_02. As a result, we detected four additional orthomyxovirus segments from the library JKCDC22.

The five orthomyxovirus segments were 1943–2778 nt in length (Fig. 2a), with a similar relative abundance of 59.1 ± 11.2 (mean ± standard deviation; assessed using mapped reads per kilobase of the viral contig length per million unclassified reads, FKPM) in the library JKCDC22, suggesting that these contigs come from a same virus. They had conserved extreme ends, i.e. 5ʹ-AAAANNNCCU-3ʹ, which is the same as those of influenza viruses (Shaw and Palese 2013). Like other orthomyxoviruses, the five segments encompass a single open-reading frame (ORF) (Fig. 2a). The HHblits analysis showed that the five segments hit against the influenza A virus (AIV) PB1, the Hubei earwig virus 1 PB2, the AIV PB3, the influenza B virus (BIV) hemagglutinin (HE), and the Wuhan Asiatic toad influenza virus (WATIV) nucleoprotein (NP) with probabilities of 100%, 98.8%, 99.1%, 96.6%, and 92.5%, respectively. Particularly, the NP-like sequence showed only 18.5% identity with an e-value of 0.27, suggesting that it was very divergent from the current known viruses. We then compared the five conserved RNA-dependent RNA polymerase (RdRp) motifs with those of other orthomyxoviruses. Motifs A, B, C, D, and E were sequentially arranged in the segment PB1, with the same conserved residues as those of other orthomyxoviruses (Fig. 2a).

The representatives of all 21 orthomyxovirus species, along with other unclassified anuran- and fish-related viruses, were used to infer the phylogenetic relationships of the five segments. The maximum likelihood (ML) trees showed that these anuran- and fish-related viruses occupied multiple basal lineages with similar topologies in their PB2, PB3, and NP (Fig. 2b–f). Interestingly, segments PB1, PB2, and PB3 of the virus did not neighbor the other anuran-related viruses (Fig. 2b–d) but were next to a sturgeon influenza-like virus (SILV) and a grass carp influenza-like virus (GCILV). The two fish-related viruses were recovered from transcriptomes of sturgeons and grass carps in China, with their remaining segments unavailable (Petrone et al. 2023). The segment HA was positioned next to AIV, BIV, and WATIV (Fig. 2e), and the segment NP was rooted at the base of fish-related viruses (Fig. 2f). But the phylogenetic positions of HA and NP were not well supported by the SH-aLRT and bootstrap tests (Fig. 2e and f), possibly due to a lack of the counterparts of SILV and GCILV. We calculated the pairwise identities of the five segments between these orthomyxoviruses (Fig. 3a–e). The virus identified here showed the highest 39.9%–56.3% identities with SILV and GCILV in the segments PB1, PB2, and PB3, much lower than the average pairwise identities of different species within the genera Quaranjavirus and Thogotovirus (43.8%–76.1%) but close to the highest between-genus identities (39.2%–60.7%) within the family (Fig. 3a–c). Accordingly, the virus is a candidate for a new genus within the family Orthomyxoviridae. Here, we coined “Chadfomvirus” for the genus name and “Chadfomvirus changdense” for the species name, reflecting a new orthomyxovirus identified in a frog of Changde, and the virus is the prototype of chadfomviruses (CFVs) within the genus.

We designed a nested RT-PCR method targeting the segment PB1 and used it to detect the virus in these tissue samples. The viromics captured CFV only in the library JKCDC22. The RT-PCR detection confirmed that the intestine pool contained two individual tissues (JKCDC22-10 and JKCDC22-19, positive rate: 10.0%) positive for CFV (Table 1). Additionally, one brain sample (JKCDN22-4) of Changde BSFs was also positive for the virus (Table 1). We tried to recover more segments by deep-sequencing the intestine sample JKCDC22-10, but besides the five segments, we did not find any other orthomyxovirus-like sequences. We conducted a transmission electron microscopy (TEM) examination of the sample JKCDC22-10 and, surprisingly, found multiple virions within different visions. They were irregularly spherical with a diameter of ∼80 nm and were covered with ∼20 nm membrane-anchored peplomers (Fig. 3f). These virions were smaller than the typical size of AIVs (∼100 nm in diameter), did not possess AIV-specific ribonucleoproteins in the interior of each particle (Shaw and Palese 2013), and shared few characteristics of the cell-isolated fish-related isavirus (ISV), mykissvirus (MKV), and sardinovirus (SDV) (Falk et al. 1997, Batts et al. 2017, Mohr et al. 2020). Viromic analysis of the intestine sample did not reveal any other viruses but CFV. Molecular detection of all viruses also revealed that the sample was only positive for CFV but not for other viruses, supporting the idea that these virions were CFV particles.

Characteristics of a recombinant frog virus 3

Frog virus 3 (FV3) is the prototype member of the genus Ranavirus within the family Iridoviridae (Chinchar et al. 2014). It has a ∼105 kb-long linear, double-stranded DNA (dsDNA) genome (ICTV 2023). FV3 was initially isolated from a renal tumor of a leopard frog, but its association with diseases has not been well established (Johnson and Wellehan 2005). The MDA viromics here successfully recovered a complete FV3 genome from the brain library (JKCDN22) of Changde BSFs. We named the virus FV3/CHN/2022/Frog/JKCDN22 (abbreviated as FV3/JKCDN22). The genome was 104 462 bp in length and showed a similar organization as that of the FV3 prototype virus (Fig. 4a). It encoded at least 80 genes, with 68 associated with replication, transcription, and metabolism, while the functions of the remaining genes cannot be assigned (Fig. 4a). At the complete genome level, FV3/JKCDN22 showed the highest nt identity (99.1%) with FV3/CHN/Frog/MWH421017 (Genbank accession number: MG791866), followed by 98.8% with FV3/CHN/Turtle/SSTIV (accession number: EU627010). The former was isolated from epithelioma papilloma cyprini cells of a Chinese BSF, and the latter was detected in a soft-shelled turtle (Trionyx sinensis) from Shenzhen, Guangdong Province, China. We used FV3/JKCDN22 as the query sequence to examine any recombination events against its three closest relatives. Results showed that FV3/JKCDN22 was a recombinant with its major parent being FV3/MWH421017 and the minor parent being the turtle-originated FV3/SSTIV (Fig. 4b). The recombination occurred in the range of 80 750 bp–83 250 bp corresponding to a putative NTPase/helicase-like gene and a hypothetical gene (Fig. 4b). The phylogenetic analyses based on the DNA polymerase (DNAp) and NTPase/helicase-like genes showed that its closest relative was FV3/MWH421017 as to the gene DNAp (Fig. 4c). But when using the NTPase/helicase-like gene, it switched to FV3/USA/Frog and FV3/CAN/Frog/KEN1, while FV3/MWH421017 was clustered with the turtle strain and another Chinese frog strain (Fig. 4c), which, again, validated that recombination occurred among Chinese frog and turtle FV3 strains. The MDA viromics showed that FV3/JKCDN22 was most abundant in the brain and lung libraries of Changde BSFs with reads per million total reads (RPMt) of 10^3.0−3.2 and per unclassified reads (RPMu) of 10^4.5−4.7 though also present in other libraries (RPMt ≤ 10^0.8 and RPMu ≤ 10^2.5) (Fig. 1a). Meanwhile, the MTT viromics detected the virus in the brain and lung libraries of Changde BSFs with RPMt of 10^0.5−1.5 and RPMu of 10^2.1−3.1 (Fig. 1a). These suggested that FV3/JKCDN22 was replicating in the lung and brain at sampling. The PCR detection validated that FV3/JKCDN22 was more prevalent in the brains (25.0%) and lungs (15.0%) of Changde BSFs (Table 1). However, it was also detected in the livers and spleens of Changde BSFs and various tissues in Changsha2 and Xiangtan BSFs (Table 1).

Other divergent viruses

Besides above identified orthomyxoviruses and FV3, the viromes revealed several plant/fungi-associated viruses, such as solemoviruses, partitiviruses, and endornaviruses (Supplementary Table S1), which were not our focus and hence were not described in this study. The remaining contigs were affiliated with adenoviruses (AdVs), nodaviruses (NVs), phenuiviruses (PhVs), astroviruses (AstVs), and picornaviruses (PVs), which infect vertebrates and/or arthropods. These contigs were 3468–9996 bp-long, covering almost the complete genomes except the three AdV contigs (Fig. 5a–e). Many families of vertebrate viruses have their ultimate origins in invertebrates, and multiple arthropod-associated viruses can infect and replicate in vertebrate hosts (Li et al. 2015, Chen et al. 2023b); therefore, we identified and investigated these viruses to understand their diversities, tropisms, and prevalence. The three AdV contigs were detected in the intestine libraries of Changsha1 (JKCSC22) and Changsha2 (WTCSC22) BSFs and the brain library (JKCDN22) of Changde BSFs. The contig AdV/JKCSC22 covered from the 5ʹ end to almost the entire DNAp gene, while the contigs AdV/JKCDN22 and AdV/WTCSC22 detached from AdV/JKCSC22 and overlapped at the Hexon gene (Fig. 5a). All three contigs were very divergent from currently known viruses, with AdV/JKCSC22 showing 41.3% aa identity with a siadenovirus identified from a psittacine. The latter two were 69.4% aa identical to each other and showed 40.6%–41.3% aa identities with a Pacific parrotlet-related siadenovirus. Phylogenetic analyses based on the DNAp and Hexon revealed that AdV/WTCSC22 and AdV/JKCDN22 formed a lineage, which, as well as AdV/JKCSC22, was basally rooted next to all other AdV representatives (Fig. 5a). These suggested that the three AdV contigs represented at least two new species of a new genus within the family Adenoviridae (Table 2).

NVs are classified into two genera, i.e. Alphanodavirus and Betanodavirus (ICTV 2023). Their genomes are composed of two single-stranded positive-sense RNA (+ssRNA) segments, which encode the RdRp and the capsid protein, respectively (ICTV 2023). We obtained two NV contigs from the intestine libraries of Changde and Xiangtan BSFs. However, interestingly, they concatenated the two segments (Fig. 5b), which were confirmed by RT-PCR amplification and Sanger sequencing. They were only 30.0% aa identical to each other at the RdRp and ≤37.8% with currently known NVs. ICTV proposes the demarcation of genera within the family Nodaviridae based on host range, phylogeny, etc. (ICTV 2023). The ML tree showed that the two contigs were distantly clustered together and were positioned between the two genera (Fig. 5b). Therefore, they represent two species within a new genus of the family Nodaviridae. PhVs include several pathogenic viruses, such as Rift Valley fever virus, SFTS phlebovirus (ICTV 2023). Their genomes contain three -ssRNA segments, with the L segment being ∼6.4 kb in length (ICTV 2023). We captured only one complete L segment from the intestine library of Changde BSFs. It surpassed the maximum of all known PhVs with a length of 7295 nt and encoded the entire RdRp (Fig. 5c). It was 76.8% aa identical to an Asiatic toad phenui-like virus, which was detected in an Asiatic toad from Sichuan Province, China. It was clustered with the Asiatic toad phenui-like virus, but was distantly related to members of the genus Phlebovirus with the highest aa identity of 35.0% (Fig. 5c). Based on its phylogenetic relationship and host range (ICTV 2023), the virus can be considered at least a new species within the family Phenuiviridae.

AstVs are a group of +ssRNA viruses that have a broad host spectrum and infect a wide range of vertebrates (ICTV 2023). AstVs are comprised of two genera, i.e. Avastrovirus and Mamastrovirus, but there is a myriad of divergent AstVs that can be candidates for new species, even new genera (ICTV 2023). We recovered two AstV complete genomes from the liver library (WTCSG22) of Changsha2 BSFs and the lung library (JKXYF22) of Xinyu BSFs. AstV/WTCSG22 showed the highest aa identity of 53.6% with Zhengjiang Chinese fire belly newt AstV 2, while AstV/JKYXF22 was 62.0% identical to a Beihai tree frog astrovirus. They were 30.3% aa identical to each other over ORF1b, but both formed independent lineages distantly related to the known aviastroviruses and mamastroviruses (Fig. 5d). ICTV proposes species criterion based on host range and less than 65% aa identity (ICTV 2023). Hence, the two AstVs are candidates for two new species within the family. Members within the order Picornavirales (PVs) share a common origin in their RdRp genes (Wolf et al. 2018). They show diverse genomic organization and have been classified into nine families, of which the family Picornaviridae includes the most abundant species members, i.e. at least 158 species within 68 genera (ICTV 2023). We recovered four PV complete genomes, with PV/JKXTC22 and PV/WTCSC22 highly identical to a cripavirus and a rice curl dwarf-associated virus, respectively. Interestingly, the contig (PV/WTCSG22) identified from the liver library of Changsha2 BSFs showed 99.8% nt identity with the contig (PV/JKXTG22) from the liver library of Xiangtan BSFs, suggesting that they were two closely related variants of the same virus. However, they were very divergent from all other known PVs, with only 31.6% aa identity to a worm lizard PV within the family Picornaviridae (Fig. 5e), suggesting that they should be considered a member of a new genus within the family.

We screened all tissues to detect these viruses using specific PCR or RT-PCR methods (Table 1). Results were largely consistent with the viromic analyses, i.e. if a tissue type had a higher relative abundance of a virus, the virus was much more prevalent in the tissue. The molecular detection also showed that these viruses varied in their prevalence across these samples, with AstV/WTCSG22 positive in BSFs of all locations. Notably, PV/WTCSG22 was detected in all tissues of BSFs sampled from Hunan Province, showing an extremely high positive rate in the livers, e.g. 100% in Changde and 60% in Xiangdan, but only positive in five BSF brains sampled in Jiangxi Province, which might suggest the predominant circulation of the virus in Hunan BSFs.

Discussion

Here, we revealed, for the first time, the virus flora residing in different organs of BSFs using complementary viromic methods. Although MTT can capture a snapshot of the total transcriptome of all replicating viruses (Zhang et al. 2019), it is weak at detecting DNA viruses and unable to capture the complete genomes of DNA viruses (Sun et al. 2022). Through MTT and MDA sequencing, we assembled a number of high-quality viral sequences, including multiple complete genomes, which enabled us to better understand their genetic diversity and evolutionary position. Although we recovered only 28 high-quality viral sequences, they showed great genetic diversities. The phylogenetic analyses, coupled with pairwise identity comparison, demonstrated that some sequences are so divergent that at least 10 new species and 4 new genera should be established for them within the families Orthomyxoviridae, Adenoviridae, Nodaviridae, and Picornaviridae (Table 2). This should be partially attributed to the fact that amphibian viruses have been insufficiently investigated before. We expect that anurans will be an important growth point for virus diversity. Undoubtedly, exploration of anuran viruses will not only help anuran conservation but also expand our knowledge sets of virus diversity and evolution.

According to their host tropisms, relative abundance patterns across different organs, and molecular detection results, these viruses can be classified into two groups with different origins. One group was associated with the frog diet and contained these plant- and insect-related viruses, i.e. unclassified ribovirus, rice yellow mottle virus, unclassified picornavirus, uukuvirus, partitivirus, two endornaviruses, and two dicistroviruses (Fig. 1a). Most of them were only present in the intestine libraries with low relative abundances. However, we also noted that the rice yellow mottle virus and partitivirus were also present in libraries of liver, brain, or spleen samples. It should be attributed to the environment-originated contamination at sampling. This issue is hard to avoid in viromic analyses, and some studies have also found partitivirus- and dicistrovirus-like sequences in the blood viromes of humans (Phan et al. 2018, Cordey et al. 2019, Jurasz et al. 2021). The other group included the remaining viruses, which are capable of infecting vertebrates. They should be considered genuine BSF-infecting agents, since the viromics and the molecular detection both revealed that they were positive in multiple organs, indicative of the replication and spread of these vertebrate viruses in multiple organs of BSFs.

The most interesting virus identified here is the new frog orthomyxovirus. The virus is phylogenetically related to two Chinese fish orthomyxoviruses, plus these BSFs consumed commercial expanded pellet diets made of soybeans, fish protein concentrates, yeast extracts, etc., raising a question about whether the virus is a frog virus or a fish one. We tend to believe that the virus is a genuine frog-origin agent for three reasons. First, these fish components have been subjected to several procedures of industrial processing, such as pulverization, dehydration, and fermentation, which would destroy viruses in fish components. Second, if the diet was the source of the virus, the CFV should be common in these BSF populations. On the contrary, the virus was only found in a single pool, and the molecular detection validated that only 2 out of the 20 intestine samples were positive for it (Table 1). Last but not least, the virus showed a high relative abundance in the virome (FPKM: 59.1 ± 11.2) and was detected in the brain of a BSF, indicating it was replicating and had entered the brain at sampling. Currently, the members within the family Orthomyxoviridae have been rapidly expanded (ICTV 2023), with several anuran- and fish-related viruses awaiting approval by ICTV, which is largely attributed to the wide application of the viromic technique in virus exploration (Shi et al. 2018, Petrone et al. 2023). Orthomyxoviruses are among the most concerning pathogens because of their association with a number of emerging infectious diseases. Undeniably, to comprehend the origin, evolution, and cross-host transmission of orthomyxoviruses, it is imperative to monitor their diversity and host range.

Genomic recombination is an important evolutionary mechanism for viruses to respond to the pressure of natural selection and invade heterogenous hosts by producing a rich pool of genetic diversity (Martin et al. 2015). Co-infection of multiple viruses in a host provides ample opportunities for virus recombination. FV3 is a nucleo-cytoplasmic large DNA virus that is characterized by its replication strategy, i.e. viral DNA is transported to the cell nucleus to guide the synthesis of genome and greater than genome length DNA, then progeny DNA is transported into the cytoplasm, where large concatemers of viral DNA are formed by recombination (ICTV 2023). Although natural recombination of iridoviruses has been rarely reported, previous scattered studies showed that the recombination of an infectious spleen and kidney necrosis virus and a red seabream virus resulted in a recombinant megalocytivirus causing a red seabream disease in Taiwan (Shiu et al. 2018) and a Brazilian bullfrog FV3 recombined with an Australian toad Bohle iridovirus (Candido et al. 2019). Here, we identified a BSF FV3 that showed similarity profiles across the complete genome frequently intersecting with frog- and turtle-associated FV3s (Fig. 4b), indicating recombination occurred between the frog- and turtle-related FV3 viruses. FV3 is known for its wide host spectrum, ranging from anurans, reptiles, and fish (Chinchar et al. 2014). Cross-host transmission and recombination accelerate the evolution and mutation of FV3. FV3 was associated with a renal tumor (Johnson and Wellehan 2005), so we cannot underestimate that cross-host recombination would change the biological and pathogenic traits of FV3. This finding highlights that cross-host transmission and recombination of viruses are not a rare event among lower vertebrates and hence pose a substantial risk to their conservation. Accordingly, continuous surveillance of FV3 in anurans as well as in its other hosts should be intensified in the future, in which we not only need to understand its prevalence but also must sequence its complete genome to monitor its recombination events with other host species-related viruses.

This study profiled the viromic composition in different tissue types of BSFs, which provides an important insight into understanding the tissue tropisms of these frog viruses. This knowledge not only helps predict their pathogenicity but also helps deduce their transmission routes. We included the intestines, lungs, livers, brains, and spleens of BSFs in this study. Their viromes and molecular detections revealed that these viruses showed different tissue tropisms. Although the intestines harbored the most diverse viruses (Fig. 1a), some of them replicated mainly in other solid organs. CFV/JKCDC22 should be an enteric virus transmitted through the fecal and oral routes since it had a higher relative abundance in intestine samples (Fig. 4b) and was also more prevalent in these tissues. We also detected it in a brain sample, suggesting that the CFV may enter to and replicate in the brain. FV3/JKCDN22 replicated mainly in the lung and the brain, though it also appeared in the intestine, liver, and spleen (Fig. 1a). These AdVs, NVs, and the PhV should be enteric viruses, but AdV/JKCDN22 was also present in the brain tissues (Fig. 1a). Interestingly, we found that AstV/WTCSG22, PV/WTCSG22, and PV/JKXTG22 are hepatophilic (Fig. 1a). This is not surprising since some AstVs within the genus Avastrovirus and hepatoviruses and tremoviruses within the family Picornaviridae are capable of causing hepatitis in aves and/or mammals (ICTV 2023). However, we did not observe any clear association between these viruses and disease in the frogs. These frogs were apparently healthy at sampling, and the necropsy did not reveal any visible lesions in these organs. Even so, we cannot underestimate their impacts on frogs, particularly those in the wild. Nowadays, the environment and climate are undergoing dramatic changes, and frogs in the wild are under great pressure to survive. Undoubtedly, viral infections pose a substantial risk to their health. Therefore, further cell isolation and experimental animal infection with these viruses will be required to fully understand their virulence and pathogenicity. Besides, the virus diversity harbored by those frogs in the wild should be extensively investigated.

Materials and methods

Ethics statement and sample collection

The procedures for sampling and processing BSFs were reviewed and approved by the Administrative Committee on Animal Welfare of the Changchun Veterinary Research Institute (Institutional Animal Care and Use Committee Authorization, permit number AMMS-11-2020-012). We recruited BSF farms from Changde, Changsha, and Xiangtan in Hunan Province and from Xinyu in Jiangxi Province in 2022. Changsha had two farms, and the remaining locations had one farm for each. These farms were 20–150 km apart. These BSFs were fed with expanded pellet diets made of soybean and fish protein concentrate, and they can also freely prey on accidentally intruded insects. We randomly sampled 20 apparently healthy BSFs from each farm and collected their intestines, lungs, livers, spleens, and brains at the local laboratory, which were immediately cryo-transported to our laboratory, where they were stored at −80℃.

Sample processing and Illumina sequencing

To understand the complete spectrum of virus diversity harbored by these BSFs, samples of the same tissue type from each farm were mixed and subjected to virome analysis using a combination of MDA and MTT techniques (Sun et al. 2022). We cut about a match-head-sized piece (∼0.05 g) of each tissue for pooling and thoroughly homogenized them using a ground-glass tissue grinder (Chunbo-Bio, Haimen, China) with DMEM buffer at a v/w ratio of 1:5. Homogenates were clarified at 10 000 × g for 10 min at 4℃. The supernatants were filtered through a 0.45-μm pore-sized membrane (Millipore, Burlington, MA) to further remove cell and bacteria debris. Filtrates were digested using DNase I (TaKaRa, Dalian, China) and RNase A (TaKaRa) at 37℃ for 1 h to eliminate free nucleic acids. The resultants were allocated into two groups for MDA and MTT preprocessing, respectively. For MDA, the total DNA was extracted using a DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) and amplified using an Illustra GenomiPhi V2 DNA Amplification kit (GE healthcare, Chicago, IL). A total of 1 μg products were subjected to DNA library preparation using a NEBNext Ultra DNA library Prep Kit for Illumina (NEB, Ipswich, MA). The library was quantified using a Qsep1 Bio-Sequence Analyzer (Houze Bio-Tech, Hangzhou, China) and then sequenced on a NovaSeq 6000 sequencer (Illumina, San Diego, CA), with each library producing at least 6 Gb data. For MTT, the total RNA was purified using Trizol (TaKaRa) as per the manufacturer’s protocol. The rRNA depletion was implemented using a RiBo-Zero Magnetic Gold Kit (Epicentra Biotechnologies, Madison, WI). After quantification using a Qubit^TM 4 fluorometer (Invitrogen, Singapore), the remaining RNA was subjected to library preparation with a NEBNext Ultra directional RNA library prep kit (NEB). The library was sequenced as described for the MDA method.

Virome annotation

All raw data were quality checked using fastp version 0.20.0 with the default settings and then subjected to removal of host and microbiomic sequences. Since the genome for R. nigromaculata is not available, nor are its phylogenetic relatives, the resultant high-quality reads of each library were mapped against the COI and cytb sequences of R. nigromaculata using Bowtie2 version 2.4.1 with the very sensitive global preset. The unclassified reads were then subjected to a metagenomic assignment of bacteria, archaea, and fungi using Kraken2 version 2.0.9-beta with the RefSeq-based database version 9 March 2022. The remaining reads of each pool were de novo assembled using megahit version 1.2.9 with the default mode. Contigs of ≥500 bp were queried against our EVRD database version 2.0 (Chen et al. 2022) using DIAMOND blastx with an e value cutoff of 1e-5. These viral-like sequences (VLSs) were then searched against the RefSeq branches of fungi, archaea, and bacteria using blastn and DIAMOND blastx at the nt and aa levels, respectively. If a sequence had a hit to nonviral organisms with an identity of ≥50% (for both nt and aa) and an alignment length of ≥500 (for nt) or ≥150 (for aa), it was considered nonviral and was removed. The remaining VLSs were further queried against the EVRD problematic viral sequences that originate from backbones or functional cassettes of vectors, viral sequences related to reagents and laboratory components, and host genomes (Chen et al. 2022), with contigs showing ≥99% similarity over an alignment of ≥400 nt to subjects being removed. To recover viral sequences with remote homology to known viruses, VLSs were in silico-translated into protein sequences using prodigal version 2.6.3 with a standard translation table and the meta mode. The predicted protein sequences were subjected to iterative profile-profile searches against UniRef30 version 2023_02 using HHblits (Gabler et al. 2020). Hits were manually checked for their targets, probabilities, and e-values. To validate the quality of VLSs, we employed two means, i.e. a similar homology profile across a contig against its reference sequence and conserved or complementary sequences in non-coding regions. These unclassified reads of each clean data were mapped against viral contigs using Bowtie2 with the end-to-end sensitive mode. Mapped reads were counted using samtools version 1.10. Read counts of ≤10 were removed to minimize index hopping. The relative abundance of each viral sequence was calculated by dividing the million total and unclassified reads by the read count of the viral sequence.

Molecular detection and validation

We employed molecular detection to investigate the prevalence of these viruses and validate these sequences with unusual structures. We developed nested RT-PCR and PCR to detect RNA and DNA viruses, respectively. Primer pairs were designed according to contigs using Primer Premier 5 (Supplementary Table S2). Total RNA and DNA of each sample were extracted using a RNeasy Mini kit (Qiagen) and a DNeasy Blood and Tissue kit (Qiagen). RNA was reverse transcribed into cDNA using a first cDNA synthesis kit (TaKaRa). The cDNA and DNA were amplified using the 2 × PCR MasterMix (Tiangen) with the following PCR programs: 30 cycles (for outer PCR) or 35 cycles (for inner PCR) of denaturation at 94℃ for 30 s, annealing at 54℃ (or adjusted for certain primer pairs) for 30 s, and extending at 72℃ for 40 s. Distilled water was used as a negative control. The products were visualized by electrophoresis using 1.0% agarose gel. The expected amplicons were directly sequenced by the Sanger method on an ABI 3730 sequencer (Comatebio, Changchun, China).

Transmission electron microscopy

A CFV-positive intestine sample was subjected to TEM examination. The intestine sample was homogenized using DMEM and centrifuged at 10 000 × g for 10 min. Ten μL supernatants were negatively stained with 2% phosphotungstic acid and examined with a TEMH-7650 transmission electron microscope (Hitachi, Tokyo, Japan).

Phylogeny reconstruction and recombination prediction

Viral sequences were queried against Genbank using blastn (for sequences highly identical to references) or blastx (for divergent sequences) to choose their references. Sequences of prototype species or genus within certain families were also retrieved. Viral sequences were aligned against their references using mafft L-INS-I or E-INS-I (Katoh and Standley 2013), the resultant multiple sequence alignments were automatically trimmed using trimAl version 1.4 (Capella-Gutiérrez et al. 2009), and the evolutionary model was determined using modelfinder with the AIC criterion. The ML phylogeny was inferred using IQ-TREE version 1.6.12 with 1000 ultrafast bootstrap replicates (Nguyen et al. 2015). The pairwise identity between sequences was determined using SDT version 1.2. To detect the recombination event among FV3s, we first queried FV3/JKCDN22 against FV3/CHN/Frog/MWH421017, FV3/CHN/1995/Frog/RGIV, and FV3/CHN/Turtle/SSTIV using RDP4 (Martin et al. 2015), with RDP, GENECONV, BootScan, MaxChi, Chimaera, SiSican, 3Seq, LARD, and PhyIPro methods. If a recombinational prediction was simultaneously supported by at least three methods with MC-corrected probability of ≤1e-3, it was considered significant. The significant recombination event was exhibited using similarity plot and bootscan analyses using SimPlot version 3.5.1.

Acknowledgements

Not applicable.

Author contributions

B.H. and Y.H. conceived and designed the study. Y.H. participated in the sample collection. C.L., Y.H., N.L., and Y.L. performed the experiments. C.L., L.Y., and B.H. analyzed the data. B.H, C.L., C.T., and Y.H. draft the paper with input from all authors. B.H. and C.T. revised the paper.

Supplementary data

Supplementary data is available at VEVOLU Journal online.

Conflict of interest:

The authors declare that they have no conflicts of interest.

Funding

This research was supported by the National Key Research and Development Program of China (2023YFF1305401). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability

All Illumina data sets generated in this study have been deposited in the NCBI Sequence Read Archive (SRA) database under BioProject accession number PRJNA1071828. The viral contigs determined in this study are available in the NCBI Genbank under accession numbers PP352007-PP357223 and PP357265. All original underlying data in this manuscript, i.e. read count table, multiple sequence alignments, phylogenetic trees, and the original TEM photos, are available from Figshare at https://doi.org/10.6084/m9.figshare.25138856.v2.

References

Author notes