The diverse liver viromes of Australian geckos and skinks are dominated by hepaciviruses and picornaviruses and reflect host taxonomy and habitat

Abstract

Lizards have diverse ecologies and evolutionary histories, and represent a promising group to explore how hosts shape virome structure and virus evolution. Yet, little is known about the viromes of these animals. In Australia, squamates (lizards and snakes) comprise the most diverse order of vertebrates, and Australia hosts the highest diversity of lizards globally, with the greatest breadth of habitat use. We used meta-transcriptomic sequencing to determine the virome of nine co-distributed, tropical lizard species from three taxonomic families in Australia and analyzed these data to identify host traits associated with viral abundance and diversity. We show that lizards carry a large diversity of viruses, identifying more than thirty novel, highly divergent vertebrate-associated viruses. These viruses were from nine viral families, including several that contain well known pathogens, such as the Flaviviridae, Picornaviridae, Bornaviridae, Iridoviridae, and Rhabdoviridae. Members of the Flaviviridae were particularly abundant across species sampled here, largely belonging to the genus Hepacivirus: fourteen novel hepaciviruses were identified, broadening the known diversity of this group and better defining its evolution by uncovering new reptilian clades. The evolutionary histories of the viruses studied here frequently aligned with the biogeographic and phylogenetic histories of the hosts, indicating that exogenous viruses may help infer host evolutionary history if sampling is strategic and sampling density high enough. Notably, analysis of alpha and beta diversity revealed that virome composition and richness in the animals sampled here was shaped by host taxonomy and habitat. In sum, we identified a diverse range of reptile viruses that broadly contributes to our understanding of virus-host ecology and evolution.

Introduction

Animal virology has traditionally focused on viruses of mammals and birds as they are the most likely natural reservoir hosts for viruses that may emerge in humans and animals of economic importance (Zhang et al. 2018; Mollentze and Streicker 2020). While this has provided major insights, it has necessarily resulted in a skewed view of viral diversity that limits our understanding of virus evolution and ecology at a broad scale, including the frequency and determinants of cross-species transmission and host range. Indeed, many viral families traditionally associated with mammals and birds have now been found in other vertebrate classes such as amphibians, non-avian reptiles (hereafter referred to as reptiles), and fish (Shi et al. 2018; Harding et al. 2022). Since far greater biological diversity exists within these other vertebrates—that comprise ∼33,000 documented species—it is reasonable to assume that they also harbor a substantial diversity of viruses (Zhang et al. 2018).

A variety of factors make lizards particularly informative for the study of viral ecology and evolution. Squamates (lizards and snakes) are a highly diverse and successful group of vertebrates, comprising 96.3 per cent of reptile diversity (Pincheira-Donoso et al. 2013), with over 10,000 extant species globally (Herrera-Flores, Stubbs, and Benton 2021). Squamates exhibit a range of morphologies, life history traits, and diverse ecologies, inhabiting every continent except Antarctica (Pincheira-Donoso et al. 2013; Pyron, Burbrink, and Wiens 2013). Australia is home to the largest number of reptile species, and hosts ∼10 per cent of the world’s squamates (Tingley et al. 2019), which form the most diverse vertebrate order in Australia (Wilson and Swan 2017). These animals are well adapted to the Australian landscape, colonizing and radiating across a diverse range of environments and ecological niches (Pianka 1969, 1973; Morton and James 1988). Australian lizards comprise old endemic Gondwanan lineages that pre-date the isolation of Australia and Antarctica, as well as more recent immigrant lineages from the north (Brennan and Oliver 2017). For example, the Gondwanan group of Pygopodoidea have a crown age of greater than 50 million years, while the genus Gehyra have a crown age in the mid-Miocene ∼20 million years ago, having colonized from Asia (Oliver and Hugall 2017). Several squamate species are at high risk of extinction within the next 20 years (Geyle et al. 2021) and infectious disease emergence can severely threaten reptile populations (Zhang et al. 2018). Determining what viruses are circulating in Australian reptiles may therefore aid in the conservation of Australian reptiles. Despite this, and despite the importance of reptiles in ecology and evolution, reptilian viromes remain extremely under sampled, with virome research having been conducted sporadically on only a small number of species (Hyndman et al. 2018; Shi et al. 2018; Chang et al. 2020; Harding et al. 2022; Lu et al. 2022; Waller et al. 2022; Liu et al. 2023). Further, viruses have only been characterized/discovered in a limited number of reptile families, relative to mammalian, avian, and fish families, and within these sparsely sampled reptile families, the upper range of the number of known viruses is 10–100-fold smaller than in birds and mammals, respectively (Carlson et al. 2022).

A key, yet rarely addressed, aspect of virus evolution and emergence is understanding the traits that determine the diversity and abundance of viruses carried by a host (Wille et al. 2019). A number of studies have linked viral diversity and abundance to specific host traits, including phylogenetic history, habitat, body mass, geographic range, community diversity, biome, location, and infection status with particular pathogens (Olival et al. 2017; Wille et al. 2018; Geoghegan et al. 2021). However, the effects of host ecology on viral diversity have only been explored in a handful of systems with varying results. To better understand the diversity of viruses carried in lizards and how host traits affect viral abundance and diversity, we used meta-transcriptomic sequencing to explore the virome of nine Australian lizard species residing in various environments and habitats across the biologically diverse monsoonal tropics of northern Australia.

Methods

Ethics statement

All work was carried out according to the Australian Code for the Care and Use of Animals for Scientific Purposes with approval from the institutional animal ethics committee (Permit ANU animal ethics approval A2016-07) and State authorities (collecting permits NT 58454 and WA SF010911).

Animal sampling

Liver samples were collected from apparently healthy Carlia amax, Carlia munda, Carlia sexdentata, Cryptoblepharus metallicus, Heteronotia planiceps, Heteronotia binoei, Gehyra nana, Gehyra arnhemica, and Oedura marmorata. Lizards were collected in autumn and winter months between April 2016 and June 2017, in arid regions of the eastern Kimberley and mesic regions of the ‘Top End’, Australia (Table 1, Supplementary Table S1). Samples were collected in Australian bioregions Arnhem Coast (ARC), Victoria Bonaparte (VIB), Central Arnhem (CEA), Daly Basin (DAB), or Darwin Coastal (DAC) (Supplementary Table S1), as defined by Interim Biogeographic Regionalisation for Australia (IBRA), version 7. Excised liver was stored in RNAlater^™ Stabilization Solution (Thermofisher Scientific, MA, USA), at room temperature while in the field and then at 4ºC for longer term storage. Sampled animals displayed no obvious signs of serious pathology.

RNA extraction

Liver tissue was homogenized using the Qiagen Tissue Lyser II with 3 mm stainless steel beads in Qiagen buffer RLT (Hilden, Germany), and RNA extracted using the Qiagen RNeasy Plus minikit (Hilden, Germany) according to the manufacturer’s protocol. Purified RNA was pooled in equimolar ratios into 11 pools, grouping RNA samples from the same lizard species and collection location, with 6 to 12 individuals per pool (Table 1, Supplementary Table S1). Pooled RNA was further purified using the RNeasy MinElute clean-up kit (Qiagen, Hilden, Germany) and quantified using the Qubit RNA Broad-range Assay with the Qubit Fluorometer v3.0 (Thermofisher Scientific).

Meta-transcriptomic sequencing

RNA pools were assessed for quality using the Agilent 2100 Bioanalyzer with the Agilent RNA 6000 Nano Assay (Agilent Technologies, CA, USA). Library construction and sequencing was performed at the Australian Genomic Research Facility (Victoria, Australia). Libraries were constructed using the TruSeq Total RNA Library Preparation protocol (Illumina, CA, USA) following rRNA removal using the Illumina Ribo-zero Gold epidemiology kit. Paired-end (100 bp) sequencing of each RNA library was performed on a HiSeq 2500 sequencing platform (Illumina, CA, USA).

Genome/transcript assembly, annotation, and abundance calculation

After trimming with Trimmomatic v0.38 (Bolger, Lohse, and Usadel 2014), reads were assembled de novo using two separate assemblers—Trinity v2.5.1 (Grabherr et al. 2011) and Megahit v1.2.9 (Li et al. 2015)—to increase the chances of correctly assembling all viral contigs. Contigs from both methods were combined and duplicates removed (retaining the longest version) using CD-HIT-EST v4.8.1 (Fu et al. 2012). Viral contigs were identified using BLASTn (Altschul et al. 1990) and DIAMOND BLASTx (Buchfink, Reuter, and Drost 2021) tools through alignment with the NCBI nucleotide (nt) database (e-value cut-off 1 x 10^–10) and non-redundant protein (nr) database (e-value cut-off 1 x 10^–5), respectively. The Geneious assembler (available in Geneious Prime^® 2022.2.2) was used to extend viral contigs where possible. Open reading frames (ORF) were identified using the Find ORFs tool within Geneious Prime, and conserved domains were identified using RSP-TBLASTN v2.12.0+ (Altschul et al. 1997) against the NCBI Conserved Domains database. The abundance of each viral contig was calculated as expected counts (from mapped trimmed reads) using the RSEM tool in Trinity v2.5.2 (Li et al. 2010). Overall abundance was calculated as expected count/total number of trimmed reads in library × 100. Novel viruses that shared greater than 90 per cent amino acid identity in the RdRp were considered to represent the same species. Likely vertebrate-infecting viruses were defined as those that belong to classically vertebrate-infecting viral families and/or those that clustered with viruses known to infect vertebrate species in phylogenetic trees.

To reduce the reporting of false positives due to index-hopping during sequencing for each viral contig present in more than one library, a viral contig was presumed to be a contaminant from another library if it met the following criteria: contig abundance was less than 0.1 per cent of the abundance of that contig in the library where that contig was most abundant. This is based on the index-hopping rate of about 0.1–2 per cent as listed by Illumina (https://sapac.illumina.com/techniques/sequencing/ngs-library-prep/multiplexing/index-hopping.html).

Alpha and beta diversity analyses

Diversity statistics were only obtained for viruses considered likely to infect vertebrate species (i.e. those likely to infect the sampled host) and viruses likely to be exogenous. Retroviruses were excluded due to the difficulty in determining whether they are endogenous or exogenous, and viral groups with disrupted ORFs were considered likely to be endogenous. We performed generalized linear models and selected the best-fit model at the family taxonomic level (among a set of possible models describing the relationship between taxonomy, habitat, environment, and the number of individuals in the library) based on the lowest AIC, as done previously (Geoghegan et al. 2021; Chen et al. 2022). In addition to default family, we also assessed the best fit model with quasipoisson distributions due to small sample sizes, and report the P-values from these models. Due to the small same size of this study, results should be interpreted with caution. Host genus- and species-level taxonomy were not considered due to small sample sizes. The Csex_M library was excluded from all analyses because it contained no exogenous, biologically relevant viruses and as the only riparian library, may have skewed results. Omar_M was similarly removed as it was the only library representing Diplodactylidae. Alpha diversity, including richness, Shannon index, Simpson Index, Shannon effective, and Simpson effective were calculated for each library using a modified version of the Rhea scripts (Lagkouvardos et al. 2017). Beta diversity was calculated using the Bray–Curtis dissimilarity matrix and virome structure was plotted as a function of non-metric multidimensional scaling (NMDS) ordination and tested using Adonis tests (PERMANOVA) using the vegan and phyloseq packages (McMurdie, Holmes, and Watson 2013; Oksanen et al. 2022). Cmun_M was a clear outlier for the beta-diversity analyses and was thus removed. Analyses were performed in R version 4.2.3 in R Studio 2022.07.1 and scripts are supplied at https://github.com/michellewille2/LizardViromes.

Phylogenetic analysis

Viral amino acid sequences were aligned with representative sequences from the same viral family obtained from NCBI, using Clustal Omega v1.2.3 available in Geneious Prime. Where necessary, large data sets were condensed to a more manageable size using CD-HIT version 4.8.1 (Fu et al. 2012), using an identity cutoff of 95 per cent. TrimAl v1.4.1 (Capella-Gutiérrez, Silla-Martínez, and Gabaldón 2009) was used to remove ambiguously aligned regions (using the gappyout setting in all cases with the exception of the Flaviviridae alignment which required stricter manually applied settings to remove uninformative columns). Alignments were visualized in Geneious Prime. Maximum likelihood trees were inferred using IQTree v2.1.3 (Nguyen et al. 2015) with model selection estimated using ModelFinder within IQTree (Kalyaanamoorthy et al. 2017). Branch supports were estimated with the Shimodaira-Hasegawa (SH)-like approximate likelihood ratio test (Guindon et al. 2010). The size and length of each alignment are provided in Supplementary Table S2 and details of assembled viral nucleotide consensus sequences that were translated for inclusion in phylogenies are provided in Supplementary Table S3.

Virus nomenclature

Novel viral species were determined based on demarcation criteria assigned by the ICTV (International Committee on Taxonomy of Viruses) for the relevant families/genera (https://ictv.global/report/genome) and were randomly assigned provisional names. Note that because these virus names are provisional, we have only assigned common names rather than a full binomial nomenclature.

Results

Viruses in Australian lizard species

We sequenced eleven pooled RNA libraries with an average read count of 35,041,464 (Table 1). Each library was generated from liver samples of six to twelve apparently healthy lizards of a single species sampled from either the mesic Top End or from the relatively arid east Kimberley regions of Australia (Victoria-Bonaparte region; Table 1). Two widespread gecko species—Heteronotia binoei and Gehyra nana—were sampled from both regions. A high diversity of vertebrate-associated viruses was detected, comprising nine different viral families, twelve genera, and thirty-one novel viruses (Fig. 1). While for most viral families only a single genus and species were documented, we identified viruses from at least two genera in each of the Flaviviridae, Picornaviridae, Rhabdoviridae, and Amnoonviridae, and two species of the same genera in the Astroviridae (Fig. 1B). Despite being present in only one library, and being represented by a single viral species, the Arenaviridae were by far the most abundant viral family (Fig. 1), with the Flaviviridae the second most abundant and detected in all three lizard families (5/11 libraries) (Fig. 1A). Members of the Picornaviridae were also abundant and detected in five of eleven libraries, although were restricted to the Gekkota species (Gekkonidae and Diplodactylidae). No vertebrate-infecting viruses were shared between different lizard species (Fig. 1B). At higher taxonomic scales, the genera Hepacivirus (Flaviviridae), Shanbavirus (Picornaviridae), and an unclassified genus of Astroviridae, were found in multiple lizard species (Fig. 1). Five viral families were found across multiple host species: the Flaviviridae, Picornaviridae, Astroviridae, Rhabdoviridae, and Amnoonviridae.

A number of other vertebrate-associated viral families were detected, but were determined as likely endogenous virus elements (EVEs) because longer contigs had disrupted open reading frames. These included members of the Adintoviridae, Retroviridae, and Chuviridae, which were each present in every library, as well as the Hepadnaviridae (present in all libraries except Omar_M and Cmet_M) and the Circoviridae (present only in Omar_M, Garn_M, Cam_M and Csex_M). The likely endogenous viruses were excluded from plots and the abundance and alpha diversity analyses because they are more likely to be associated with host expression than active viral replication, and because our focus was on exogenous viruses. Of note, library CSex_M (comprising Carlia sexdentata) did not contain any viruses likely to infect vertebrates, aside from some likely EVEs, but did contain some viruses that are not associated with vertebrates (Supplementary Fig. S1). Across all libraries, we detected several viral groups that were unlikely to be associated with vertebrate infection, representing the Permutotetraviridae, Tectiviridae, Mimiviridae, Autolykiviridae, Baculoviridae, unclassified Ortervirales, unclassified Riboviria, and tombus-like, solemo-like, narna-like, partiti-like, and toti-like viruses. None of these viruses were abundant (Supplementary Fig. S1) and were considered likely to be viruses of commensal organisms in the lizard livers (i.e. narna-like) or contaminants. Only vertebrate-associated viral families are included in the analyses described below.

Abundance and diversity of lizard viruses

To help determine the factors influencing virome diversity, we compared the abundance and alpha diversity of vertebrate-infecting viruses between libraries from different lizard species and sampling regions within Australia. Specifically, we considered host taxonomy (at the family level) as well as ecological variables, including habitat, and environment (Omar_M was excluded from these analyses as it was the only library from the Diplodactylidae family, and Csex_M was excluded from these analyses as it was the only riparian library and no biologically relevant viruses were found in this library). Given variability in the number of individuals per pool, we first assessed the effect of the number of samples per pool/library to account for any bias introduced by sampling strategy. The number of samples per pool correlated with the richness (P = 0.037), Shannon diversity (P = 0.027), Shannon Effective (P = 0.03), Simpson diversity (P = 0.03), and Simpson Effective diversity (P = 0.034), but not the abundance (P = 0.190) of vertebrate viruses (Supplementary Fig. S2). As such, the number of samples per pool was considered in all statistical models and included as cofactors in final models addressing all diversity measures aside from abundance.

Host taxonomy and habitat were consistently important parameters in the best-fit models for all diversity measures (except Shannon effective and Simpson effective for host taxonomy), indicating that these factors are important modulators of viral diversity in this sample of lizards. Taxonomy had a statistically significant effect on viral richness (P = 0.019) and Shannon index (P = 1.69e-06), but was not significant in modulating abundance or the Simpson index (Supplementary Fig. S3). Habitat significantly modulated richness (P = 2.1e-4), Shannon index (P = 4.312e-11), Shannon effective (P = 0.002), and Simpson Effective (P = 0.011) diversity, but did not affect abundance or Simpson index (Supplementary Fig. S4). While taxonomy and habitat had a clear effect on virus diversity in these animals, a larger sample is required to perform individual comparisons between groups to determine (for example) which specific habitat or host family harbor a greater number of viruses compared to others.

We also explored beta diversity in relation to host traits of the animals sampled here (Supplementary Fig. S5). A similar trend to alpha diversity was observed, in which host taxonomy (P = 0.044) and habitat (P = 0.029) were statistically significant in modulating beta diversity. Overall, we documented a general trend of host taxonomy and habitat as potential modulators of viral diversity in the lizard populations sampled here.

Evolutionary relationships of novel viruses

Flaviviridae

Strikingly, most of the Flaviviridae identified in the lizards studied here fell within the genus Hepacivirus. The Australian lizard hepaciviruses were all novel and fell into three clades in a host-specific manner (Fig. 2). Hepaciviruses from Carlia amax formed a monophyletic group with a viviparous lizard hepacivirus from Zootoca vivipara (Lacertidae) sampled in the UK (Scotland) and a Gecko hepacivirus previously sampled from Oedura in Australia (although the RdRp sequence from the Gecko hepacivirus was only ninety-two amino acid residues, necessarily impacting phylogenetic accuracy). In addition, the long branch lengths suggest substantial unsampled genetic diversity within this clade. While this clade has strong support, its specific phylogenetic placement within the hepaciviruses is poorly resolved such that its closest relatives are unclear (Fig. 2). Gehyra nana had a high diversity of hepaciviruses, with viral sequences clustering into two groups: one distinct lineage that broadly clustered with mammalian viruses and a second that grouped with reptile and bird viruses (Fig. 2). Viruses in the avian/reptile lineage grouped with viruses from other Gekkota species sampled in China, including a Goniurosaurus luii (Gekkota; Eublepharidae, AVM87252.1), Hemidactylus bowringii (Gekkota; Eublepharidae, AVM87554.1), and a Teratoscincus roborowskii (Gekkota; Sphaerodactylidae, AVM87251.1). This is consistent with the Asian origins of Gehyra (Oliver and Hugall 2017). Viruses from the Oedura marmorata sampled in this study formed their own clade within the avian/reptile host lineage, basal to those from China and G. nana.

Interestingly, the Flaviviridae identified in the Cryptoblepharus metallicus library did not cluster with the other Australian lizard Flaviviridae, but rather with Tanama bat virus and Lumpfish flavivirus (Fig. 2), which are part of a broad ‘flavi-like’ group that are phylogenetically distinct from members of the genus Orthoflavivirus. The long branch lengths between these three viruses suggest that there is considerable unsampled viral diversity in this clade, and that this constitutes a new virus that we have tentatively named Bebevirus.

Picornaviridae

Members of the Picornaviridae were found in all Gekkota species analyzed, with the exception of H. planiceps (although a picornavirus contig from O. marmorata was excluded from the phylogeny as there were no RdRp contigs of sufficient length). The picornaviruses isolated from lizards in this study fell into two main groups. Viruses from H. binoei and G. nana grouped together and shared a common ancestor with bat picornavirus (AIF74248.1; 53.5–62 per cent amino acid [aa] identity) and Chameleon picornavirus 1 (DAZ91100; 50–58.7 per cent aa identity) (Fig. 3), a virus found in multiple lizard species including Kinyongia boehmei (Chamaeleonidae), Podarcis muralis (Lacertidae), and Timon pater (Lacertidae). Bat picornavirus is a member of the genus Shanbavirus, and (based on identity in the RdRp) the H. binoei and G. nana viruses, as well as Chameleon picornavirus 1, should also belong to the shanbaviruses (https://ictv.global/report/chapter/picornaviridae/picornaviridae). Three distinct Shanbavirus species were identified in this study, provisionally named Halpertvirus (G. nana host), Schrutevirus (G. nana host), and Beeslyvirus (H. binoei host). Notably, Beeslyvirus was found in H. binoei from both mesic and arid environments and was the only virus found in more than one library.

The second group contained the G. arnhemica picornavirus sequences that were closely related and hence assigned as the same virus species, here provisionally named Kapoorvirus. This virus was most closely related to worm lizard picornavirus (57.6 per cent aa identity in the RdRp) from Blanus cinereus (Blanidae) and more broadly related to members of the genus Tremovirus (Fig. 3), a group associated with reptiles and birds and known to cause encephalomyelitis in avian hosts. The Kapoorviruses shared 34.3–41.6 per cent amino acid identity in the RdRp with the tremoviruses included in the phylogeny, which placed them on the border of the level of divergence required to signify a new genus (64 per cent divergence). As such, we have tentatively placed them within the genus Tremovirus.

Astroviridae

Members of the Astroviridae were detected in two Heteronotia libraries from two different species—H. binoei and H. planiceps. A distinct astrovirus was identified in each library, designated as Hetplan gecko astrovirus (from H. planiceps) and Hetbin gecko astrovirus (from H. binoei), and the two viruses grouped together in the RdRp phylogeny (Fig. 3). These viruses then grouped more broadly with a larger well-supported clade of viruses from lizards (Fig. 3), clustering in a fashion reflective of host family and genus. Specifically, the viruses from Heteronotia grouped closely with viruses from other members of the Gekkonidae, including Hemidactylus bowringii and Gekko similignum, with each also clustering by host genus, all of which fell basal to a clade of viruses from Agamidae (Phrynocephalus erythrurus and Phrynocephalus theobaldi), mirroring the host phylogeny (Fig. 3). Genus demarcation criteria for the Astroviridae is unclear and is most strongly associated with host species. As a consequence, the viruses found here, along with the published lizard astroviruses, would most likely form their own genus. It is also notable that the lizard astrovirus cluster grouped more broadly with viruses found in amphibians and birds.

Rhabdoviridae

Viral contigs related to the Rhabdoviridae were identified in two libraries (Fig. 1). The H. binoei rhabdovirus fell with a group of viruses not classified to any genus, including those isolated from ray-finned fishes (Actinopterygii; from liver, gut, and gills), and a spotted paddle-tail newt (Pachytriton brevipes; from gut) (Fig. 4). This virus is highly divergent (49 per cent aa identity in the RdRp to the closest relative, Wenling dimarhabdovirus 9) and would therefore represent a new species that we have named Heteronotia binoei rhabdovirus. This virus, together with Wenling dimarhabdovirus 9, Fujian dimarhabdovirus, and Wenling dimarhabdovirus 1, should form at least one new genus based on the observed level of sequence diversity. This group of unclassified viruses then fell as a sister-group to a large clade of rhabdoviruses including lyssavirus and many other pathogenic viruses (Fig. 4). Only short Rhabdoviridae contigs could be obtained for the RdRp region from the G. nana library (Gnan_M), the longest being 124 residues. A phylogeny was estimated using a smaller region of the RdRp (230 aa after trimAl trimming) where this contig aligned to confirm the phylogenic groupings observed in the larger alignment (774 aa after trimAl trimming). Regardless of the RdRp alignment used, the G. nana rhabdovirus grouped with the Sripuviruses (Fig. 4). It shared 75 per cent identity in the RdRp protein to its closest relative—Almpiwar virus, isolated from skinks in Northern Queensland Australia. As all RdRp contigs from Gnan_M (including the one in the phylogeny) exhibited less than 90 per cent amino acid identity to the closest relative, this virus likely represents a new species of Sripuvirus, tentatively named Gehyra nana rhabdovirus.

Amnoonviridae

Several contigs related to Lauta virus (MT386081), a member of the Amnoonviridae, were identified in C. amax and H. binoei sampled from the same area of the north-east Top End. However, the contigs obtained were short and hence difficult to classify. The amino acid sequence identity between the longest contigs (270 aa and 238 aa) from the two libraries is only 26.7 per cent, while the identity of the longest contigs from C. amax and H. binoei with the nearest published relative (MT386081, Lauta virus) were 31.3 per cent and 35.8 per cent, respectively. The nearest known relative, found in an Australian gecko (Gehyra lauta) liver in 2013 (Ortiz-Baez et al. 2020), is itself highly divergent from other published viruses, and was only detected using a protein structure prediction approach. All Amnoonviridae contigs detected in this study were from the RdRp segment as these viruses are so distinct from known viruses that it was challenging to identify other genomic regions, particularly as these viruses are likely to be segmented. Phylogenetic analysis of the RdRp of contigs longer than 100 aa from each library suggested that there were four distinct species of lizard virus within Amnoonviridae, two in each library (Fig. 4). Three of these species clustered with Lauta virus, creating a lizard-infecting clade, while the fourth was basal (yet quite distant) to a clade of fish viruses, including Tilapia lake virus (Fig. 4). However, the clustering of these contigs may be unreliable since they are short (160–270 aa) and not the complete length of the RdRp. Additionally, contigs from the same library did not overlap (or only overlapped minimally) in the alignment and therefore their relationship is difficult to determine.

Viral families in which a single species was identified

A single virus species was detected in the lizard meta-transcriptomes from the Caliciviridae, Arenaviridae, Bornaviridae, and Iridoviridae. A member of the Caliciviridae was detected in the H. binoei (Hbin_A) (Fig. 1), but only partial contigs were obtained which aligned to most of the capsid, parts of the RdRp, 2C-like protein and the proteinase. Alignment of a 449 aa region of the RdRp with published sequences grouped this virus with caliciviruses from other squamates; Cyclophiops major (Colubridae), Phrynocephalus theobaldi (Agamidae), and a frog, Sylvirana guentheri (Ranidae) (Fig. 3). The H. binoei calicivirus was highly divergent from its closest relative—Guangdong greater green snake calicivirus (AVM87221)—sharing only 35.4 per cent identity in the aligned region of the RdRp. We therefore suggest that it comprises a new species that we have tentatively named Newselinavirus. Based on alignment of the near complete major capsid protein (527 aa), this virus has >60 per cent amino acid sequence difference to other caliciviruses and therefore would also constitute a different genus (https://ictv.global/report/chapter/caliciviridae/caliciviridae).

A member of the Arenaviridae was detected in H. planiceps (Hplan_A) with the highest abundance of any virus detected in this study (>0.2 per cent of total reads for that library, Fig. 1) with a near complete genome obtained. Both the S and L segments were detected, both of which had nearly complete coding regions. The S segment was at least 1.5 times more abundant than the L segment (expected count of 96,978 and 63,717 for the S and L segments, respectively) or almost 2.6 times more abundant when adjusted for segment length (fragments per kilobase of transcript per million mapped reads of the S segment was 131,036, and the L segment was 51,039). The coverage of the nucleoprotein gene and glycoprotein gene in the S segment was similar, with a mean coverage of 2,825 and 3,011, respectively, but with a dip in coverage between the two ORFs. The mean coverage for the L segment was lower, at 1,111. Based on the RdRp, this virus clusters with the Hartmaniviruses (Fig. 4), a genus known to infect snakes with unknown pathogenicity, although with only 23.6 per cent aa identity to its closest relative in the RdRp region. However, based on the genus demarcation criteria for Arenaviridae (i.e. members of the same genus share >35 per cent nucleotide identity in the L gene; https://ictv.global/report/chapter/arenaviridae/arenaviridae), this virus falls just within the genus Hartmanivirus. Like other members of this genus, the virus detected here, tentatively named Hetplan gecko arenavirus, lacked the gene encoding the zinc-binding matrix protein (Hepojoki et al. 2015).

A near-complete genome of a member of the Bornaviridae was found in the C. munda library (Cmun_M). Phylogenetic analysis of the RdRp showed that it formed a clade with two snake viruses from the genus Orthobornavirus which also contained Loveridges garter snake virus 1 (YP_009055063) and viruses with mammalian and avian hosts (Fig. 4). As it shared only 53.5 per cent nucleotide identity across the entire genome with its nearest relative, the virus detected here would constitute its own species (Kuhn et al. 2015), tentatively named Carlia munda bornavirus, and would likely belong to the genus Orthobornavirus. It is not clear whether the two snake viruses that grouped with Carlia munda bornavirus are pathogenic (Pfaff and Rubbenstroth 2021), but other members of the Orthobornavirus clade cause neurotropic disease in mammals and avian hosts, while members of the Carbovirus genus cause neurological disease in snakes (Hyndman et al. 2018).

Transcripts from multiple genes related to the Iridoviridae (double-strand DNA viruses) were detected in H. binoei (Hbin_M) (Fig. 1), including the major capsid, myristylated membrane protein, phosphotransferase, ATPase, DNA helicase, dUTP, and NTPase I. This suggested that an Iridoviridae was actively replicating in one or more lizards sampled in this library. Phylogenetic analysis of the major capsid protein revealed that this virus clustered within the subfamily Betairidovirinae that was originally considered invertebrate-specific but which is increasingly being associated with vertebrate infection (Russo et al. 2021) (Fig. 5). Interestingly, this virus grouped with erythrocytic necrosis virus, a pathogen of fish (72.4 per cent aa identity in the capsid), and Rhinella marina erythrocytic-like virus found in cane toads (85.4 per cent aa identity) and would likely fall into the same genus (https://ictv.global/report/chapter/iridoviridae/iridoviridae), but as a novel species that we have named Heteronotia binoei irido-like virus (Fig. 5). A new genus containing these and erythrocytic viruses of ectothermic hosts has been putatively proposed, but not yet formally ratified or named (Emmenegger et al. 2014). Three other erythrocytic viruses from Squamata hosts—Thamnophis sauritus erythrocytic virus (EV), Lacerta monticola EV, and Pogona vitticeps EV—fell within this erythrocytic virus lineage based on sequence of a short PCR product of a homolog to a region of the DNA-dependent DNA polymerase (Grosset et al. 2014; Russo et al. 2021), although the transcript encoding this particular protein was not found in our data, precluding a phylogenetic analysis. However, it is likely the virus discovered here is a relative to the other Squamata erythrocytic viruses and supports the presence of Squamata hosts for this lineage.

Discussion

To date, relatively little attention has been paid to the viromes of squamates (lizards and snakes) (Shi et al. 2018; Harding et al. 2022). In Australia, the Squamata are the most species-rich vertebrate assemblage (Wilson and Swan 2017), and the varied ages of Squamate radiations (Brennan and Oliver 2017), alongside their varied population structures and their adaption to multiple habitats and environments (Pianka 1969, 1973), make them an informative group to study in the context of viral evolution and ecology. Using a meta-transcriptomic approach, we characterized the virome of nine lizard species from three families and five genera, finding a wide diversity of novel viruses, including those from viral families/genera that commonly cause disease in humans and other animals, including the Flaviviridae (Hepacivirus), Bornaviridae, Rhabdoviridae, and Picornaviridae. Our findings expand many known viral families, adding entire new clades. This includes a new lizard-specific clade of the Amnoonviridae, a family previously containing only fish viruses and a single lizard-infecting virus (Ortiz-Baez et al. 2020).

The arenavirus detected in Heteronotia planiceps exhibited a remarkably high abundance. The S segment was 2.6 times more abundant than the L segment, relatively consistent with the expected ratio of S segment to L segment RNAs seen in other arenavirus virions and during infection (Meyer, De La Torre, and Southern 2002; Haist, Ziegler, and Botten 2015). The coverage of the nucleoprotein gene and glycoprotein gene on the S segment was similar but with a dip in coverage between the two ORFs, perhaps indicating that this virus was likely in active replication since these two ORFs are ambi-sense and their mRNAs are therefore synthesized separately, with heterogeneous expression of the intergenic region (Meyer and Southern 1993). The significance of the high abundance of the arenavirus in this library is not known, but it is notable that divergent Arenaviridae have also been detected in very high abundance in multiple fish species across several studies, including the pygmy goby arenavirus (Geoghegan et al. 2021), eastern mosquitofish arenavirus (Costa et al. 2021), and Wenling frogfish arenavirus 1 (Shi et al. 2018). Notably, in the cited studies, additional Arenaviridae were detected in lower abundance in the same or other host species, indicating that arenaviruses are not always detected at high abundance, and perhaps the abundance is relevant to the stage of infection or the specific virus.

Viral diversity and abundance differed markedly between the hosts sampled here, and consistent with other studies on virome and pathogen ecology, environmental variables had some effect on virome diversity across our sampled hosts—particularly host habitat. The effect of various habitats on alpha diversity in viruses has previously been demonstrated in fish (Geoghegan et al. 2021), rodents (Tirera et al. 2021), and bats (Bergner et al. 2020), while virome composition has also been shown to differ significantly in different habitats for bats, shrews, and rats (Chen et al. 2023). As different habitats have a range of environmental factors that modulate host behavior and ecology, it is unsurprising that host habitat may also impact virome composition. For example, foraging behavior is likely to play a role in the transmission of low pathogenic avian influenza, where avian species that forage in shallow water are more likely to influence avian influenza ecology (Wille et al. 2023).

Host taxonomy also impacted richness, Shannon diversity, and beta diversity of the viruses sampled here, consistent with some degree of host-specific virus evolution, including co-divergence. The effect of host taxonomy on viral diversity is similar to that seen in fish (Geoghegan et al. 2021) and mammalian viromes (Olival et al. 2017), although a different pattern was observed in some bird viromes (Wille et al. 2018, 2019), suggesting that it is dependent on the taxa in question and that more expansive sampling is required. Over evolutionary time the likelihood of successful cross-species virus transmission is predicted to decline because of host genetic differences in viral receptors, the cellular machinery required for replication, immune response, and other factors that affect viral infection and spread (Ferris, Heise, and Baric 2016). This may in part explain the link between host taxonomy and viral diversity and composition. Alternatively, this association could simply reflect virome-level virus-host co-divergence over longer evolutionary timescales. Indeed, in most cases, the lizard viruses discovered here grouped together and/or with viruses found in other members of the reptilian class Lepidosauria, and phylogenies generally followed a broad pattern of host-virus co-divergence. For example, in the case of the Astroviridae, virus sequences from Squamata clearly group by host genera, and viruses from the Gekkonidae fall basal to those from the Agamidae, matching the host phylogeny (Oliver and Hugall 2017). Similarly, in the case of the relatively well-sampled hepaciviruses and picornaviruses, more closely related Squamata hosts tended to carry more closely related viruses, with viruses from the same host species forming a monophyletic group and viruses from the same family and infraorder similarly tending to do the same.

As well as providing some evidence for virus-host co-divergence, we documented many instances of virus-host jumping or multi-host viruses at higher taxonomic scales. Indeed, in some cases, lizard viruses clustered more closely with amphibian, mammalian, or fish viruses, rather than with other reptile or bird viruses. Three distinct clades of lizard hepaciviruses were identified, all of which group separately to a clade of turtle hepaciviruses. Although the hepacivirus phylogeny is difficult to resolve, the presence of four clades of reptile hepaciviruses and the placement of two reptile clades within the mammalian hepacivirus clade is suggestive of host jumping between reptiles and mammals at some time in the distant evolutionary past. Similarly, the Iridoviridae and Rhabdoviridae phylogenies are highly incongruent with their host phylogenies, suggesting frequent host jumping. In the case of the two lizard Rhabdoviridae, one groups with fish viruses while the other clusters in the Sripuvirus clade with viruses from reptiles, amphibians, and invertebrates. The closest relative of the Sripuvirus detected here, Almpiwar virus, was isolated from skinks in Northern Queensland Australia, and neutralizing antibodies against Almpiwar virus were also found in the sera of crocodiles, a wild bird, and multiple mammalian hosts, including a human (McAllister et al. 2014). Although, the presence of neutralizing antibodies against Almpiwar virus in mammals was relatively rare, it does suggest that this group of viruses could have a broader host range than reptiles alone. In addition, there is evidence that Almpiwar virus, along with other members of this genus—including Charleville virus, which was also isolated in Australia from sandflies and Gehyra (the same host genus as the virus detected here)—are arthropod-borne (McAllister et al. 2014; Vasilakis et al. 2019). The lizard Iridoviridae found here was most closely related to an amphibian and fish virus, with a relatively short genetic distance between them. This is perhaps unsurprising as recent host jumping seems to be commonplace within the Iridoviridae; some members of this family can infect fish, reptiles, and amphibians (Brenes et al. 2014), and some invertebrate-infecting Iridoviridae can also apparently infect vertebrates (Papp and Marschang 2019). In sum, this work supports a model of frequent host-jumping throughout the evolutionary history of most virus families on a backdrop of virus-host co-divergence (Geoghegan, Duchêne, and Holmes 2017; Geoghegan and Holmes 2017; Shi et al. 2018).

Our phylogenetic analyses also revealed some geospatial clustering. For virus families with multiple detections, there is little overlap of virus species from different environments. This suggests that lizard viruses are evolving within their host populations in a relatively isolated manner, as might be expected given the generally strong phylogeographic structuring across species ranges (Moritz et al. 2016, 2018; Potter et al. 2018). It is also of interest that phylogenetic groups that contain viruses sampled from both mesic and arid locations tend to have mesic associated viruses in the basal positions. The hepacivirus phylogeny is particularly enlightening in this respect and it suggests the movement of viruses from one environment to another: both clades of G. nana hepaciviruses have viruses collected from mesic environments located in basal phylogenetic positions to those viruses sampled from animals in arid environments. This is consistent with the theory that Australia’s mesic terrestrial biota is mostly ancestral (Byrne et al. 2011). Additionally, more consistent mesic conditions may favor the retention of viral diversity, such that younger taxa dominate in the more climatically variable arid regions of the Kimberley. This could also be associated with the more variable demographic histories of the drier Kimberley compared to the mesic Top End populations (Potter et al. 2018).

The hepacivirus phylogeny also reflects host historical biogeography. As >90 per cent of Australian reptiles are endemic (Chapman 2009), it is expected that viruses infecting Australian reptiles would have minimal opportunities to spread between countries. Indeed, even the currently limited sampling of lizard viruses demonstrates little movement between countries. However, in one clade of hepaciviruses, viruses from Australian O. marmorata hosts were separated from viruses from Australian G. nana by viruses sampled from China, with the Chinese sampled viruses falling as sister taxa to the Australian G. nana viruses. Interestingly, O. marmorata are thought to have Gondwanan origins, while the G. nana are believed to have immigrated from Asia around the Eocene-Oligocene transition (Oliver and Hugall 2017). Thus, this clade of hepaciviruses could reflect the historical biogeography of their hosts, with evidence of a lineage of viruses introduced to Australia within immigrant ancestors. As the evolutionary histories of the viruses studied here frequently aligned with the biogeographic and paleoclimatic histories of the hosts, this supports the role of virus evolution in informing animal host ecological histories where sampling is strategic and sufficiently dense (Wilfert and Jiggins 2014). Clearly, however, increased sampling across taxa and biogeographic regions would be useful in confirming this hypothesis.

Although the disease association, if any, of the viruses collected here is unknown, we did identify viruses related to known pathogens, including an Orthobornavirus, multiple Hepaciviruses, and a member of the Iridoviridae. It is noteworthy that seemingly healthy lizards can carry a very high abundance and richness of viruses in the liver, particularly the Gehyra species (high viral richness) and H. planiceps (high abundance). Viral emergence can have devastating effects on reptile species, such as the Bellinger River snapping turtle that is now critically endangered following mass mortalities due to a novel nidovirus outbreak (J. Zhang et al. 2018). Although it is unclear where the novel nidovirus (with kidney tropism) originated, the closest known viral relatives were observed in pythons and lizards (usually associated with respiratory disease) (J. Zhang et al. 2018). Given that lizards make up 59 per cent of non-avian reptiles (Pincheira-Donoso et al. 2013), and one-fifth of reptile species are threatened (Cox et al. 2022), this expansion of the lizard virosphere and elucidation of viral genomic sequences may help inform and prepare for potential threats to reptile species.

While this work is the first structured examination of the lizard virome, the sample size is relatively small, restricted to specific Australian taxa, and there are several caveats. Previous studies have found variation in viromes depending on the age of the individuals sampled (Bergner et al. 2020) and the time of year that sampling occurred (Raghwani et al. 2022). These variables were not controlled in this study. In addition, a larger sample size would enable better comparison between matched species sampled from different environments and a greater sample across the squamate phylogeny would expand these results. The viral families detected here are reflective of the fact that liver samples were sequenced, and as viromes commonly differ by the tissue of sampling, studying a wider range of tissues would be beneficial. Despite these limitations, this study demonstrates that lizards carry a large diversity of viruses, often in high abundance and potentially species-specific. As such, they are not only interesting models of vertebrate evolution and ecology, but also serve as good host models by which to study virus ecology and evolution. This work further demonstrates that virus evolution may be a useful tool for understanding or corroborating host biogeographic and paleoclimatic histories.

Data availability

All raw data (fastq files) generated for this study are available in the NCBI SRA database under BioProject PRJNA1078569, BioSample accessions SAMN40006542-52, and SRA accessions SRR28055538-48. Consensus sequences of assembled viral contigs presented in phylogenies are available in GenBank under accession numbers PP711172-PP711215. Scripts for statistical analyses are published at https://github.com/michellewille2/LizardViromes.

Supplementary data

Supplementary data is available at Virus Evolution Journal online.

Acknowledgements

We acknowledge the logistical support provided by Niccy Aitken and Leo Tedeschi at Australian National University while completing the wet-lab component of this study. We acknowledge the Sydney Informatics Hub and the University of Sydney’s high-performance computing cluster Artemis for providing the high-performance computing resources that contributed to the research results reported within this paper.

Funding

This work was supported by an Australian Research Council Australian Laureate Fellowship awarded to E.C.H. (grant number FL170100022).

Conflict of interest:

None declared.

References

Author notes