Virome analysis of New Zealand’s bats reveals cross-species viral transmission among the Coronaviridae

Abstract

The lesser short-tailed bat (Mystacina tuberculata) and the long-tailed bat (Chalinolobus tuberculatus) are Aotearoa New Zealand’s only native extant terrestrial mammals and are believed to have migrated from Australia. Long-tailed bats arrived in New Zealand an estimated two million years ago and are closely related to other Australian bat species. Lesser short-tailed bats, in contrast, are the only extant species within the Mystacinidae and are estimated to have been living in isolation in New Zealand for the past 16–18 million years. Throughout this period of isolation, lesser short-tailed bats have become one of the most terrestrial bats in the world. Through a metatranscriptomic analysis of guano samples from eight locations across New Zealand, we aimed to characterise the viromes of New Zealand’s bats and determine whether viruses have jumped between these species over the past two million years. High viral richness was observed among long-tailed bats with viruses spanning seven different viral families. In contrast, no bat-specific viruses were identified in lesser short-tailed bats. Both bat species harboured an abundance of likely dietary- and environment-associated viruses. We also identified alphacoronaviruses in long-tailed bat guano that had previously been identified in lesser short-tailed bats, suggesting that these viruses had jumped the species barrier after long-tailed bats migrated to New Zealand. Of note, an alphacoronavirus species discovered here possessed a complete genome of only 22,416 nucleotides with entire deletions or truncations of several non-structural proteins, thereby representing what may be the shortest genome within the Coronaviridae identified to date. Overall, this study has revealed a diverse range of novel viruses harboured by New Zealand’s only native terrestrial mammals, in turn expanding our understanding of bat viral dynamics and evolution globally.

1. Introduction

Bats harbour a wide diversity of viruses, some of which have the potential to spillover to other hosts, including humans (Li et al. 2005; Wang and Anderson 2019). Viruses such as lyssaviruses, coronaviruses, Ebola, Hendra, and Nipah are all carried by bats, with significant public health and economic consequences (Mackenzie, Field, and Guyatt 2003; Leroy et al. 2005; Wang and Bats 2007; Ithete et al. 2013; Schatz et al. 2013; Plowright et al. 2015; Cui, Li, and Shi 2019; Epstein et al. 2020). Accordingly, bats have become the focus of many virological surveys to help identify novel viruses and assess the risk they pose for zoonoses (Tan et al. 2021). Despite this, investigations into the viruses present in the native bat species of New Zealand have been limited.

New Zealand is home to only two native extant terrestrial mammals; the lesser short-tailed bat (Mystacina tuberculata) and the long-tailed bat (Chalinolobus tuberculatus) (O’Donnell et al. 2021). A third species, the greater short-tailed bat, M. robusta, has not been sighted since 1967 and is likely extinct (O’Donnell et al. 2021). The lesser short-tailed bat can be further classified into three subspecies: the northern lesser short-tailed bat (M. tuberculata aupourica), the central lesser short-tailed bat (M. tuberculata rhyacobia), and the southern lesser short-tailed bat (M. tuberculata tuberculata). An additional two bat species have been identified in the New Zealand fossil record: Vulcanops jennyworthyae and M. miocenalis (Hand et al. 2015, 2018). Viral sequences from the Papillomaviridae, Polyomaviridae, Caliciviridae, Hepeviridae, Poxviridae, Parvoviridae, Adenoviridae, Picornaviridae, and Coronaviridae families have previously been uncovered in lesser short-tailed bats located on a remote, predator-free, offshore island (Hall et al. 2014; Wang et al. 2015), while viral transcripts from the Coronaviridae have previously been uncovered in mainland long-tailed bats and lesser short-tailed bats using a pan-coronavirus polymerase chain reaction (PCR) approach (Tortosa et al. 2023). Building on this previous work, we aimed to use total RNA sequencing to document complete viromes of these species and, from this, determine whether viruses have exhibited cross-species transmission between New Zealand’s bat species.

Bats are believed to be an important reservoir of a multitude of viruses in part because they roost in large numbers, exhibit species co-habitation, and are able to fly, readily transferring viruses among geographic regions. In addition, bats are thought to be able to tolerate a high abundance of viruses due to unique components of their immune system (Subudhi, Rapin, and Misra 2019; Letko et al. 2020; Van Brussel and Holmes 2022). As of 2021, viruses spanning a total of twenty-four RNA viral families and eleven DNA viral families have been uncovered in hosts belonging to the order Chiroptera (Van Brussel and Holmes 2022). Of the bat-associated viral sequences published in GenBank, 85% were RNA viruses, of which 30% and 24% were from the families Coronaviridae and Rhabdoviridae, respectively (Van Brussel and Holmes 2022).

Around 84 million years ago, the continental crust from which New Zealand was formed separated from Gondwana (Mortimer et al. 2019). All terrestrial mammals present in New Zealand originated from overseas with their arrival aided by winds, ocean currents, and human movement (O’Donnell et al. 2021). The long-tailed bat is believed to have arrived from Australia during the Pleistocene around two million years ago (O’Donnell et al. 2021). It is a member of the Vespertilionidae family, diverging from the last common ancestor approximately 17 million years ago, such that it could have migrated to New Zealand earlier than previously thought (Dool et al. 2016). The Vespertilionidae, which comprises of more than 360 species, is considered one of the most widely dispersed mammalian families in the world (O’Donnell et al. 2021). In comparison, the lesser short-tailed bat is the sole surviving member of the Mystacinidae and based on fossil records migrated from Australia during the early Miocene period approximately 16–18 million years ago (Worthy and Holdaway 1994; Hand et al. 1998). Around this time, the remainder of the mystacinid lineage became extinct in Australia (O’Donnell et al. 2021). Mitochondrial phylogenetic analysis indicates that the lesser short-tailed bat last shared a common ancestor with the Noctilionoidae superfamily around 35 million years ago (Den Bussche RA and Hoofer 2000), again suggesting that this species may have arrived earlier than fossil records propose.

The viral diversity of two species of bats from Australia, Gould’s wattled bat (C. gouldii), and chocolate wattled bat (C. morio), which are closely related to New Zealand’s long-tailed bats, has previously been analysed. Viruses from the Coronaviridae, Adenoviridae, and Paramyxoviridae families have previously been identified in both these species (Prada et al. 2019). Of note, the coronaviruses found in Gould’s wattled bats were closely related to the coronaviruses previously detected in New Zealand’s long-tailed bats, suggesting that virus–host codivergence over extended time scales has likely played a role in their evolution (Prada et al. 2019; Tortosa et al. 2023).

Due to its isolation in New Zealand, the lesser short-tailed bat has developed several atypical behaviours distinct from any other known bat species. In particular, the lesser short-tailed bat is considered the most terrestrial bat in the world, such that it has adapted to spending large periods of time on the forest floor foraging for food (Daniel 1976; O’Donnell et al. 2021). Additionally, the diet of lesser short-tailed bats is considered one of the broadest of any bat species, including nectar, fruit, pollen, terrestrial, and flying invertebrates, as well as potential fungi. This diet allows them to remain active for most of the years unlike most other non-tropical insectivore bat species (Daniel 1976; McNab and O’Donnell 2018; O’Donnell et al. 2021). In contrast, long-tailed bats are strictly insectivores and spend more time flying in comparison to lesser short-tailed bats (McNab and O’Donnell 2018; O’Donnell et al. 2021). During winter, long-tailed bats enter a state of short-term torpor as the number of flying insects, their main food source, reduces (Sedgeley 2003; O’Donnell et al. 2021).

Herein, we sampled guano from both long-tailed and lesser short-tailed bats. We investigated the viromes of these species to determine if viruses have jumped between these hosts and test whether the marked differences in their ecology and behaviour shaped the composition of their virome.

2. Materials and methods

2.1. Bat guano sample collection

Guano from lesser short-tailed bats and long-tailed bats was collected as part of several annual monitoring programmes carried out by the New Zealand Department of Conservation. Guano was collected using tarps, harp traps, and catch bags from eight locations/roosts across New Zealand (Fig. 1). A total of 219 individual guano samples were collected during two sampling periods over summer: February 2020 and January–February 2021, around which time female bats are lactating and feeding non-volant young (January) and the young of the year begin flying (January–February) (O’Donnell 2002a; Sedgeley 2006). New Zealand bats that have reached sexual maturity (one year for males and two–three years for females) mate during early autumn (March), while females give birth to a single pup during the summer (November–December) (O’Donnell 2002a; Sedgeley 2006). More information regarding sample locations, species, as well as the number of individual guano samples within each group, is provided in Supplementary Table 1. Guano samples were submerged into 1 ml of RNAlater and briefly stored at 4°C until they were sent to the University of Otago, Dunedin, where they were stored at −70°C or −80°C until total RNA was extracted.

2.2. Bat guano total RNA extraction

Frozen bat guano samples stored in RNAlater were first defrosted. Approximately 25–50 mg of the defrosted guano was placed in 15 ml RNase-free round bottom tubes. In a fume hood, 500μl of cold TRIzol^TM was added and samples were homogenised for one minute using a TissueRuptor (Qigaen). The homogenate was then centrifuged at 4°C at 12,000 ×g for five minutes. Total RNA from the clear supernatant was then extracted following the TRIzol manufacturers protocol with minor alterations. Briefly, a secondary chloroform phase separation step was added to remove any residual phenol contamination while an additional ethanol wash step was added to remove any residual guanidine contamination. Extracted RNA was quantified using a Nanodrop. RNA from 217 of the 219 guano samples was at suitable concentrations for downstream processing. Equal volumes of RNA from 5-41 individuals were pooled into 12 groups based on bat species and sample location to a total pooling volume ranging from 38.1 μl to 64.5 μl (see Supplementary Table 1).

2.3. Bat guano RNA sequencing

Extracted RNA was subject to total RNA sequencing. Libraries were prepared using the Illumina Stranded Total RNA Prep with Ribo-Zero Plus (Illumina) and 16 cycles of PCR. Paired-end 150 bp sequencing of the RNA libraries was performed on the Illumina NovaSeq 6000 platform using a single S4 lane.

2.4. Virome assembly and virus identification

Paired reads were trimmed and assembled de novo using Trinity v2.11 with the ‘trimmomatic’ flag option (Haas et al. 2013). Sequence similarity searches against the NCBI nucleotide (nt) database (2021) and the non-redundant (nr) protein database (2021) using BLASTn and Diamond (BLASTx), respectively, were used to annotate assembled contigs (Buchfink, Xie, and Huson 2015). Contigs were categorised into higher kingdoms using the BLASTn ‘sskingdoms’ flag option. Non-viral blast hits including host contigs with sequence similarity to viral transcripts (e.g. endogenous viral elements) were removed from further analysis during manual screening. A maximum expected value of 1 × 10⁻¹⁰ was used as a cut-off to filter putative viral contigs. Viral contigs that have previously been identified as viral contaminants from laboratory components were also removed from further analysis (Asplund et al. 2019). Based on the BLASTn and Diamond results (database accessed November 2022), putative viral contigs were further analysed using Geneious Prime 2022.2.2 to find and translate open reading frames (ORFs).

2.5. Protein structure similarity searching for viral discovery

Similar to previous work (Waller et al. 2022), we used a protein structure similarity search to identify highly divergent viral transcripts that did not share significant amino acid sequence similarity to other known transcripts. Such ‘orphan contigs’ (Ortiz-Baez et al. 2020) were translated into ORFs using the EMBOSS getorf program (Rice, Longden, and Bleasby 2000), with the minimum nucleotide size of the ORF set to 1,000 nucleotides, the maximum nucleotide size of the ORF set to 50,000, and the ‘methionine’ flag option set to only report ORFs with a methionine amino acid starting codon. Reported ORFs were submitted to Phyre2, which uses remote homology detection to build 3D models to predict and analyse protein structure and function (Kelley et al. 2015). Virus sequences with predicted polymerase structures with a confidence value of >90%  were aligned with representative amino acid sequences from the same viral family obtained from NCBI RefSeq using MAFFT v7.490 (L-INS-I algorithm). Conserved domains were visually confirmed before phylogenetic trees were estimated using the same method outlined in section 2.7.

2.6. Estimating viral transcript abundance estimations

Viral abundances were estimated using Trinity’s ‘align and estimate’ tool. RNA-Seq by Expectation–Maximization (RSEM) (Li and Dewey 2011) was selected as the method of abundance estimation, Bowtie2 (Langmead and Salzberg 2012) as the alignment method, and the ‘prep reference’ flag enabled. To mitigate the impact of contamination due to index-hopping, viral transcripts with expected abundances of less than 0.1% of the highest expected abundance for that virus across other libraries were removed from further analysis. Total viral abundance estimates for viruses from vertebrate hosts (i.e. bats), vertebrate-associated hosts (i.e. where the true host of the virus remained unclear), or viruses from both vertebrate hosts and vertebrate-associated hosts across viral families and orders were compiled across all libraries. Estimated abundances were standardised to the number of paired reads per library.

2.7. Virus phylogenetic analysis

Translated viral protein polymerase sequences (i.e. RNA-dependent RNA polymerase, RdRp) or capsid sequences were aligned with representative protein sequences from the same taxonomic viral family or order obtained from NCBI RefSeq as well as the closest BLAST hit using MAFFT v7.450 (L-INS-I algorithm) (see Supplementary Table 2 for lengths of sequence alignments) (Katoh et al. 2002). Poorly aligned regions were removed using trimAL v1.2rev59 with the gap threshold flag set to 0.9 (Capella-Gutierrez, Silla-Martinez, and Gabaldon 2009). IQ-TREE v1.6.12 was used to estimate maximum likelihood phylogenetic trees for each viral species/family/order (Nguyen et al. 2015). The LG amino acid substitution model was selected with 1,000 ultra-fast bootstrapping replicates and the -alrt flag specifying 1,000 bootstrap replicates for the SH-like approximate likelihood ratio test. Phylogenetic trees were annotated using Figtree v1.4.4 (Rambaut 2023). To compare genes within alphacoronaviruses, phylogenetic trees were first estimated for both the RdRp (ORF1b) as well as the spike protein, and a tanglegram was created using the ape (Paradis, Claude, and Strimmer 2004) and dendextend (Galili 2015) packages with the tanglegram function in RStudio v4.3.1.

2.8. Analysis of alpha diversity on virome composition

All statistical analysis plots were created using RStudio v4.3.1 with the tidyverse ggplot2 package (Wickham et al. 2019). Using the diversity analysis function which is part of the vegan package (Oksanen et al. 2020), the Shannon index was selected as the index method to measure alpha diversity of standardised abundance estimates of viral families/orders across bat species and location (North Island or South Island). A Welches t-test was used to determine whether there was a significant difference in virome alpha diversity across bat species and location (North Island and South Island). The cor and cor.test function in R was used to test the correlations between the pairwise differences in Shannon index and the pairwise distance between roost site (km).

2.9. Sequencing depth across viral transcripts

Bowtie2 was used to map raw reads against viral transcripts (Langmead and Salzberg 2012). Read depth at each nucleotide position across the viral transcripts was measured using samtools depth function (Li et al. 2009). Plots showing read depth at each nucleotide position across the viral transcripts were made in RStudio v4.3.1 with the tidyverse ggplot2 package (Wickham et al. 2019).

2.10. Viral nomenclature

A virus was arbitrarily considered a novel species if it shared <90% amino acid similarity within the most conserved region (i.e. RdRp/polymerase or capsid) (King, MJ, and EB 2011; Wang et al. 2017), unless otherwise stated. For novel virus sequences, we have provided a proposed virus name (prior to formal verification by the International Committee on Taxonomy of Viruses (ICTV)). Viruses that were most closely related to a virus identified in a vertebrate metagenome sample, such that their true host was uncertain, as well as invertebrate host-associated viruses, were classified as ‘bat metagenome derived’ viruses as these viruses were unlikely to directly infect bats. We have also proposed the general name, pekapeka, where viruses were found in both bat species (Māori refer to New Zealand’s bats generally as pekapeka).

3. Results

Total RNA isolated from 219 bat guano samples were pooled into 12 libraries based on bat species, sampling location/roost, and sampling year and were subject to whole-RNA sequencing (Supplementary Table 1). The number of sequencing reads generated from the 12 bat guano pools varied between 58-105 million paired-end reads per library (Fig. 1b, Supplementary Table 1). As some of the bat roosts were in close proximity (Fig. 1a), we further confirmed the host assignment of guano samples. To do so, we analysed the standardised abundance of M. tuberculata and C. tuberculatus contigs within each library (Fig. 1c). While a single library assigned as long-tailed bats (B12) contained contigs matching lesser short-tailed bats, their overall abundance was just 0.11%. All other libraries contained contigs from the correct host. Viral contigs contributed 0.02–18.9 % of the total standardised abundances across each of the bat guano pooled samples, while contigs that could not be assigned to a source based on sequence similarity (i.e. orphan contigs) made up 1.75–60.24 % of the total standardised abundances (Fig. 1d).

3.1 Viral abundance and diversity

Analysis of bat guano samples revealed that the viromes of long-tailed bats were more diverse than those of lesser short-tailed bats, with samples of the former containing bat viral contigs spanning seven viral families (Picornaviridae, Astroviridae, Paramyxoviridae, Papillomaviridae, Poxviridae, Rhabdoviridae, and Coronaviridae) (Fig. 2). Surprisingly, no bat-specific viral contigs were identified in the lesser short-tailed bat samples in this study, although short transcripts of bat papillomaviruses, coronaviruses, hepeviruses, polyomaviruses, and caliciviruses have previously been uncovered in lesser short-tailed bats in New Zealand using Illumina MiSeq2000 and PCR approaches (Wang et al. 2015; Tortosa et al. 2023). Both bat species examined here contained bat ‘metagenome derived’ viral contigs (Fig. 2). These represent viral contigs in which the true host is uncertain but were likely dietary- or environmental-associated viruses. Long-tailed bat samples contained metagenome-derived viral contigs spanning the Caliciviridae and Adenoviridae families, and data from both bat species contained metagenome-derived viral contigs from Basavirus sp., Flaviviridae, Hepeviridae, Nodaviridae, Peribunyaviridae, Picobirnavirdae, and Reovirales (Fig. 2).

3.2 Novel bat viruses

3.2.1 Coronaviridae

We identified several alphacoronavirus full genomes from long-tailed bat guano sampled from the Grand Canyon Cave (Fig. 3a) and identified two species of alphacoronaviruses. One alphacoronavirus ORF1b transcript, which contains the RdRp, identified here shared 90.22% amino acid similarity with Alphacoronavirus sp. WA2028 from Gould’s wattled bats from Australia (QGX41956.1), suggesting that the same virus species infects hosts from Australia and New Zealand based on ICTV species definitions (Fig. 3a, Supplementary Table 2) (King, MJ, and EB 2011; Prada et al. 2019). We termed this virus Pekapeka alphacoronavirus 2. Phylogenetic analysis of the ORF1b region of the Pekapeka alphacoronavirus 2 revealed that this transcript was also likely the same viral species as two previously identified viruses, Mystacina coronavirus New Zealand/2013 and Chalinolobus tuberculatus alphacoronavirus, although these previously published viral sequences were very short in length (186 and 129 amino acids, respectively) and hence were not a significant match due to inadequate (approximately 10%) query coverage. Pekapeka alphacoronavirus 2 shared amino acid sequence similarity of 98% and 100% with Mystacina coronavirus New Zealand/2013 and Chalinolobus tuberculatus alphacoronavirus, respectively. The former virus was previously uncovered in lesser short-tailed bats sampled from the remote, offshore New Zealand island, Whenua Hou in 2013 (Hall et al. 2014), while the latter was uncovered and classified in long-tailed bats from the Grand Canyon Cave and sampled in 2020 (Fig. 3a) (Tortosa et al. 2023).

A full genome assembled for Pekapeka alphacoronavirus 2 revealed that this genome was markedly short in length, comprising only 22,416 nucleotides (Fig. 4). Upon further investigation of its genome organisation, the short coronavirus, which was relatively highly abundant with a standardised abundance of 6.5×10⁻⁴, did not encode the non-structural proteins (NSP) 1 and 2, while NSP3 appeared to be truncated in comparison to other bat alphacoronaviruses (Fig. 4). Additionally, while ORF1b shared >90% amino acid similarity with ORF1ab from Alphacoronavirus sp. WA2028 (QGX41956.1) previously identified in a Gould’s wattled bat from Australia (Prada et al. 2019), other ORFs shared only 32–65% amino acid sequence similarity. Given this large deletion in ORF1a, we suggest that this coronavirus be classified as a new species, Pekapeka alphacoronavirus 2 (Supplementary Fig. 2 shows a graph depicting the read depth mapped across the Pekapeka alphacoronavirus 2 genome).

We identified two further alphacoronavirus contigs from the B5 and B6 libraries. These contigs shared >90%amino acid sequence similarity with each other indicating that they are the same viral species. We termed this virus Pekapeka alphacoronavirus 1. Both Pekapeka alphacoronavirus 1 sequences shared 87% amino acid sequence similarity with the closest known genetic relative, Alphacoronavirus sp., WA1087 again identified in Gould’s wattled bat in Australia (QGX41944.1) (Prada et al. 2019) (Fig. 3a, Supplementary Table 2). Phylogenetic analysis showed that two previously identified viruses, Chalinolobus tuberculatus alphacoronavirus in long-tailed bats (Tortosa et al. 2023) and Mystacina tuberculata alphacoronavirus in lesser short-tailed bats (Tortosa et al. 2023), were likely the same viral species. These previously identified virus sequences were short in length (133 amino acids) compared to the full Pekapeka alphacoronavirus 1 genome identified here in the B6 library. While these short sequences shared 100% amino acid identity to Pekapeka alphacoronavirus 1, the query coverage of just 5% meant they were not considered significant matches. Nevertheless, assuming these are the same viral species, we propose the revised name, Pekapeka alphacoronavirus 1, reflecting that this virus infects both New Zealand’s bat species.

To further explore these newly discovered alphacoronaviruses, we investigated the evolutionary history of their ORF1b and spike proteins through a cophylogenetic analysis. Interestingly, the spike proteins of Pekapeka alphacoronavirus 1 and 2 shared >99% amino acid sequence similarity, reflecting very recent common ancestry, despite only sharing 84% amino acid sequence similarity within the ORF1b polymerase (Fig. 4b). Of more note, rather than falling as sister-group to Pekapeka alphacoronavirus 1, Pekapeka alphacoronavirus 2 (B5) shared 98% amino acid sequence similarity within ORF1b to Mystacina coronavirus New Zealand/2013, although these two viruses shared only 38% amino acid sequence similarity in their spike proteins (Fig. 3a, 4b). Such incongruence between the spike and polymerase phylogenies indicates that Pekapeka alphacoronavirus 2 (B5) has a recombinant evolutionary history (Fig. 4b).

3.2.2 Astroviridae

Three novel viral polymerase contigs were identified in long-tailed bat guano sampled from the Grand Canyon Cave that shared sequence similarity to viruses within the Astroviridae (positive-sense, single stranded (+ssRNA)) (Fig. 3b). Two of these contigs shared >90% amino acid sequence similarity with each other and spanned a similar region of the polymerase protein, indicating that these contigs were a singular viral species, tentatively named Chalinolobus tuberculatus astrovirus 1 (Supplementary Table 2). Chalinolobus tuberculatus astrovirus 1 contigs fell within the genus Mamastroviruses (Fig. 3b) and shared 71% amino acid sequence similarity with the closest known genetic relative, Bat astrovirus BtSY2 (WBV74323.1), which was identified in pomona roundleaf bat (Hipposideros pomona) in China (Wang et al. 2023). The third astrovirus contig shared 78.3% amino acid similarity with the closest known genetic relative, Bat astrovirus (QTE76055.1), previously identified in a common vampire bat (Desmodus rotundus) in Peru (Bergner et al. 2021).

3.2.3 Paramyxoviridae

Six Paramyxoviridae contigs were uncovered in long-tailed bat guano metatranscriptomes sampled from the Grand Canyon Cave. All six Paramyxoviridae contigs identified here fell within the subfamily Orthoparamyxovirinae (Fig. 3c). Two of these contigs shared >90% amino acid sequence similarity with each other and were provisionally named Chalinolobus tuberculatus paramyxovirus 1 (Fig. 3c). Chalinolobus tuberculatus paramyxovirus 1 shared approximately 73% amino acid sequence similarity with the closest known genetic relatives, Bat paramyxovirus (AIF74192.1) and Miniopterus schreibersii paramyxovirus (AGU69461.1), which were previously uncovered in greater tube-nosed bats (Murina leucogaster) and common bent-wing bats (Miniopterus schreibersii) from China (Supplementary Table 2) (Wu et al. 2016).

Another contig that was identified here shared 86% amino acid sequence similarity with the closest known genetic relative, Bat paramyxovirus (AIF74185.1), which was uncovered in Daubenton’s bats (Myotis daubentonii) from China (Wu et al. 2016). This contig was provisionally named Chalinolobus tuberculatus paramyxovirus 2. A further two contigs shared greater than 90% similarity with each other and were consequently named Chalinolobus tuberculatus paramyxovirus 3. This virus shared approximately 87% amino acid similarity with the closest known genetic relatives, Paramyxovirus PREDICT_PMV-52 (ALS55487.1) and Bat paramyxovirus (AIF74192.1), previously uncovered in western bent-winged bats (Miniopterus magnater) from Thailand and greater tube-nosed bats from China, respectively (Wu et al. 2016). The final transcript belonging to the Paramyxoviridae we uncovered shared 92% amino acid sequence similarity with Bat paramyxovirus, which was previously uncovered in greater tube-nosed bats from China (Wu et al. 2016). The close genetic match suggests that the same virus species is also present in New Zealand.

3.2.4 Papillomaviridae

A viral transcript sharing amino acid sequence similarity to the Papillomaviridae (a double-stranded DNA (dsDNA) virus) was identified in bat guano from long-tailed bats sampled from Mackay Creek in the Eglington Valley. This contig, termed Chalinolobus tuberculatus papillomavirus, shared 88.46% amino acid sequence similarity with the closest known genetic relative, Eptesicus serotinus papillomavirus 2, which was identified in serotine bats (E. serotinus) from Spain (Fig. 3d, Supplementary Table 2) (García-Pérez et al. 2014). The novel Chalinolobus tuberculatus papillomavirus along with Eptesicus serotinus papillomavirus 2 fell within the Firstpapillomavirinae subfamily.

3.2.5 Rhabdoviridae

We identified a contig from the Rhabdoviridae (negative-sense, single stranded (-ssRNA)) in long-tailed bat guano sampled from Grand Canyon Cave. This viral contig shared 62.63% amino acid sequence similarity with the closest known genetic relative, Wenzhou Rhinolophus pusillus ledantevirus 1 sampled from least horseshoe bats (Rhinolophus pusillus) in China (Supplementary Table 2). Termed Chalinolobus tuberculatus rhabdovirus; this virus fell within the Alpharabdovirinae subfamily alongside other species from the Alpharabdovirinae identified in Chinese rufous horseshoe bats (R. sinicus) (Luo et al. 2021), as well as eloquent horseshoe bats (R. eloquens) and striped leaf-nosed bats (Macronycteris vittatus) from Kenya (Fig. 3e) (Kading et al. 2013; Walker et al. 2015).

3.2.6 Picornaviridae

Two contigs from the Picornaviridae were identified in bat guano from long-tailed bats sampled from the Grand Canyon Cave. Both contigs shared 100% amino acid sequence similarity with each other and were provisionally named Chalinolobus tuberculatus picornavirus. The Chalinolobus tuberculatus picornavirus fell within the Paavivirnae subfamily alongside the closest known genetic relatives, Parechovirus sp. QAPp32 (amino acid sequence similarity 61.21%) identified in the intestinal sample of the common pipistrelle bat (Pipistrellus pipistrellus) from China, as well as Ferret parechovirus (amino acid sequence similarity 55.46%) identified in a rectal swab from ferrets in the Netherlands (Smits et al. 2013) (Fig. 3f, Supplementary Table 2).

3.3 Novel bat metagenome-derived viruses

All of the viruses described in section 3.3 (Fig. 5) are likely dietary, commensal, or environmental in origin as they were most closely related to metagenome derived viruses or invertebrate viruses.

3.3.1 Basavirus sp

Three Basavirus contigs were identified in bat guano from lesser short-tailed bats sampled from Pureora Forest Park, two of which shared 100% amino acid sequence similarity and were termed Mystacina tuberculata metagenome-derived basavirus 1 (Fig. 5a, Supplementary Table 2). The Mystacina tuberculata metagenome-derived basavirus 1 was most closely related to Picornavirales sp., previously identified in a soil sample from China (Chen et al. 2022) (amino acid sequence similarity 63%) and Fisavirus 1, previously identified in the intestinal contents of a common carp (Cyprinus carpio) from Hungary (amino acid sequence similarity 59%) (Reuter et al. 2015). The third contig shared 56% amino acid sequence similarity with the closest genetic relative, Basavirus sp., which was identified in the guano of a lesser Asiatic yellow bat (Scotophilus kuhlii) from Vietnam (Fig. 5a, Supplementary Table 2) (Oude Munnink et al. 2017). Although these viruses were related to other bat viruses, the true host of these basavirus species is unclear since such divergent members of the Picornavirales are often found within the fecal viromes of vertebrates (Oude Munnink et al. 2017).

3.3.2 Flaviviridae

We identified four flaviviruses provisionally named, Pekapeka metagenome-derived flavivirus, and Chalinolobus tuberculatus metagenome-derived flaviviruses 1–3 (Fig. 5b, Supplementary Table 2). Viral contigs from Pekapeka metagenome-derived flavivirus were identified in both lesser short-tailed and long-tailed bat guano samples (Fig. 5b, Supplementary Table 2). As a result, we termed this virus Pekapeka metagenome-derived flavivirus reflecting that this virus was identified in both of New Zealand’s bat species. Pekapeka metagenome-derived flavivirus fell within the unclassified pesti-large genome flavivirus clade and was highly divergent as indicated by its long branch length on the phylogenetic tree (Fig. 5b). The closest BLASTp hit to the Pekapeka metagenome-derived flavivirus was to a Hymenopteran flavi-related virus, identified in a species of bee (Thyreus orbatus) from Italy (amino acid sequence similarity 22–23%) (Paraskevopoulou et al. 2021) (Supplementary Table 2).

Chalinolobus tuberculatus metagenome-derived flavivirus 1 also shared 23% amino acid sequence similarity with the closest genetic relative, Hymenopteran flavi-related virus (Paraskevopoulou et al. 2021) (Supplementary Table 2). Chalinolobus tuberculatus metagenome-derived flavivirus 1 fell within the pesti-large genome flaviviruses clade infecting invertebrate hosts (Fig. 5b). Two additional viral transcripts, Chalinolobus tuberculatus flaviviruses 2 and 3, shared 66.2% sequence similarity and 65.91% sequence similarity with their closest their genetic relatives and fell alongside members of the Flavivirus genera that had previously been identified in invertebrate species in Turkey and Senegal (Öncü et al. 2018; Gil et al. 2021) (Fig. 5b, Supplementary Table 2).

3.3.3 Hepeviridae

Six novel hepeviruses were identified in the lesser short-tailed and long-tailed bat guano metatranscriptomes and were termed Mystacina tuberculata metagenome-derived hepevirus 1–4 and Chalinolobus tuberculatus metagenome-derived hepevirus 1–2 (Fig. 5c, Supplementary Table 2). All six novel hepeviruses described shared <90% amino acid sequence similarity with their closest known genetic relatives, indicating that these viruses were novel (Supplementary Table 2). The closest genetic relatives to these viruses were to avian and reptilian metagenome-associated viruses (Truchado et al. 2020; French et al. 2022; Shan et al. 2022; Waller et al. 2022) (Fig. 5c).

A further hepevirus was discovered which shared >90% amino acid sequence similarity with the closest known genetic relative, Tuatara cloaca-associated hepevirus 1, previously discovered in cloacal swabs from tuatara (Sphenodon punctatus) (Waller et al. 2022).

3.3.4 Nodaviridae

We identified eight nodaviruses within the bat guano metatranscriptomes of both long-tailed bats and lesser short-tailed bats (Fig. 5d). One of these viruses had previously been discovered, sharing >96% amino acid sequence similarity with Bat nodavirus (Dacheux et al. 2014). This virus had been identified in the brain of a Serotine bat but was classified as an invertebrate host-related virus (Supplementary Table 2). The remaining seven nodaviruses were classified as novel and were provisionally named Mystacina tuberculata metagenome derived nodavirus 1–5 and Chalinolobus tuberculatus metagenome derived nodavirus 1–2. Their closest known genetic relatives were a number of vertebrate host metagenome-associated viruses (Dacheux et al. 2014; Geoghegan et al. 2021; Chen et al. 2022; Shan et al. 2022; Zhu et al. 2022).

3.3.5 Peribunyaviridae

Two peribunyaviruses were uncovered in the bat guano metatranscriptomes of long-tailed bats sampled from the Grand Canyon Cave and lesser short-tailed bats sampled from Pureora Forest Park (Fig. 5e, Supplementary Table 2). The novel viruses were provisionally termed Mystacina tuberculata metagenome-derived peribunyavirus and Chalinolobus tuberculatus metagenome-derived peribunyavirus. The two Chalinolobus tuberculatus metagenome-derived peribunyavirus contigs identified in this study shared 54% and 41% amino acid sequence similarity with their closest known genetic relatives, Hubei insect virus 1 and Gonipterus platensis bunyan-like virus, respectively, previously identified in arthopods and a Gum-tree weevil (Gonipterus platensis) (Shi et al. 2016; Silva et al. 2020) (Supplementary Table 2).

A protein structural similarity-based approach was used to identify the highly divergent viral contig termed Mystacina tuberculata metagenome-derived peribunyavirus. The top PHYRE2 structural match was to a La Crosse virus RNA-directed RNA polymerase l, a member of the Peribunyaviridae family, modelled with a confidence of 92%, a percentage identity of 28% and a coverage of 39%. Upon analysing the amino acid sequence of Mystacina tuberculata metagenome-derived peribunyavirus along with other reference Peribunyaviridae polymerase sequences, including the La Crosse virus, the Mystacina tuberculata metagenome-derived peribunyavirus shared multiple conserved residues within the highly conserved RdRp motifs B–E, further supporting the view that this contig is viral (Fig. 6). Upon phylogenetic analysis, Mystacina tuberculata metagenome-derived peribunyavirus fell within a sister clade to Soft tick bunyavirus, which was previously identified in soft ticks (Argas vespertilionis) collected from bat guano in Japan (Oba et al. 2016) (Fig. 5e).

3.3.6 Picobirnaviridae

Nine novel viruses within the Picobirnaviridae were identified in the bat guano metatranscriptomes of lesser short-tailed bats sampled from Pureora Forest Park and were termed Mystacina tuberculata metagenome-derived picobirnavirus 1–4 (Fig. 5f, Supplementary Table 2). The closest known genetic relatives were viruses sampled from vertebrate host metagenomes (Luo et al. 2018; Ramesh et al. 2021; Chen et al. 2022).

An additional two novel viruses were identified in the bat guano metatranscriptomes of long-tailed bats sampled from the Grand Canyon Cave and were named Chalinolobus tuberculatus metagenome-derived picobirnavirus 1 and 2 (Fig. 5f). Their closest genetic relative was Limbe picobirna-like virus (amino acid sequence similarity approximately 56%) previously identified in a straw-coloured fruit bat (Eidolon helvum) but was classified as an invertebrate host virus (Yinda et al. 2018) (Supplementary Table 2).

3.3.7 Reovirales

We identified five novel reoviruses within the bat guano metatranscriptomes of lesser short-tailed bats sampled from Pureora Forest Park and long-tailed bats sampled from the Grand Canyon Cave (Fig. 7). All five novel reoviruses-termed Mystacina tuberculata metagenome-derived reovirus 1–4 and Chalinolobus tuberculatus metagenome-derived reovirus are likely members of the Spinareoviridae (Fig. 7). We also uncovered a viral contig that shared 97.2% amino acid sequence similarity with Avian associated reo-like virus 3 previously identified in the metagenome of a New Zealand fantail (Rhipidura fuliginosa) (Supplementary Table 2) (French et al. 2022).

3.4 Factors shaping the diversity of bat viromes

We next asked whether the alpha diversity of bat viromes was influenced by bat species and/or sampling location (North or South Island). This analysis revealed that neither bat species nor location significantly influenced alpha diversity of bat viromes (P = 0.2 and 0.6, respectively) (Fig. 8a and b). Even though long-tailed bats harboured bat viruses across seven viral families (Fig. 2), all but one of these viruses was found across just two libraries sampled from the Grand Canyon Cave. Therefore, the absence of bat-specific viruses within other libraries rendered this comparison statistically insignificant. Additionally, there was a weak positive yet not significant correlation between the pairwise differences in Shannon index of bat viruses across bat libraries and the pairwise distance (in kilometers) between roost sites (Pearsons correlation = 0.19, P = 0.13) (Fig. 8c).

Similarly, when we considered bat viromes comprised of bat viruses as well as bat metagenome-derived viruses, neither bat species nor location influenced the alpha diversity (P = 0.3 and 0.1, respectively) (Supplementary Fig. 1a). In contrast, when considering only the metagenome-derived virome, alpha diversity was significantly higher in lesser short-tailed bats compared to long-tailed bats (P = 0.04) (Supplementary Fig. 1b), likely because lesser short-tailed bats consume a broader diet than long-tailed bats and may therefore be exposed to a wider range of dietary-related viruses (Daniel 1976; O’Donnell et al. 2021). Similarly, bats that were sampled from the North Island were significantly more diverse in their metagenome-derived viromes in comparison to bats that were sampled from the South Island (P = 0.02) (Supplementary Fig. 1b). However, the majority of samples from the North Island were from lesser short-tailed bats (four out of six) in comparison to only one (out of four) sampled from the South Island.

4. Discussion

We analysed the viromes of Aotearoa New Zealand’s only native terrestrial mammalian species, the lesser short-tailed bat and the long-tailed bat, to determine patterns of cross-species virus transmission and whether their contrasting host ecologies and life histories have influenced their virome composition. Our analysis resulted in the discovery of many bat viruses in guano samples from long-tailed bats, yet no bona fide bat viruses in lesser short-tailed bats were identified. At face value, this suggests that long-tailed bats may be a more important reservoir compared to lesser short-tailed bats. Notably, viruses spanning the families Picornaviridae, Astroviridae, Paramyxoviridae, Papillomaviridae, Poxviridae, Rhabdoviridae, and Coronaviridae were identified in guano metatranscriptomes of long-tailed bats, including, to our knowledge, the shortest known coronavirus genome documented to date. These novel bat viruses were all closely related to other previously identified bat viruses. Prior to this study, only short transcripts from two viruses, both alphacoronaviruses, had been identified in long-tailed bat guano (Tortosa et al. 2023). Consequently, this study has expanded our knowledge of the viruses that circulate within New Zealand’s long-tailed bat populations.

The Grand Canyon Cave, home to a population of long-tailed bats located in the central North Island of New Zealand, was where all but one of the novel bat viruses were discovered. Long-tailed bats often move between roosts and primarily roost during the day in small colonies distributed among multiple tree-roosts in native forests (O’Donnell 2000), yet are also known to roost at night in large numbers in the Grand Canyon Cave, particularly during spring (O’Donnell 2002b; O’Donnell et al. 2021). The maximum number of night roosting long-tailed bats that were counted during any one time in the Grand Canyon Cave was 358, while a total of 533 individuals were banded over three seasons indicating that bats from multiple different day roosting colonies may come together to roost in the cave at night (O’Donnell 2002b). While the reason behind night roosting in long-tailed bats is unclear (O’Donnell 2002b), the Grand Canyon Cave likely provides an ideal site to facilitate viral transmission and thus high viral richness as demonstrated here. Examining the role that bat roosting behaviour plays in the transmission of viruses within New Zealand’s bat species would add to a better understanding of the potential risk of viral spillover.

In previous studies, short virus transcripts from the Coronaviridae, Papillomaviridae, Picornaviridae, Polyomaviridae, Caliciviridae, and Hepeviridae have been uncovered in guano from lesser short-tailed bats located on Whenua Hou (Codfish Island), a small (14 km²), predator-free, offshore island in southern New Zealand (Hall et al. 2014; Wang et al. 2015). Whenua Hou is home to a natural population of more than 2,000 short-tailed bats that occupy more than 120 tree roosting sites across the island (Sedgeley 2006; O’Donnell 2021). While we also identified novel viral contigs belonging to the Hepeviridae in this host, phylogenetic analysis placed these contigs as most likely to be dietary, commensal, or environmental in origin as they were closely related to viruses identified in the metagenomes of avian and reptilian species. Here, we were unable to detect any bat viruses in guano sampled from mainland populations of lesser short-tailed bats.

We identified two alphacoronaviruses in long-tailed bat guano that were closely related to alphacoronaviruses previously identified in lesser short-tailed bats. Although there is no evidence that lesser short-tailed bats and long-tailed bats occupy the same roosts (Sedgeley 2003), roost sites can be within close proximity where both species may interact (O’Donnell 2021; O’Donnell et al. 2021). That both bat species harbour the same viral species, such that these viruses have jumped between these hosts, provides evidence of this interaction. The large geographical distance between where these viruses where sampled (>1,000 km) indicates that these viruses may be present across multiple roosts throughout New Zealand, although we were unable to detect them across other locations sampled here. Longitudinal sampling of bats across New Zealand may shed new light on the persistence of viruses within these populations, particularly as the rate of viral shedding amongst bats is known to vary throughout the year (Joffrin et al. 2022).

In this study, bat guano was sampled during a single time point, either January 2020 or January–February 2021, during which time both long-tailed and short-tailed bats produce offspring (O’Donnell 2002a; Sedgeley 2006). Communal roosts of long-tailed bats during this time typically consist of breeding females and their pups while adult males and non-breeding females make up on average 15% and 22% of the inhabitants, respectively (O’Donnell and Sedgeley 1999; O’Donnell 2002a). Similarly, female short-tailed bats form a communal maternity roost during this time; however, males and non-breeding females are also known to use communal maternity roosts as day roosts (O’Donnell et al. 2021). During the lactation period in late summer, male short-tailed bats are also known to gather around the maternity roost, a behaviour known as lek breeding (Lloyd 2001). Consequently, the changes in communal roost composition, temporal variations in bat activity, and changes in viral shedding over time were not assessed here. Longitudinal sample collection would be useful to elucidate how these unique breeding behaviours and social interactions impact virome composition.

Long-tailed and lesser short-tailed bats are highly divergent, belonging to different families. As such, the incongruence observed between viral and host topologies supports the cross-species transmission of coronaviruses. Both the long-tailed bat alphacoronaviruses identified here were closely related to those present in Gould’s wattled bat from Australia (Prada et al. 2019). Like New Zealand’s long-tailed bats, these bats are members of the Vespertilionidae. It is likely, then, that long-tailed bats were reservoirs of coronaviruses before migrating to New Zealand and transmitted these viruses to lesser short-tailed bats since their arrival (Dool et al. 2016; O’Donnell et al. 2021). This new evidence suggests that estimates of the time to most recent common ancestor of orthocoronaviruses may require further calibration (Hayman and Knox 2021), although such analyses are challenged by uncertainty over the timing of bat migrations to New Zealand (approximately 1–17 million years ago (mya) for long-tailed bats and approximately 16–35 mya for lesser short-tailed bats) estimated from fossil records and species divergence times (Worthy and Holdaway 1994; Hand et al. 1998; Den Bussche RA and Hoofer 2000; Dool et al. 2016; O’Donnell et al. 2021). Clearly, obtaining full-length genomes of alphacoronaviruses previously identified in lesser short-tailed bats is an important step to arriving at a more comprehensive understanding of their evolutionary history.

While coronaviruses typically range from 26 to 32 kb in length (Lu et al. 2020), we identified an alphacoronavirus of only 22,416 nucleotides, containing a large deletion in ORF1a of both NSP1 and NSP2, as well as a truncated NSP3. NSP1 is involved in inhibiting the host immune response, decreasing the overall expression of genes in the host cells while redirecting host machinery to produce viral proteins (Yuan et al. 2021). While gammacoronaviruses and deltacoronaviruses are known to lack NSP1, all currently known alpha and betacoronaviruses contain NSP1 (Papineau et al. 2019; Yuan et al. 2021; Bermudez, Miles, and Muller 2023). The function of NSP2 is currently unclear (Gupta et al. 2021), although deletion of NSP2 in two other betacoronaviruses—Severe acute respiratory syndrome coronavirus 1 and Murine hepatitis virus—still resulted in viable viruses in cell culture, although viral growth and RNA synthesis was reduced (Graham et al. 2005). NSP3 is involved in a number of roles including polyprotein processing, acting as a scaffold protein and is an essential component of the replication and transcription complex (Lei, Kusov, and Hilgenfeld 2018). To the best of our knowledge, this alphacoronavirus may be the shortest coronavirus genome identified to date. Despite its sequence similarity within the RdRp, its shortened genome and divergence in other genes suggest that it should be classified as a novel virus species.

We inferred cophylogenies to compare the evolutionary histories of the alphacoronaviruses genes identified here. Notably, spike proteins from Pekapeka alphacoronavirus 1 and 2, both found in long-tailed bats, shared >99% amino acid similarity despite being different viral species (i.e. <90% amino acid sequence similarity within the RdRp) indicative of relatively recent common ancestry. In contrast, in the ORF1b polymerase, Pekapeka alphacoronavirus 2 was more closely related (98% amino acid identity) to Mystacina coronavirus New Zealand/2013 previously identified in short-tailed bats than to Pekapeka alphacoronavirus 1. Hence, the evolution of these alphacoronaviruses in New Zealand reveals the presence of both cross-species transmission among long- and short-tailed bats, and recombination.

Two novel viruses from the Flaviviridae were uncovered in the guano of both bat species. Flaviviruses have been known to infect a wide range of hosts and cause disease in humans, livestock, and wildlife (Simmonds et al. 2017). Both novel viruses identified here fell within the ‘pesti-large genome’ clade (Mifsud et al. 2023). The pesti-large genome clade is associated with invertebrate hosts (Mifsud et al. 2023), making it likely that these viruses are dietary in origin (Daniel 1976; O’Donnell et al. 2021). One virus, Pekapeka metagenome-derived flavivirus, was especially divergent, expanding the diversity of this clade and highlighting the potential to uncover highly divergent novel viruses within hosts that have been severely understudied.

In sum, this study has increased our understanding of the viruses that circulate within New Zealand’s bat populations and indicates that coronaviruses have likely jumped between these bat populations since long-tailed bats arrived at least two million years ago. A future focus on understanding the influence of roosting behaviour and seasonality on virome composition and further sampling New Zealand’s mammalian hosts will be important to expand our knowledge of viral diversity in New Zealand bats and determine whether the bat viruses documented are present in other terrestrial mammalian hosts.

Data availability

The raw sequencing reads generated in this project are available on the Aotearoa Genomic Data Repository, DOI number https://doi.org/10.57748/cyhd-ad62 while the virus sequences have been submitted to GenBank and the accession numbers can be found in Supplementary Table S2. Alignments and code for the statistical analysis can be found at https://github.com/stephwaller/NZ-Bat-Virome-Paper.git.

Supplementary data

Supplementary data is available at VEVOLU Journal online.

Acknowledgements

We would like to thank the Department of Conservation for coordinating the sampling of the bat guano from multiple sites across New Zealand. We would also like to thank local iwi for supporting our research.

Funding

This project was funded by a New Zealand Royal Society Rutherford Discovery Fellowship (RDF-20-UOO-007) awarded to J.L.G. and a University of Otago Doctoral Scholarship awarded to S.J.W.

Conflict of interest:

None declared.

References

Author notes