Meta-Transcriptomic Discovery of a Divergent Circovirus and a Chaphamaparvovirus in Captive Reptiles with Proliferative Respiratory Syndrome.

Viral pathogens are being increasingly described in association with mass morbidity and mortality events in reptiles. However, our knowledge of reptile viruses remains limited. Herein, we describe the meta-transcriptomic investigation of a mass morbidity and mortality event in a colony of central bearded dragons (Pogona vitticeps) in 2014. Severe, extensive proliferation of the respiratory epithelium was consistently found in affected dragons. Similar proliferative lung lesions were identified in bearded dragons from the same colony in 2020 in association with increased intermittent mortality. Total RNA sequencing identified two divergent DNA viruses: a reptile-infecting circovirus, denoted bearded dragon circovirus (BDCV), and the first exogeneous reptilian chaphamaparvovirus-bearded dragon chaphamaparvovirus (BDchPV). Phylogenetic analysis revealed that BDCV was most closely related to bat-associated circoviruses, exhibiting 70% amino acid sequence identity in the Replicase (Rep) protein. In contrast, in the nonstructural (NS) protein, the newly discovered BDchPV showed approximately 31%-35% identity to parvoviruses obtained from tilapia fish and crocodiles in China. Subsequent specific PCR assays revealed BDCV and BDchPV in both diseased and apparently normal captive reptiles, although only BDCV was found in those animals with proliferative pulmonary lesions and respiratory disease. This study expands our understanding of viral diversity in captive reptiles.

1. Introduction

Squamates (scaled reptiles) are one of the most diverse groups of vertebrate fauna in Australia [1]. Bearded dragons (Pogona spp.), including central bearded dragons (Pogona vitticeps), are native reptiles that inhabit wide geographic areas throughout the arid and semiarid regions of Australia. Bearded dragons are increasingly popular exotic pets, and are now one of the most common companion lizard species in Australia, Europe, Asia, and North America [2]. Bearded dragons are also ideal for laboratory studies: their well-established history in captivity makes them suitable for investigations of reptile physiology, and they have unusual patterns of chromosomal, genomic, and temperature-based sex-determination, including sex reversal with specific behavioral consequences [3,4]. Despite their significance, our knowledge of reptile pathogens, particularly viruses, remains limited. Indeed, to date, the only viruses known to infect bearded dragons are double-strand DNA viruses from the families Iridoviridae [5,6,7,8] and Adenoviridae [9,10], single-stranded DNA viruses from the Parvoviridae [11,12,13] and RNA viruses from the Paramyxoviridae [14,15].

The diagnosis of microbial infections in reptiles is complex and often associated with multiple synergistic and predisposing factors. It is likely that undetected viruses contribute to a number of infectious disease outbreaks in captive reptiles, including bearded dragons. To this end, we investigated two unusual mortality events of bearded dragons in a research colony in Australia in 2014 and 2020. In 2014, approximately 40 central bearded dragons died after emergence from brumation in a colony of over 400 animals located in Canberra, Australia. Affected animals were either found dead or were found listless and died within 24 h. In 2020, the same colony experienced low-grade post-brumation mortality followed several months later by star gazing and poor body condition in juveniles and acute death in an adult bearded dragon.

Brumation is an extreme hypometabolic state used by some reptiles to cope with low or unpredictable food availability and unfavorable seasonal conditions for the duration of winter. A number of specific protective measures have been identified during this hibernation in bearded dragons, including increased neuroprotection in the brain, maintenance of heart function through hypertrophy, and upgrading antioxidant capacity and mitochondrial maintenance by skeletal muscle atrophy [16]. Emergence from brumation may also represent a period of increased disease susceptibility in reptiles, as shown by downregulated transcription of several genes responsible for microbial pathogen defense, cellular and oxidative stress, and cell differentiation and growth [17].

Herein, we employed PCR testing and meta-transcriptomic approaches combined with gross and microscopic pathology to investigate the possible involvement of viral pathogens associated with unknown mortality events in a captive colony of bearded dragons.

2. Materials and Methods

2.1. Sample Collection and Processing

Tissues were collected from two disease outbreaks in a research colony in Canberra, Australia. Five dead bearded dragons were submitted for postmortem examination during a mortality event in 2014. In 2020, 12 live, sick animals (including 8 central bearded dragons (Pogona vitticeps), 3 jacky lizards (Amphibolurus muricatus), and 1 grassland earless dragon (Tympanocryptis pinguicolla)) were examined, sedated with 20 mg/kg of im alfaxalone (Alfaxan® CD-RTU 10 mg/mL, Jurox, Australia), delivered into the biceps to rapidly induce a surgical plane of anesthesia, followed by 1 ml/kg of iv pentobarbitone (Lethobarb, 325 mg/mL, Virbac, Australia) to effect euthanasia. Gross postmortem examinations were conducted immediately, and fresh portions of cervical spinal cord, liver, heart, spleen, and kidney were collected aseptically and frozen at −80 °C. A range of tissues were fixed in 10% neutral buffered formalin, processed in ethanol, embedded with paraffin, sectioned, stained with hematoxylin and eosin, and mounted with a cover slip prior to examination by light microscopy. Giemsa, Ziehl–Neelsen, and periodic acid-Schiff stains were applied to a subset of embedded tissue samples to identify bacteria, and to exclude the presence of mycobacteria and fungi. Aerobic, anaerobic, and fungal cultures were undertaken on lung and liver samples aseptically collected from the five dragons that died in 2014.

Proliferation of respiratory epithelium was graded on a scale of 0–4, where 0 indicated no discernable lesion and 4 represented severe and extensive proliferation of infundibular and mesobronchial epithelium. Samples were collected under the Opportunistic Sample Collection Program of the Taronga Conservation Society Australia’s Animal Ethics Committee (Approval #3b/12/18, 20 December 2018) and under scientific license #SL100104 issued by the NSW Office of Environment and Heritage and valid 2 May 2011 to 30 April 2021. Samples were submitted to Murdoch University for PCR testing of pneumotropic viruses, and to the University of Sydney for meta-transcriptomic analysis and further PCR testing.

For PCR testing for pneumotropic viruses, total nucleic acid was extracted from tissues using the MELTTM Total Nucleic Acid Isolation System (Ambion, TX, USA) according to the manufacturer’s instructions. Total nucleic acid was eluted into 30 µL elution buffer. For metatranscriptomics and PCR testing for parvoviruses and circoviruses, the RNA extracted from the spinal cord, liver, lung, and kidney of diseased animals was processed using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany). Initially, frozen tissue was partially thawed and submerged in RLT lysis buffer containing 1% β-mercaptoethanol and 0.5% reagent DX before being homogenized using a hand-operated TissueRupture (Qiagen). To pool samples in equal proportions, RNA concentrations and integrity were validated using a NanoDrop spectrophotometer (ThermoFisher Scientific, MA, USA) and a TapeStation (Agilent, CA, USA). Illumina TruSeq stranded RNA libraries were prepared on the pooled samples following rRNA depletion using a RiboZero Gold rRNA removal kit (Epidemiology). Finally, 150 bp paired-end sequencing of the rRNA-depleted RNA library was generated on an Illumina NovaSeq 6000 service system at Australian Genome Research Facility (AGRF), Melbourne.

2.2. PCR Testing for Pneumotropic Viruses

Nucleic acid from the kidney, liver, and lung was pooled and tested for sunshineviruses, orthoreoviruses, and ferlaviruses using one-step reverse transcription (RT)-PCR. For each virus genus, 1 µL extracted nucleic acid was added to 0.8 µL SuperScript® III RT/Platinum® Taq Mix (Invitrogen, Victoria, Australia), 10 µL 2× Reaction Mix, and 1 µM (final concentration) of each primer, and was then made up to a final volume of 20 µL. For nested PCRs, 1 µL PCR product was used as template for the second round of amplification using Platinum® PCR Supermix (Invitrogen) in a final volume of 20 µL. For sunshineviruses, the primer pair SunshineS2–SunshineAS2 was used [18]. For orthoreoviruses, the primer pairs 1607F–2608R and 2090F–2334R were used for the first and second rounds, respectively [19]. For ferlaviruses, the primer pairs L5–L6 and L7–L8 were used, respectively [20]. Sunshine Coast virus [21], a reptile orthoreovirus (kindly provided by Dr Rachel Marschang) and a ferlavirus (American Type Culture Collection VR-1408) were used as positive controls. PCR products were visualised with agarose gel electrophoresis.

2.3. Pathogen Discovery Using Meta-transcriptomics

Our meta-transcriptomic approach to pathogen discovery was based on those previously employed by our group [22,23]. RNA sequencing reads were trimmed for quality using Trimmomatic [24] before de novo assembly with Trinity, version 2.5.1 [25]. Assembled sequence contigs were annotated using both nucleotide BLAST searches against the NCBI nt database and Diamond Blastx against the NCBI nr database [26], with e-value cut-offs of <10−6 and <10−5, respectively. Open reading frames were predicted from the potential viral contigs using Geneious v11.1.5 [27], with gene annotation and functional predictions made against the Conserved Domain Database (CDD) [28]. To evaluate the virus abundance and coverage, reads were mapped back to the genome with BBmap [29].

2.4. PCR Testing for Chaphamaparvoviruses and Circoviruses

Virus-specific PCR primers were designed based on identified transcripts from the RNA-seq data (Table 1). Accordingly, the extracted RNA (5 µL) from each tissue of individual cases was reverse transcribed using the SuperScript IV VILO cDNA synthesis system (Invitrogen). The cDNA generated from each sampled tissue (2 µL) was used for viral specific PCRs targeting regions identified by RNA-Seq. For BDCV, the primer pairs BDCV_F2 and BDCV_R1 were used (cycling conditions = 98 °C × 1 m, 35 × (98 °C × 10 s, 62 °C × 10 s, 72 °C × 45 s)). For BDchPV, the primers BDChPV_F4 and BDChPV_R3 were used (cycling conditions 98 °C × 1 m, 35 × (98 °C × 10s, 62.2 °C × 10 s, 72 °C × 60 s)). All PCRs were performed using Platinum SuperFi DNA polymerase (Invitrogen) with a final concentration of 0.2 μM for both forward and reverse primers. PCR products were visualized with agarose gel electrophoresis and Sanger sequenced at the AGRF.

Table 1

PCR primers used in this study.

Primer Name	Target	Sequence 5′-3′	Site Start	Site End	Tm (°C)	GC(%)
BDCV_F1	BDCV_cap	GTATTGTTTCAGAGATTCGCCG	65	86	57.8	45.5
BDCV_R1		ATTGTAATGAAGGCGTTGTTCC	1144	1165	57.6	40.9
BDCV_F2	BDCV_cap	GCCGAAGTGTCTGGTAATACTA	974	995	57.3	43.5
BDCV_R2		GATTGGTTAAGAGTCGGGATTTC	440	462	57.6	40.9
BDCV_F3	BDCV_rep	AATACCTGACCAGTCTGCCA	1551	1570	57.3	50.0
BDCV_R3		CGATGGCGTAGGAACAGATG	659	640	59.3	55.0
BDCV_F4	BDCV_rep	GGTGAAGTCGATTGTTCCTTTG	312	333	58.4	45.5
BDCV_R4		GTGGGTATTGTTTCAGAGATTCG	1633	1611	58.9	43.5
BDchPV_F1	BDchPV_NS	CTGAACAATATCTCTTCGTTGGAG	217	240	57.4	41.7
BDchPV_R1		TCCCAGACACCCATATAACAATC	1490	1512	57.7	43.5
BDchPV_F2	BDchPV_NS	GCTCTCCTGAATATTGTGATATGTC	954	978	57.4	40.0
BDchPV_R2		TTCCCTACCCACTTCTTTCTG	2255	2275	56.9	47.6
BDchPV_F3	BDchPV_VP	CTAGCTCAGACATCAGATCCAG	2004	2025	57.5	50.0
BDchPV_R3		GTCGGTGGAGGATTGTCTATAG	3346	3367	57.5	50.0
BDchPV_F4	BDchPV_VP	AGCTTTTCCCAAATGACACTG	2848	2868	57.0	42.9
BDchPV_R4		CAGAAGACACAACAGTAGGTAAAG	4087	4110	57.2	41.7

BDCV: Bearded Dragon CircoVirus, BDchPV: Bearded Dragon ChaphamaParvoVirus, rep: replicase, cap: capsid, NS: Non-Structural protein, VP: Viral capsid Protein.

2.5. Phylogenetic Analysis

Phylogenetic trees were estimated based on MAFFT [30] alignments of the conserved nonstructural region of the viruses identified here, with all analyses utilizing representative members of each virus family taken from NCBI/GenBank. Maximum likelihood (ML) phylogenies were inferred using IQ-TREE (version 1.4.3) [31], employing the best-fit LG+F+Γ4 model of amino acid substitution on the nonstructural (NS) protein (717 amino acids) and the LG+F+I+Γ4 model on the replicase (rep) protein (274 amino acids). Statistical support for individual nodes was estimated via bootstrap resampling (1000 replicates). Data were visualized using Figtree 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/). The genome sequence of BDCV and BDchPV are available on NCBI/GenBank (accession numbers: MT732118 and MT73219). Raw sequencing reads are available at the Sequence Read Archive (SRA) under accession PRJNA644669.

2.6. Mining the Sequence Read Archive (SRA)

To help determine if BDCV and BDchPV were present in other bearded dragon or other reptile species, we screened publicly accessible high-throughput sequencing data available on the NCBI SRA database (https://www.ncbi.nlm.nih.gov/sra). Accordingly, a large RNA-seq data set was obtained using the NCBI SRA toolkit (version 2.9.2) from the bearded dragon genus Pogona (NCBI taxid:52201). Retrieved FASTQ reads were then subjected to a blast analysis using Diamond v.0.9.25 [32] against customized databases containing the core genes from reference circoviruses and parvoviruses, employing an e-value of cut-off of 1 × e−5.

3. Results

3.1. Clinical and Histopathological Findings of Diseased Captive Reptiles

The mass mortality and morbidity of up to 40 adult central bearded dragons in 2014, and intermittent ill health and star gazing in hatchlings from the same research colony in 2020, were both accompanied by unusual and severe proliferation of pulmonary epithelium. In total, 17 cases were investigated in this study. The severe proliferation of pulmonary epithelium was found in 7 cases, including 5 bearded dragons from the 2014 outbreak, and 1 bearded dragon and 1 Jacky lizard from 2020 mortality events. The infundibular mucosa, which is normally squamous, exhibited metaplasia ranging from columnar to pseudostratified columnar and ciliate. The mesobronchial epithelium was similarly thickened, pseudostratified, and dysplastic with multifocal loss of cilia. Amorphous basophilic inclusions were evident within enlarged nuclei of infundibular epithelial cells in two out of seven affected animals. Eosinophilic intracytoplasmic inclusion bodies were evident within the respiratory epithelium of three out of seven affected animals. Rarely binucleate epithelial cells and syncytia were identified within mesobronchial epithelium.

Similar inclusions were evident within hepatocytes of three out of five dragons and renal tubular epithelium of two out of five dragons from the initial 2014 mortality event. Pasteurellaceae bacteria of unknown species were isolated from the lung and liver of each of the five dragons examined in 2014, and Aeromonas hydrophila was isolated in the lung of two of these animals; however, these were not considered contributory to pulmonary proliferation nor inclusion body formation. A summary of the signalment, gross, and histological findings in each animal in the initial 2014 outbreak is summarized in Table S1, and the histological lesions are shown in Figure 1.

Figure 1

Photomicrographs of pulmonary lesions in central bearded dragons from the 2014 and 2020 outbreak. Thickening of the normally squamous infundibular epithelium. (Haematoxylin and eosin (HE) stain 40×) (A) and the segmental nature of this change observed in an affected animal from the 2020 outbreak (HE 20×) (B) Metaplasia of the normally squamous infundibular epithelium to pseudostratified columnar and multifocally ciliate. Normal lung appears on the right and infundibular epithelium becomes proliferative on the left (HE 1000×) (C) and multifocal epithelial cells with large, irregularly shaped, basophilic nuclei (arrowheads, HE 1000×) (D) Bizarre atrial epithelial cells containing eosinophilic cytoplasmic inclusions (arrowheads, (HE 1000×) (E).

3.2. PCR Testing for Pneumotropic Viruses

All samples were PCR-negative for sunshineviruses, ferlaviruses, and orthoreoviruses.

3.3. Metatranscriptomic Pathogen Discovery

Sequencing of a pooled rRNA-depleted RNA-seq library resulted in 81,188,754 raw reads. After filtering, 54,804,486 paired-end trimmed reads were generated, which were then de novo assembled into 266,118 contigs. Analyses of these read data revealed the presence of two DNA viruses, denoted here as bearded dragon circovirus (BDCV) and bearded dragon chaphamaparvovirus (BDchPV), with read abundances of 0.0006% and 0.003%, respectively. Virus-targeted PCR assays (Table 1) and Sanger sequencing were used to recover the full genome sequences of these two viruses. Virus-specific PCR was also used on spinal cord, lung, liver, and kidney to provide insight into the distribution of these viruses. No other reptilian associated circoviruses and parvoviruses were detected using SRA screening.

3.4. Genome Characterization of a Novel Bearded Dragon Circovirus

We identified an abundant circoviruslike contig (1693 bp, 304 transcripts per million (TPM)) from our pooled lung, kidney, and liver samples. The typical features of the genus Circovirus are nonenveloped, icosahedral, single-stranded circular DNA (ssDNA) genomes of approximately 1.5–2 kb in size. The full genome of BDCV was recovered through PCR assays and Sanger sequencing. The BDCV genome comprises 1761 bp of circular DNA with a 43.8% GC content, encoding rep and cap genes of 308 and 218 amino acid residues, respectively.

The replicase protein of circoviruses introduces an endonucleolytic nick within the typical stem-loop structure with nonamer sequence 5′-AGTATTAC-3′ in the intergenic region of the genome, thereby initiating the rolling circle replication (RCR). Comparison with other members of the Circoviridae revealed the presence of the classical conserved motifs, including RCR motifs, as well as SF3 helicase motifs (Table 2 and Figure 2).

Table 2

Rep protein (amino acid) motifs detected in the bearded dragon circovirus.

Name of Region		Sequence in the Circoviridae	Sequence in Bearded Dragon Circovirus
RCR motif	RCR motif I	FT(L/I)NN	FTINN
	RCR motif II	PHLQG	PHIQG
	RCR motif III	YC(S/x)K	YCSK
SF3 helicase motif	Walker A	G(P/x)(P/x) GxGK(S/t)	GEPGSGKS
	Walker B	uuDDF	ILDDF
	Walker C	uTSN	ITSN

Figure 2

Genome characterization and phylogeny of bearded dragon circovirus. (A) Genome organization of bearded dragon circovirus. The outermost circles of the metadata ring represent the GC content (%, green), nucleotide polymorphism (%, orange), and read mapping coverage (yellow) of the genome. In the innermost circles, the proteins encoded by the replicase (Rep) and capsid protein (Cap) are labelled in pink and cyan, respectively. (B) Phylogenetic analysis of the Rep gene (nonstructural replicase protein) of circoviruses, including members of the Cyclovirus and Circovirus genera. Bootstrap values > 70% are presented for key nodes (1000 replicates). The tree was midpoint rooted for clarity only. Scale bar shows the number of substitutions per site. (C) Comparison of nonamer sequences with closely related viruses in the Circoviridae.

To determine the evolutionary history of BDCV with respect to other circoviruses, we performed a phylogenetic analysis of the Rep protein. Notably, BDCV clustered robustly (99% bootstrap support; 77.85% amino acid pairwise identity in Rep genes) with a bat-associated virus lineage that includes Tadarida brasiliensis circovirus 1 (accession number YP_009170674.1) and Rhinolophus ferrumequinum circovirus 1 (YP_009506275.1), recently identified in metagenomic analyses of bat gut samples from Brazil and China, respectively. At the nucleotide level, the BDCV genome showed 69.33% and 71.53% sequence similarity with Tadarida brasiliensis circovirus 1 and Rhinolophus ferrumequinum circovirus 1, respectively. Since the species demarcation threshold for circoviruses is <80% genomewide pairwise identity, we suggest that BDCV likely represents a new species within the Circoviridae.

3.5. Identification of a Highly Divergent Bearded Dragon Chaphamaparvovirus

We identified two abundant chaphamaparvoviruslike transcripts in the RNA-seq library. The Parvoviridae are a family of small, nonenveloped, ssDNA animal viruses with a linear genome of 4–6 kb in length. Chaphamaparvoviruses are a recently defined genus in which a variety of component virus species have now been identified [33,34,35,36]. Four sets of bridging PCR assays were designed to recover a near complete virus genome. This revealed that the near complete virus genome comprised 4181 nt with two distinct ORFs that encode the nonstructural protein (NS, 633 aa) and the structural protein (VP, 600 aa). We further utilized specific primers to amplify the targeted NS (BDchPV_F2 and R1) and VP region (BDchPV_F4 and R3) for virus screening in different organs (Table 1).

Two conserved motifs were identified in the NS protein, including the putative endonuclease metal coordination motif “HLH” and the helicase motif “GPASTGKS” at amino acid positions 110–112 and 310–317, respectively (Figure 3). Additionally, the potential accessory ORF1 (positions 74–358, 95 aa) and ORF3 (p15-like) (positions 136–543, 136 aa), previously identified in most other amniote-associated chaphamaparvoviruses, were both detected in BDchPV. The conserved PLA2 and G-rich motifs widely present in other members of Parvoviridae were absent from the VP protein, as is the case for other chaphamaparvoviruses [37,38].

Figure 3

Genome of bearded dragon chaphamaparvovirus in comparison to related viruses. The amino acid sequence size of the ORFs of each virus is shown. Yellow and orange boxes refer to the NS and VP genes, respectively. The classic motif features HXH (blue asterisk) and GPXNTGKX (green asterisk) on NS genes are labelled.

Phylogenetic analysis based on the predicted amino acid sequence of the complete NS protein revealed that BDchPV fell within the Chaphamaparvovirus lineage, clustering with viruses associated with fish and crocodiles. Tilapia parvovirus is a recently identified Chaphamaparvovirus isolated from the feces of farmed tilapia and crocodiles in China [36]. However, the amino acid sequences of the novel BDchPV identified here shared only 32.7% and 31.6% pairwise identity in NS gene and VP gene with this virus (accession no: MN162688.1), respectively (Figure 4).

Figure 4

Genome characterization and phylogenetic relationships of bearded dragon chaphamaparvovirus. (A) Genome depiction and reads mapping coverage of bearded dragon chaphamaparvovirus. (B) Phylogenetic analysis of the nonstructural protein gene (NS gene) of parvoviruses, including other chaphamaparvoviruses and members of the Parvovirinae, Densovirinae, and Hamaparvovirine. Bootstrap values >70% were presented for key nodes (1000 replicates). The tree was midpoint rooted for clarity only. Scale bar shows the number of substitutions per site.

3.6. Prevalence of BDCV and BDchPV through PCR Screening

The nucleic acid extracted from archived bearded dragon cases in the 2014 (n = 5) and captive reptile cases in 2020 (n = 12) outbreaks were screened using BDCV and BDchPV specific primers, respectively (Table 3). The circovirus was detected from kidney and lung samples from one bearded dragon (case: 10043.3) in the original respiratory disease outbreak in 2014, and of three out of twelve liver samples from the second outbreak (2020) were BDCV positive. In total, BDchPV was detected in for liver samples, one brain sample, and one lung sample, including two original respiratory cases in 2014 and three additional cases in 2020 (Table 3). Collectively, the case prevalence of both viruses from the two outbreaks was 29.4% (five out of seventeen) for BDCV and 29.4% (five out of seventeen) for BDchPV. Sanger sequencing results from PCR products showed 95%–99% nucleotide identity to the index cases.

Table 3

PCR testing of captive reptiles from the 2014 and 2020 outbreaks.

Year	Case no.	Species	Clinical Sign	Organ	PCR Results
					BDCV	BDchPV
2014	10043.1	CBD	Respiratory	Kidney	-	-
				Lung	-	-
	10043.2	CBD	Respiratory	Liver	-	-
	10043.3	CBD	Respiratory	Kidney	+	-
				Lung	+	+
	10043.4	CBD	Respiratory	Liver	-	+
	10043.5	CBD	Respiratory	Liver	+	-
2020	13270.1	CBD	Mild neurological	Brain	-	-
				Liver	-	+
	13270.2	CBD	Normal	Brain	-	-
				Liver	-	-
	13270.3	CBD	Mild nonspecific illness	Brain	-	-
				Liver	+	+
	13270.4	CBD	Neurological	Brain	-	+
				Liver	+	+
	13270.5	CBD	Normal	Brain	-	-
				Liver	-	-
	13270.6	CBD	Neurological	Brain	n.a.	n.a.
				Liver	-	-
	13270.7	CBD	Neurological	Brain	n.a.	n.a.
				Liver	-	-
	13270.8	CBD	Neurological	Brain	n.a.	n.a.
				Liver	-	-
	13270.9	JL	Normal	Brain	-	-
				Liver	+	-
	13270.10	JL	Neurological	Brain	-	-
				Liver	-	-
	13270.11	JL	Neurological	Brain	-	-
				Liver	-	-
	13270.12	GED	Moribund	Brain	-	-
				Liver	-	-
Total: BDCV positive cases: 5/17; BDchPV positive cases: 5/17

CBD: central bearded dragon (Pogona vitticeps); JL: jacky lizard (Amphibolurus muricatus); GED: grassland earless dragon (Tympanocryptis pinguicolla); n.a.: nonapplicable; green: negative; red: positive.

4. Discussion

Our knowledge of the viruses of herpetofauna is expanding rapidly. Herein, we describe the meta-transcriptomic discovery and characterization of two complete genomes of DNA viruses, as well as their prevalence from mass mortality events of captive reptiles.

Circoviruses are small DNA viruses that infect a wide range of vertebrates, with most pathogenic consequences reported in association with immunocompromised or immunosuppressed hosts, including canine [39], psittacine [40], porcine [41], and mink [42] species. In 2019, a circovirus was detected in the liver and gut of black-headed pythons (Aspidites melanocephalus) with spinal osteopathy [43]. The virus was then identified in tissues of a Boelen’s python (Morelia boeleni) and two annulated tree boas (Corallus annulatus). That investigation claimed to be the first to describe a reptilian circovirus, although other circo-like viruses had previously been described in testudines and squamates [44,45]. In all cases, the association between the circovirus and disease was either weak or absent. Endogenous circoviruses have also been detected in snake genomes [46]. It is notable that the bearded dragon circovirus newly identified here clustered with a bat-associated circoviruses lineage. However, these bat viruses were identified through metagenomic sequencing of environmental related samples (i.e., fecal viromes) such that their true hosts remain unclear.

Parvoviruses are associated with disease in a variety of host species, ranging from canines [47] to livestock [48,49,50,51], rodents [36], and humans [52]. Reptile parvoviruses were identified in several snake species (snake adeno-associated virus) [11,53,54,55,56,57] and some lizard species including bearded dragons [11,12]. To date, most of the exogenous reptilian parvoviruses identified belong to the genus Dependoparvovirus, and it has been thought that these viruses require helper viruses (usually adenoviruses) for replication [13], although more recent work is challenging this notion [58]. Chaphamaparvoviruses (ChPVs) are a newly identified genus of parvoviruses [59] identified in rodents, birds, pigs, bats, Tasmanian devils, dogs, cats, primates, and even invertebrates [33,34,35,36,37,60]. Although a ChPV was previously identified in crocodile faeces [36], it is possible that this originated in the tilapia fish fed to these crocodiles. Our documentation of a novel and exogenous reptile chaphamaparvovirus in central bearded dragons extends the host range of these parvoviruses.

Interestingly, coinfection with circoviruses and parvoviruses is increasingly reported in different animal species, often in association with immunosuppression and exacerbated disease severity. For example, previous studies have indicated that Porcine circovirus 2 (PCV-2) associated diseases are augmented by concurrent viral infections such as Porcine parvovirus, and these viruses may serve as important cofactors in the pathogenesis of Porcine multisystemic wasting syndrome (PMWS) [61]. Similarly, Canine circovirus 1 (CaCV-1) and Canine parvovirus 2 (CPV-2) infection are linked to recurrent outbreaks of bloody diarrhea and sudden death in puppies [62], and a novel lorikeet chaphamaparvovirus coinfected with psittacine circovirus (beak and feather disease virus) in wild birds has been described [63]. Recently, a similar coexistence of Tasmanian devil-associated circovirus and a chaphamaparvovirus was identified in a Tasmania devil metagenomic virome study, but with no disease association [64].

Despite our description of these pathogens, their roles in the disease syndromes described here remain uncertain and require additional investigation. Indeed, the examination of feeding and husbandry practices surrounding brumation is warranted, particularly since the animals that did not undergo brumation appeared to be unaffected during the mortality event. Although the association between the two viruses and proliferative lung disease was not clearly established, the role of these viruses as predisposing factors for disease or reduced fitness merits further investigation.

Reptiles comprise a considerable proportion of the vertebrate biomass of many ecosystems and play an important role in the food-web, and the delivery of ecosystem services through pollination, seed dispersal, and pest control. Australia has the highest diversity of lizards of any country, yet our understanding of the biological threats to these animals is scarce. It has been difficult to obtain an accurate understanding of the health status of wildlife that are free-living, maintained for research, education, or species recovery, using traditional diagnostic techniques. This study illustrates the value of metagenomic approaches to disease investigation as a complement to traditional histology and pathogen culture. The identification of bearded dragon circovirus and the bearded dragon chaphamaparvovirus provides insight into reptilian viral diversity and reptile health, and undoubtedly merits broader surveillance.