As part of a long-term study characterizing the seasonal dynamics of potentially zoonotic viruses in wild fruit bats in Madagascar, monthly captures of Malagasy pteropodid bats were carried out at species-specific roost sites in the Districts of Moramanga and Manjakandriana, Madagascar between 2018 and 2019 (P. rufus: Ambakoana roost, −18.513 S, 48.167 E; E. dupreanum: Angavobe cave, −18.944 S, 47.949 E; Angavokely cave = −18.933 S, 47.758 E; R. madagascariensis: Maromizaha cave, −18.9623 S, 48.4525 E). In brief, bats were captured in nets hung in the tree canopy (P. rufus) or over cave mouths (E. dupreanum, R. madagascariensis) at dusk (17:00–22:00) and dawn (03:00–07:00), removed from nets, and processed under manual restraint following methods that have been previously described (61, 76, 77). All animals were identified to species, sex, and age class (juvenile vs. adult), and fecal, throat, and urine swabs were taken from each individual, collected into viral transport medium, and frozen on site in liquid nitrogen. Post-sampling, swabs were transported to −80°C freezers for long-term storage in the Virology Unit at Institut Pasteur de Madagascar. In total, 2156 bats (P. rufus: 167 females, 184 males; E. dupreanum: 495 females, 421 males; R. madagascariensis: 416 females, 473 males) were captured across the course of our long-term field study, though only fecal samples from a subset of individuals (see ‘Results’) were analyzed in part with the coronavirus research outlined here. Of those captures, 84 P. rufus were juveniles and 267 adults, 108 E. dupreanum were juveniles and 810 adults, and 126 R. madagascariensis were juveniles and 767 adults. Additional details relevant to our time series are included in the results.

This study was carried out in strict accordance with research permits obtained from the Madagascar Ministry of Forest and the Environment (permit numbers 019/18, 170/18, 007/19) and under guidelines posted by the American Veterinary Medical Association. All field protocols employed were pre-approved by the UC Berkeley Animal Care and Use Committee (ACUC Protocol # AUP-2017-10-10393), and every effort was made to minimize discomfort to animals.

RNA was extracted from a randomly selected subset of fecal, throat, and urine swab samples in the Virology Unit at the Institut Pasteur de Madagascar, with each sample corresponding to a unique individual from the field dataset. Samples undergoing mNGS corresponded to individuals captured in Feb-Apr, Jul-Sep and December 2018 or in January 2019. Water controls were extracted in conjunction with samples on each unique extraction day. Extractions were conducted using the Zymo Quick DNA/RNA Microprep Plus kit (Zymo Research, Irvine, CA, USA), according to the manufacturer's instructions and including the step for DNAse digestion. Post-extraction, RNA quality was checked on a nanodrop to ensure that all samples demonstrated 260/280 ratios exceeding 2 and revealed quantifiable concentrations. Resulting extractions were stored in freezers at −80°C, then transported on dry ice to the Chan Zuckerberg Biohub (San Francisco, CA, USA) for library preparation and metagenomic Next Generation Sequencing (mNGS).

Four randomly selected samples from each of three bat species underwent additional quantification using an Invitrogen Qubit 3.0 Fluorometer and the Qubit RNA HS Assay Kit (ThermoFisher Scientific, Carlsbad, CA, USA). After quantification, all total RNA samples, along with water samples from Madagascar extractions, were manually arrayed into 96 well plates to automate high throughput mNGS library preparation. Based on the initial quantitation, a 2 μL aliquot from each plated sample was diluted 1:9 on a Bravo liquid handling platform (Agilent, Santa Clara, CA, USA). A 5 μL aliquot from each diluted sample was arrayed into a 384 well plate for input into the mNGS library prep. Samples derived from fecal, throat, and urine swab samples were arrayed on distinct 384 well plates for separate sequencing runs. Additional unrelated total RNA samples (a dilution series of total RNA isolated from cultured HeLa cells) and a set of local lab water samples were included on each 384 well plate to serve as library preparation controls. Input RNA samples in the 384 well plate were transferred to a GeneVac EV-2 (SP Industries, Warminster, PA, USA) to evaporate the samples to enable miniaturized mNGS library preparation with the NEBNext Ultra II RNA Library Prep Kit (New England BioLabs, Beverly, MA, USA). Library preparation was performed per the manufacturer's instructions, with the following modifications: 25 pg of External RNA Controls Consortium Spike-in mix (ERCCS, Thermo-Fisher) was added to each sample prior to RNA fragmentation; the input RNA mixture was fragmented for 8 min at 94°C prior to reverse transcription; and a total of 14 cycles of PCR with dual-indexed TruSeq adapters was applied to amplify the resulting individual libraries. An initial equivolume library pool was generated, and the quality and quantity of that pool was assessed via electrophoresis (High-Sensitivity DNA Kit and Agilent Bioanalyzer; Agilent Technologies, Santa Clara, CA, USA), real-time quantitative polymerase chain reaction (qPCR) (KAPA Library Quantification Kit; Kapa Biosystems, Wilmington, MA, USA), and small-scale sequencing (2 x146 bp) on an iSeq platform (Illumina, San Diego, CA, US). Subsequent equimolar pooling of individual libraries from each plate was performed prior to performing large-scale paired-end sequencing (2 × 146 bp) run on the Illumina NovaSeq sequencing system (Illumina, San Diego, CA, USA). The pipeline used to separate the sequencing output of the individual libraries into FASTQ files of 146bp paired-end reads is available on GitHub at https://github.com/czbiohub/utilities.

Raw reads from Illumina sequencing were host-filtered, quality-filtered, and assembled on the CZID (v3.10, NR/NT 2019-12-01) platform, a cloud-based, open-source bioinformatics platform designed for microbe detection from metagenomic data (78), using a host background model of “bat” compiled from all publicly available full-length bat genomes in GenBank at the time of sequencing. Samples were deemed “positive” for coronavirus infection if CZID successfully assembled at least two contigs with an average read depth >2 reads/nucleotide that showed significant nucleotide or protein BLAST alignment(s) (alignment length >100 nt/aa and E-value <0.00001 for nucleotide BLAST/ bit score >100 for protein BLAST) to any CoV reference present in NCBI NR/NT database (version 12-01-2019). To verify that no positives were missed from CZID, all non-host contigs assembled in CZID underwent directed, offline BLASTn and BLASTx (79) against a reference database constructed from all available full-length nucleotide and protein reference sequences for Alpha- and Betacoronavirus available in NCBI Virus (last access: August 15, 2021). Step-by-step instructions for our offline BLAST protocol can be accessed in our publicly available GitHub repository at: https://github.com/brooklabteam/Mada-Bat-CoV/.

Three full genome-length Nobecovirus contigs returned from CZID (two from R. madagascariensis and one from P. rufus) were aligned with Nobecovirus homologs from NCBI (see ‘Phylogenetic Analysis’) and annotated in the program Geneious Prime (2020.0.5). We then used NCBI BLAST and BLASTx to query identity of our full length recovered genomes and their respective translated proteins to publicly available sequences in NCBI (79). We queried identity to reference sequences for four previously described Nobecovirus strains (accession numbers: MG762674 (Rousettus bat coronavirus HKU9), NC_030886 (Rousettus bat coronavirus RoBat-CoV GCCDC1), MK211379 (Rhinolophus affinis coronavirus BtRt-BetaCoV/GX2018), and NC_048212 (Eidolon helvum bat coronavirus), as well as to the top BLAST hit overall. Finally, we aligned representative sequences from each major Nobecovirus clade and visually examined the region of p10 orthoreovirus insertion from the RoBat-CoV GCCDC1 lineage in the newly described sequences from Madagascar.

Contigs returned from CZID were combined with publicly available coronavirus sequences in NCBI to perform phylogenetic analysis. We carried out four major phylogenetic analyses, building (a) a full-genome Betacoronavirus maximum likelihood (ML) phylogeny, (b) a Betacoronavirus ML phylogeny corresponding to a conserved 259 bp fragment of the RNA-dependent RNA polymerase (RdRp) gene encapsulated in the CoV Orf1b, (c) four amino acid ML phylogenies derived from translated nucleotides corresponding to the spike (S), envelope (E), matrix (M), and nucleocapsid (N) proteins from a subset of full length genomes, and (d) a time-resolved Bayesian phylogeny corresponding to all available full genome Nobecovirus sequences in NCBI virus, which we computed in BEAST2 (80). Detailed methods for the construction of each phylogeny are available at https://github.com/brooklabteam/Mada-Bat-CoV/.

Briefly, our full genome ML phylogeny was comprised of 122 unique NCBI records, corresponding to all available full genome sequences with bat hosts under NCBI taxon IDs, Betacoronavirus (694002), unclassified Betacoronavirus (696098), Betacoronavirus sp. (1928434), unclassified Coronaviridae (1986197), or unclassified Coronavirinae (2664420) (107 records), in addition to all full genome Betacoronavirus (694002) reference sequences with a non-bat host (14 records), plus one Gammacoronavirus outgroup (accession number NC_010800). The full genome phylogeny additionally included three full length Madagascar Nobecovirus sequences returned from CZID (two from R. madagascariensis and one from P. rufus), which are described in this paper for the first time, for a total of 125 unique sequences.

Our Betacoronavirus RdRp ML phylogeny consisted of an overlapping subset of a 259 bp fragment in the center of the RdRp gene that has been previously described in Madagascar fruit bats (12) (7 records), in addition to the same RdRp fragment extracted from near-full length Nobecovirus sequences on NCBI Virus (17 records) and full length reference sequences for other Betacoronavirus subgenera available in NCBI Virus (17 records). This phylogeny also included two NCBI Virus RdRp Nobecovirus fragments, in addition to seven Madagascar Nobecovirus sequences encompassing the RdRp fragment of interest, which were returned from the assembly in CZID (four from R. madagascariensis, two from P. rufus, and one from E. dupreanum). Finally, we included the RdRp fragment of our Gammacoronavirus outgroup, for a total of 51 unique sequences.

Our amino acid phylogenies consisted of S, E, M, and N gene extractions from a subset of the same representative set of near-full genome length sequences used in the RdRp analysis: the same 17 full-length Betacoronavirus reference sequences, 17 near full-length Nobecovirus sequences, and the one Gammacoronavirus outgroup, in addition to our three full genome Madagascar sequences derived from R. madagascariensis and P. rufus. Gene extractions were carried out using annotation tracks reported with each accession number in NCBI or, in cases where annotations were unavailable, genes were manually annotated and extracted in Geneious Prime based on alignment to homologs. After nucleotide extraction, genes were translated prior to alignment.

Finally, our Bayesian time-resolved Nobecovirus phylogeny consisted of all available full genome Nobecovirus sequences in NCBI virus; because recombination events can muddle molecular clock estimation in phylogenetics, we constructed two different versions of this timetree, one including 18 Nobecoviruses sequences only for which no past history of recombination has been described, and a second which added two additional sequences from the Nobecovirus GCCDC1 clade known for its p10 orthoreovirus insertion.

After compiling sequences for each phylogenetic analysis, sequence subsets for the full-length, RdRp, and four amino acid phylogenies were aligned in MAFFT v.7 (81, 82) using default parameter values. Alignments were checked manually for quality in Geneious Prime, and the RdRp alignment was trimmed a fragment (259 bp) conserved across all sequences in the subset. All sequence subsets and alignment files are available for public access in our GitHub repository: https://github.com/brooklabteam/Mada-Bat-CoV/.

After quality control, alignments were sent to Modeltest-NG (83) to assess the best fit nucleotide or amino acid substitution model appropriate for the data. All alignments for ML analysis (full genome, 259 bp RdRp fragment, and amino acid protein sequences) were then sent to RAxML-NG (84) to construct the corresponding phylogenetic trees. Following best practices outlined in the RAxML-NG manual, twenty ML inferences were made on each original alignment and bootstrap replicate trees were inferred using Felsenstein's method (85), with the MRE-based bootstopping test applied after every 50 replicates (86). Bootstrapping was terminated once diagnostic statistics dropped below the threshold value and support values were drawn on the best-scoring tree.

We constructed the Bayesian timetree using the Bayesian Skyline Coalescent (87) model in BEAST2 (80), assuming a constant population prior, and the best fit nucleotide substitution model as indicated by ModelTest-NG. Sampling dates corresponded to collection date as reported in NCBI Virus; in cases where only year was reported, we assumed a collection date of 15-July for the corresponding year. We tested trees using both an uncorrelated exponentially distributed relaxed molecular clock (UCED) and a strict clock but ultimately reported results from the strict clock assumptions, as similar results were inferred from both. Markov chain Monte Carlo (MCMC) sample chains were run for 600 million iterations, convergence was checked using TRACER v1.7 (88), and trees were averaged after 10% burn-in using TreeAnnotater v2.6.3 (89) to visualize mean posterior densities at each node.

All resulting phylogenies (both ML and Bayesian) were visualized in R v.4.0.3 for MacIntosh, using the package ‘ggtree’ (90).

Full length Nobecovirus sequences derived from CZID were analyzed for any signature of past recombination. First, the ORF1a, ORF1b, S, NS3, E, M, N, and NS7 genes from the P. rufus Nobecovirus sequence, the longest R. madagascariensis Nobecovirus sequence (OK067320), and two full genome representative sequences from the HKU9 (NC_009021) and E. helvum African lineages (NC_048212) were extracted, translated, and concatenated. Concatenated, translated sequences were then aligned in MAFFT v.7 (81, 82) using default parameter values. Nobecovirus sequences corresponding to the RoBat-CoV GCCDC1 (27, 65) and BtRt-BetaCoV/GX2018/BtCoV92 (62, 63) genotypes were omitted from recombination analyses because inserted genes and/or genetic material upstream from the nucleocapsid in the corresponding genomes interfered with the alignment.

After alignment, genomes were analyzed for amino acid similarity in the program pySimplot (91), using the P. rufus and, subsequently, the R. madagascariensis genome as query sequences, the HKU9 and Eidolon helvum African Nobecovirus clades as references, and the corresponding Madagascar sequence as the alternative. Analyses were carried out using a window size of 100 aa and a step size of 20 aa.

Next, all three full length nucleotide sequences of Madagascar Nobecovirus genomes were aligned with grouped full genome sequences corresponding to the two disparate Nobecovirus lineages: the HKU9 lineage (EF065514-EF065516, HM211098-HM211100, MG693170, NC_009021, MG762674) and the E. helvum African lineage (MG693169, MG693171-MG693172, NC_048212). As before, alignment was conducted in MAFFT v.7 (81, 82) using default parameter values.

After alignment, genomes were analyzed for recombination in the program SimPlot (v.3.5.1). Nucleotide similarity plots were generated using the P. rufus and, subsequently, the R. madagascariensis genomes as query sequences, the HKU9 and Eidolon helvum African Nobecovirus clades as references, and the corresponding Madagascar sequence as the alternative. Bootscan analyses were conducted on the same alignment, using the same query and reference inputs. Both nucleotide similarity and Bootscan analyses were carried out using a window size of 200 bp and a step size of 20 bp.

All three annotated full-length genome sequences (two from R. madagascariensis, one from P. rufus), plus four additional RdRp gene fragment sequences (two from R. madagascariensis, one from P. rufus, and one from E. dupreanum) were submitted to NCBI and assigned accession numbers OK020086-OK020089 (RdRp fragments) and OK067319-OK067321 (full genomes).

RNA from 285 fecal, 143 throat, and 196 urine swab samples was prepped into libraries and submitted for Illumina sequencing. In 28/285 (9.82%) fecal specimens and in 2/196 (1.00%) urine specimens, at least two contigs with an average read depth >2 reads/nucleotide, and nucleotide or protein-BLAST alignments to any CoV reference sequence in NCBI were identified via CZID analysis. Because the prevalence detected in the urine samples was low, it is likely attributable to field contamination with fecal excrement upon urine swab collection, as bats often excrete both substances simultaneously under manual restraint. None of the 143 throat swabs assayed demonstrated evidence of CoV infection.

Prevalence in feces varied slightly across species, with 4/44 (9.1%) P. rufus specimens, 16/145 (11.0%) E. dupreanum specimens, and 8/96 (8.3%) R. madagascariensis specimens sequencing CoV positive. Juveniles demonstrated higher CoV prevalence than adults for P. rufus and E. dupreanum but not for R. madagascariensis. Juvenile vs. adult prevalence was 3/15 (20%) vs. 1/29 (3.5%) for P. rufus, 5/13 (38.5%) vs. 11/132 (8.3%) for E. dupreanum, and 0/13 (0%) vs. 8/83 (9.6%) for R. madagascariensis (Figure 1A). Prevalence varied seasonally across all three species, peaking coincidentally in adult and juvenile populations for P. rufus and E. dupreanum, with the highest prevalence for all three species observed during the wet season months of February-April when late-stage juveniles are present in the population, following each species' annual birth pulse (Figure 1B).

Three full or near-full CoV genome length contigs were recovered from CZID for Nobecoviruses derived from R. madagascariensis (two genomes: 28,980 and 28,926 bps in length) and P. rufus (one genome: 29,122 bps in length). In all three genomes, we successfully identified ORF1ab (including RdRp) and structural proteins S (spike), E (envelope), M (matrix), and N (nucleocapsid), in addition to accessory genes NS3, NS7a, and NS7b (Figure 2A). In keeping with convention outlined in (65), the accessory genes, NS7a and NS7b, were so named based on nucleotide alignment and amino acid identity to homologous proteins in previously described Nobecoviruses.

BLAST analysis of the full genome indicated that the P. rufus Nobecovirus sequence is highly divergent, demonstrating only 72.87–73.54% identity to all previously described Nobecovirus clades, with the top blast association to E. helvum Nobecovirus lineages (Supplementary Table 1). Additionally, Nobecovirus genomes derived from R. madagascariensis demonstrated high identity (~95%) to E. helvum Nobecovirus lineages circulating in Africa. BLASTx analysis of individual genes from viruses derived from both Madagascar species demonstrated the highest identity with previously described Nobecovirus sequences in the Orf1ab region (which includes RdRp) for both P. rufus and R. madagascariensis viruses (76.1% identity for P. rufus Nobecovirus to E. helvum bat coronavirus and ~99% identity for R. madagascariensis Nobecovirus to E. helvum bat coronavirus). By contrast, both P. rufus and R. madagascariensis Nobecovirus genomes demonstrated substantial divergence from all known homologs in the S gene, where only 46.93–86.9% identity was observed. The P. rufus Nobecovirus was similarly divergent in the N gene, though R. madagascariensis Nobecoviruses demonstrated high (~91%) identity to CoV genotypes from E. helvum in this region.

In general, BLASTx queries of NS7 accessory proteins in both R. madagascarienis and P. rufus Nobecovirus demonstrated ~40–91% amino acid identity to already-characterized Nobecovirus proteins (Supplementary Table 1). Genomes derived from R. madagascariensis appeared slightly more complex than those derived from P. rufus, allowing for annotation of one additional accessory gene, NS7c, which has been characterized in recombinant Nobecovirus sequences of the RoBat-CoV GCCDC1 lineage (27, 65). Curiously, BLASTx query of the NS7a accessory protein in the P. rufus genome showed no identity to any previously described Nobecovirus protein; rather, the highest scoring protein alignment (31.25% identity, 1e-06 E-value) of the NS7a translation encompassed 40% of the query (query coverage was located at the 3' end of the query length) and corresponded to an arachnid Low-Density Lipoprotein Receptor-Related Protein 1 (LRP-1) (Supplementary Table 2). As LRP-1 is involved in the mammalian innate immune response (92), we hypothesized that this putative novel ORF could be a viral gene involved in immune antagonism. To check the integrity of our de novo assembly in NS7a, we mapped the deduplicated raw reads from mNGS to the full genome P. rufus Nobecovirus contig generated by CZID (78) (Supplementary Figure 1). We confirmed >200x read coverage across the region corresponding to the putative NS7a accessory protein, with good representation of both forward and reverse-facing reads across the length of the protein, as well as the intergenic regions preceding and succeeding it. We were also able to identify a putative Transcription Regulatory Sequence (TRS) preceding this gene (Table 1), further validating our confidence that P. rufus Nobecovirus NS7a represents a real though highly divergent protein.

The p10 orthoreovirus insertion within the RoBat-CoV GCCDC1 Nobecovirus lineage was not observed in either Nobecovirus genomes from R. madagascariensis or P. rufus. Nonetheless, examination of the multiple sequence alignment of representative sequences of all Nobecovirus clades in this region demonstrated the presence of some variable genetic material downstream from the N gene and upstream from the NS7a gene in the divergent P. rufus Nobecovirus genome (Figure 2B). Nobecoviruses clustering in the BtRt-BetaCoV/GX2018 - BtCoV92 lineage also carry a unique coding sequence in this region, highlighting the dynamic nature of the 3' end of the CoV genome (63).

In addition to the identification of both canonical and novel ORFs described above, we also observed non-coding TRS elements preceding all the major proteins in all three Nobecovirus genomes (Table 1). Many of these correspond to the 5'-ACGAAC-3' six bp core motif common to many Betacoronaviruses, including SARS-CoV and previously described in Nobecoviruses of the GCCDC1 and GX2018/BtCoV92 lineages (62, 65, 93). For most genes, these TRS elements were located a short distance upstream from the corresponding gene (Table 1). Elements identified in the two R. madagascariensis genomes were largely comparable, suggesting that these two sequences could represent slight variations in the same virus lineage. Some putative TRS elements, including that preceding P. rufus NS7a, showed variation from the 5'-ACGAAC-3' core motif, with some recapitulating the 5'-AAGAA-3' motif common to SARS-CoV-2 (94). TRS variations may be indicative of variation in gene expression across individual bats and/or species.

Phylogenetic analysis of full length Betacoronavirus genomes confirmed that both P. rufus and R. madagascariensis genomes cluster in the Nobecovirus subgenus of the Betacoronaviruses, with the divergent P. rufus forming its own distinct clade and both R. madagascariensis genomes grouping with the previously described E. helvum reference sequence from Cameroon (64) (Figure 3A). We observed distinct groupings of five main Nobecovirus lineages in our phylogeny: (a) the largely Asian-derived HKU9 sequences, (b) the African E. helvum-derived sequences (now including new R. madagascariensis Nobecovirus genomes), (c) the recombinant GCCDC1 genomes, (d) the BtRt-BetaCoV/GX2018 and BtCoV92 genomes described respectively from China and Singapore, and (e) the divergent P. rufus genome contributed here from Madagascar. Intriguingly, the P. rufus genome groups ancestral to all other Nobecoviruses, followed by the E. helvum/R. madagascariensis African lineage, with the Asian genotypes forming three distinct (and more recent) clades corresponding to genotypes HKU9, GCCDC1, and GX2018 – BtCoV92. This finding suggests that the evolutionary origins of the Nobecovirus clade likely predate the recent divergence of P. rufus from sister Pteropus spp. taxa (though see results for Bayesian timetree below), which are primarily Asian in their distribution. To date, no other Nobecovirus sequences have yet been described from any Pteropus spp. bats; additional sampling across other Pteropus lineages will be needed to confirm our hypothesis of a Pteropus genus origin for this Betacoronavirus clade. Further phylogenetic analysis of a 259 bp fragment of the RdRp gene reconfirmed these groupings and suggested the presence of at least two distinct genetic variants within the P. rufus lineage (Figure 3B). One RdRp fragment derived from feces of the third Malagasy fruit bat, E. dupreanum, grouped within the E. helvum – R. madagascariensis African Nobecovirus lineage, consistent with previous reporting (12). Characterization of the full length genome of this virus will be needed to clarify whether it represents a genetic variant of or a distinct genotype from the R. madagascariensis virus. Phylogenetic analysis of the RdRp fragment allowed for inclusion of one partial Nobecovirus sequence derived from E. helvum bats in Kenya (HQ728482), which also grouped within the E. helvum – R. madagascariensis African clade, confirming the distribution of this genotype across West and East Africa and into the South-Western Indian Ocean Islands. Notably, one partial Cameroonian E. helvum sequence (MG693170) clustered with HKU9 sequences from Asia, rather than within the E. helvum – R. madagascariensis African clade. These findings suggest that both “African” and “Asian” Nobecovirus lineages are likely broadly geographically distributed.

Amino acid phylogenies computed from translated protein alignments of the S, E, M, and N Betacoronavirus structural genes (Supplementary Figure 2) further confirmed evolutionary relationships suggested in Figure 3. S, M, and N gene phylogenies demonstrated distinct groupings of five main Nobecovirus lineages outlined above, while in the E gene phylogeny, the P. rufus sequence grouped adjacent to the single Cameroonian-derived E. helvum sequence within the HKU9 clade. However, we note that bootstrap values were extremely low in this E gene phylogeny, suggesting that additional sampling is needed to confirm these evolutionary relationships.

Results from our Bayesian timetree analysis (Figure 4; Supplementary Figure 3) recapitulated ML support for the five distinct Nobecovirus subclades but indicated a time to MRCA for the entire Nobecovirus clade of ~300 years ago (1695; 95% HPD: 1643–1748) (Figure 4), with the African Eidolon lineage (which included our R. madagascariensis sequences) branching off ~200 years ago (1827; 95% HPD: 1797–1857). Surprisingly, both of these viral divergence times post-date the estimated radiation of Madagascar host bats from sister species (58, 59), suggesting that Malagasy bat populations may support substantial viral exchange with bats from surrounding islands and the African continent. Additional sampling of Nobecovirus lineages in Asian-distributed Pteropus spp. bats will be crucial to ascertaining whether the evolutionary origins of the Nobecovirus subgenus could date back any further than that suggested by samples highlighted here. Evolutionary relationships were preserved in phylogenies which included GCCDC1 orthoreovirus recombinant sequences, though inclusion of these sequences diminished certainty surrounding evolutionary divergence times, resulting in broad HPD distributions (Supplementary Figure 3).

SimPlot analysis confirmed the evolutionary distinctiveness of the P. rufus Nobecovirus genome, which showed <70% amino acid similarity and <50% nucleotide similarity to HKU9, E. helvum, and R. madagascariensis genotypes across the majority of its genome length (Figures 5A,B). Consistent with BLAST results, the P. rufus Nobecovirus genome demonstrated the highest similarity to previously described sequences in the Orf1b region, which includes RdRp. The R. madagascariensis Nobecoviruses, by contrast, showed >90% amino acid and nucleotide similarity to the E. helvum African lineage throughout Orf1ab, but both P. rufus and R. madagascariensis sequences diverged from all other reference genomes in the first half of the spike protein, which corresponds to the S1 subunit and includes the receptor binding domain that mediates viral entry into host cells (95). Further divergence for both P. rufus and R. madagascariensis Nobecoviruses was observed in the N structural protein and in the NS7 accessory genes. Bootscan analysis further confirmed these findings, showing that the P. rufus Nobecovirus clusters with HKU9 lineages across Orf1ab, NS3, E, and M genes but demonstrates evidence of recombination with E. helvum – R. madagascariensis African lineages in the S (particularly S1), N, and NS7 genes (Figure 5C). Similarly, bootscanning demonstrated that R. madagascariensis Nobecoviruses group with the E. helvum lineage across Orf1ab, NS3, E, and M but show evidence of recombination with HKU9 and P. rufus Nobecovirus in S (again, particularly S1), N, and NS7 genes (Figure 5C), thus highlighting the dynamic nature of these regions of the Nobecovirus genome.