The first complete mitochondrial genome sequence of Cynopterusbrachyotis (Chiroptera, Pteropodidae) from the Philippines

Abstract The technical limitations of capillary sequencing in providing insights on phylogeny have been greatly aided in recent years by the implementation of next generation sequencing platforms which can generate whole mitochondrial genome (mitogenome) sequences. In this study, enriched mitochondrial DNA of Cynopterusbrachyotis from Mindanao, Philippines was sequenced using the Illumina MiSeq platform. A total of 653,967 clean paired-end reads was assembled using a MIRA-MITObim pipeline, resulting in a consensus mitogenome sequence length of 17,382 bases and a GC content of 41.48%, which is consistent with other published mitogenomes in fruit bats. The assembled C.brachyotis mitogenome was annotated using the MITOS online server and was able to resolve all mitochondrial genes, except for one transfer RNA gene (trnT) which may be further resolved by additional capillary sequencing of the region. Sequence analysis showed that the Philippine C.brachyotis is only 90%-91% homologous with other Cynopterus spp., based on its full mitogenome sequence. Phylogenetic analysis of fruit bat mitogenomes, deposited in online repositories, revealed that the Philippine C.brachyotis in this study has diverged from Asian Cynopterus, namely Cynopterusbrachyotis and Cynopterussphinx from other parts of Asia (100% bootstrap support) with the latter two forming a separate clade. This divergence at the species level was consistent with phylogentic inference using cytochrome oxidase 1 (CO1) and cytochrome B (cytb) gene markers. Our results strengthen the previously reported hypothesis that the Cynopteruscf.brachyotis in the Philippines is distinct from its Asian counterparts and should be, therefore, treated as a new species.


Introduction
Cynopterus brachyotis (Müller, 1838) is a fruit bat widespread in Southeast Asia, ranging from Sri Lanka, Taiwan, Vietnam, Indonesia, Malaysia and the Philippines (Jayaraj et al. 2012). The species is believed to have reached the Philippine islands within the past 1 million years (Rowe 2018). Currently, it is abundant in different types of habitats in the country with differing degrees of disturbances, ranging from urban residential, agricultural to forested areas (Rowe 2018). Due to high population abundance, widespread geographical extent and the relatively stable population, the species is considered of least concern under the updated IUCN Red List (Csorba et al. 2019, Rosell-Ambal et al. 2019).
The status of C. brachyotis taxonomy in the Philippines is uncertain and challenging. Genetic studies suggest that the Philippine population has high levels of divergence in comparison to other C. brachyotis populations in other Southeast Asian countries (Bacus et al. 2021, Luczon et al. 2019. Moreover, other authorities suggested that the Philippine Cynopterus should be assigned as Cynopterus luzoniensis, a taxon believed to exist throughout the Philippine archipelago and in the Sulawesi Island of Indonesia (Campbell et al. 2004, Rosell-Ambal et al. 2019, Rowe 2018, Simmons 2005. However, these hypotheses should be further evaluated as the exact identification of C. luzoniensis in Sulawesi may be different from "C. luzoniensis" in the Philippines or in some islands in the country (Campbell et al. 2004, Jayaraj et al. 2012. Due to this confusion on the distribution and extent of the taxa mentioned above, we will be consistently using "Cynopterus cf. brachyotis" in this paper to represent the Philippine Cynopterus to avoid confusion until taxonomy is resolved. The analyses of mitochondrial genes are crucial to understand the challenging and confusing taxonomy and evolutionary history of bats in the Philippines. However, information on the mitochondrial genome (mitogenome) of Philippine bats is clearly fragmented despite the wide application of this data for nomenclature, evolution and phylogeny of many species. Fortunately, the technical limitations of capillary sequencing to provide phylogenetic evidence have been greatly aided by the use of Next Generation Sequencing (NGS) platforms that can generate whole mitogenome sequences.
As part of the larger initiative to sequence the mitochondrial DNA of native species in the country, we provide the complete mitochondrial genome (mitogenome) of Cynopterus cf. brachyotis from Mindanao Island. Moreover, phylogenetic trees, based on the sequences generated from this study and from those deposited in GenBank, were provided to show the evolutionary relationship of the Philippine population to other Cynopterus spp. from other countries in Asia. Results from this study will contribute to the taxonomic status and identity of the Philippine Cynopterus which, in turn, can be used as additional supporting information for other phylogenetic, evolution and ecological studies of bat species in the world.

Samples
The sample used in this study was a cardiac tissue from a voucher specimen, previously identified as Cynopterus brachyotis using morphological characteristics and measurements (e.g. forearm length, head length, body length, tail length and ear length) with the aid of the taxonomic key developed by Ingle and Heaney (1992). Moreover, a ~ 202 bp barcode for the mitochondrial cytochrome oxidase 1 (CO1) gene (Walker et al. 2016) was determined in a separate study (Bacus et al. 2021) that further confirmed the identity of the specimen. An IACUC approval from the University of the Philippines Manila (Protocol No.: 2018-109) and a Gratuitous permit from the Philippine Department of Environment and Natural Resources Region XI (WGP No.: 2018-07) were obtained for the collection of voucher specimens.

Library Preparation and Illumina Sequencing
The two PCR amplicons were quantified using Qubit dsDNA HS assay and normalised in a resuspension buffer to a final volume of 2.5 µl and a final concentration of 10 ng/µl (25 ng total yield). Both amplicons were pooled in a new PCR tube. DNA tagmentation, amplification and library clean-up were performed using the Nextera XT Library Prep Kit according to the manufacturer's instructions. Cleaned libraries were then sent to the Philippine Genome Center Diliman-DNA Sequencing Core Facility for normalisation and next-generation sequencing using the Illumina MiSeq platform and a paired-end read format at 2 x 150 bp for 500 cycles.

Mitogenome Assembly and Analysis
Raw sequencing reads were subjected to initial quality assessment with the FastQC software (Andrews 2010). Raw reads were trimmed of adapters and low-quality bases at both ends using Trimmomatic (Bolger et al. 2014). The filtered reads were used to generate an initial mapping assembly using MIRA ) and were iteratively mapped back to the initial mapping assembly and assembled using MITObim (Hahn et al. 2013). The reference sequence that was used for the assembly was the NCBI Reference Sequence entry for Cynopterus brachyotis mitochondrion, complete genome (GenBank Accession ID: NC_026465). The different genes/regions in the resulting assembly were annotated using the MITOS web server (Bernt et al. 2013). Both the entire assembled genome and the identified genes/regions were annotated individually using the nucleotide BLAST tool at the NCBI website (NCBI Resource Coordinators et al. 2018). Default settings were used in the parameters for running all of the aforementioned procedures. The assembled mitogenome sequence was submitted to NCBI GenBank and subsequently assigned the accession number OK163841.

Phylogenetic Analysis
The assembled mitogenome sequence of Cynopterus brachyotis from Mindanao, Philippines was aligned with available fruit bat mitogenome sequences in the NCBI GenBank. The best phylogenetic model was determined using MEGA 7, based on the lowest BIC and AIC values (Kumar et al. 2016). A phylogenetic tree was subsequently inferred using the GTR+G+I DNA substitution and site heterogeneity model using a Maximum-Likelihood estimation with 1000 bootstrap replicates in MEGA 7 (Kumar et al. 2016). A similar analysis was performed for the full cytochrome oxidase 1 (CO1) and cytochrome B (cytb) genes. For the CO1 gene, a multiple sequence alignment of 1,533 bp was used for the phylogenetic analysis using the Maximum-Likelihood method with the HKY+G+I model in MEGA 7 with 1000 bootstrap replicates. For the cytb gene, a multiple sequence alignment of 1,134 bp was used for the phylogenetic analysis using the Maximum-Likelihood method with the GTR+G+I model in MEGA 7 with 1000 bootstrap replicates. A distant vertebrate group, fish, was selected as the outgroup.

Mitogenome
A total of 2,180,947 reads was initially generated by next generation sequencing of the mitochondrial DNA-enriched tissue sample (Table 1). The quality control steps reduced this number of reads to 1,784,834. From this readpool, a total of 653,967 reads was used for the assembly of the enriched mitochondrial DNA by the MIRA-MITObim workflow. The generated mitochondrial DNA assembly had a total length of 17,382 bases with a GC content of 41.48%. The average depth of coverage of this assembly was calculated to be 10,158.30.  (Table 2). Further annotation of the assembled sequence by MITOS identified a total of 36 genes: 21 tRNA genes, 2 rRNA genes and 13 protein-coding genes (Fig. 1). All known vertebrate mitochondrial genes have been found in this mitogenome, except trnT. Performing online nucleotide BLAST on individual gene regions identified by MITOS resulted in mostly either C. brachyotis or C. sphinx top hits with percent identities within the range of 90%-97%. It is also noted that the top hit for trnM(atg) has no significant similarity found with any NCBI-BLAST database entry, a few top hits are non-Cynopterus related entries (e.g. Megaerops niphanae and Myonycteris torquata) and the top hits for trnY(tac) and trnD(gac) are non-bat related entries, such as Maxomys surifer and Equus asinus, respectively (Table 3). BLAST  Summary of mitochondrial genome assembly statistics. Table 2.
Top five NCBI-BLAST hits of the assembled mitogenome sequence.

Phylogenetics
A phylogenetic tree, based on the complete mitogenome sequences of selected fruit bats, was generated (Fig. 2). The cynopterine fruit bats formed a monophyletic clade with 100% bootstrap support. Interestingly, the Cynopterus brachyotis from Mindanao, Philippines appeared to have diverged from the rest of the Asian Cynopterus ( C. brachyotis & C. sphinx) with a robust bootstrap support (100%). This branching pattern also appeared to be consistent with phylogenetic trees generated using the two commonly-used markers, CO1 and cytb genes (Fig. 3).

Discussion
An almost complete mitogenome (99.86%, contains a single 25-base "N" gap in the final assembled sequence) was generated using an Illumina NGS and MIRA-MITObim assembly workflow. Nucleotide BLAST of the assembled mitogenome and individual genes confirmed its identity to be Cynopterus sp., the first reported mitogenome of this genus from the Philippines. Further annotation of individual mitochondrial genes by MITOS determined that almost all vertebrate mitochondrial genes were present on the assembled mitogenome, except for trnT. The length (17,382 bases) and GC content (41.48%) of the assembled mitogenome were similar and hence concordant with other published fruit bat mitogenomes (Jahari et al. 2020, Jahari et al. 2021, Yue 2019. These steps and resulting statistics confirm that a true, high-quality bat mitogenome has been assembled from the sequence reads and the identity of this assembled mitogenome can be attributed to Cynopterus sp.. However, the relatively low similarity of the whole mitogenome (90.86% and 90.56% to C. sphinx and C. brachyotis, respectively) and individual gene regions (90%-97%) to both C. brachyotis and C. sphinx entries suggest that the assembled mitogenome originated from a species that is genetically unique from the aforementioned. The presence of non-Cynopterus and non-bat BLAST hit results and even a hit result with no significant similarity, for a number of the mitochondrial genes due to low similarity to Cynopterus-related database entries, provide further evidence to this proposition.
Phylogenetic analysis of the mitogenomes and individual CO1 and cytb genes confirmed that the sequenced bat in this study formed a genetically-distinct clade from the Asian C. brachyotis and C. sphinx. Similar findings have been reported in recent years for the Philippine Cynopterus cf. brachyotis through sequence analysis of the standard (500 bp) and even shorter (200 bp Cynopterus brachyotis is hypothesised to be a species complex, based on sequence analysis of the mitochondrial cytb and control region sequences (Campbell et al. 2004). In the Philippines, two distinct lineages are present (Campbell et al. 2004), suggesting multiple colonisation events from multiple origins. The Sunda C. brachyotis lineage is present in Palawan Island and the rest of the Philippines form a distinct lineage which is proposed for taxonomic reassignment to "C. luzoniensis". However, this reassignment is not consistently accepted by some authors. Thus, we recommend further investigation using more genetic samples from different islands in the Philippines that includes representatives from the Sulawesi population to fully resolve the taxonomic issues revolving around the Philippine Cynopterus. Moreover, inconsistencies on the taxonomic identity involving this species complex suggests the importance of including morphological and ecological data in resolving taxonomic problems (Campbell et al. 2004).
Phylogenetic analysis, using both whole mitogenome and individual gene markers, such as CO1 and cytb, produced similar phylogenetic trees with high branch support values, suggesting that these gene markers can be used to reliably resolve taxonomic groups in fruit bats. While this holds true at the genus level, careful consideration should be given when using single gene markers for species delineation as exemplified by the discordance between the CO1 and cytb phylogenetic trees for the placement of C. sphinx within the cynopterine clade. Such contrasting branching patterns might arise especially when only a few representative samples are used in the analysis. This study focused on the mitochondrial genome due to its utility in taxonomic resolution. Mitochondrial DNA evolves faster compared to nuclear genetic markers, generating a higher degree of sequence variation which makes it useful to resolve lower taxonomic levels for organisms (Chan et al. 2021, Hassanin et al. 2013. Furthermore, the relatively small size of the mitogenome (~ 16,000 base pairs long) makes it a practical alternative as a molecular marker. Despite these characteristics, the sole use of the mitochondrial genome in taxonomy remains to have certain limitations. Alternatively, nuclear markers represent neutral and paternal inheritance. In particular, single nucleotide polymorphisms (SNPs) and microsatellites exhibit high polymorphism amongst populations or individuals (Galbusera et al. 2000, Jiang et al. 2016). In addition, the highly-conserved nature of nuclear rRNA gene sequences makes it suitable for resolving higher taxonomic levels for organisms (Chan et al. 2021). Other studies have shown that combining mitochondrial and nuclear DNA analysis provided better taxonomic resolution (Abreu et al. 2020, Lambret-Frotté et al. 2012. A similar approach can be used in future studies to provide a more comprehensive and accurate perspective to ascertain the taxonomic resolution of the Mindanao or Philippine Cynopterus cf. brachyotis. Maximum Likelihood phylogeny of fruit bat A) CO1 and B) cytb gene obtained from GenBank. The tree was generated using MEGA 7 with 1000 bootstrap replicates. For the CO1 gene, a multiple sequence alignment of 1,533 bp was used for the phylogenetic analysis using the Maximum-Likelihood method with the HKY+G+I model in MEGA 7 with 1000 bootstrap replicates. For the cytb gene, a multiple sequence alignment of 1,134 bp was used for the phylogenetic analysis using the Maximum-Likelihood method with the GTR+G+I model in MEGA 7 with 1000 bootstrap replicates. The scale bar denotes the number of substitutions per site which is reflected in the branch lengths (i.e. 0.5 value means 5 nucleotide changes for every 10 nucleotides). Bootstrap values for clade support are written on the nodes of the branches. An outgroup consisting of fish mitogenomes was included in the phylogenetic tree inference.
The first complete mitochondrial genome sequence of Cynopterus brachyotis ...

Conclusions
To date, the discussion on Cynopterus taxonomy in the Philippines is still fraught with conflicting evidence and, as a result, perspectives as well. Hence, this study was conducted so as to address part of these conflicts by attempting to provide further genetic evidence on Philippine Cynopterus through the sequencing of the mitogenome of Cynopterus brachyotis. This study was, indeed, able to generate an assembled Cynopterus brachyotis mitogenome from PCR-enriched bat tissue DNA, which was subsequently proven to be a high-quality and complete mitogenome sequence as determined through further quality assessment and annotation steps. The annotation steps revealed the relatively low similarity of the assembled mitogenome and individual genes with C. brachyotis and C. sphinx mitochondrial genomes and genes, respectively, suggesting that the mitogenome originated from a species that is genetically different to C. brachyotis (as originally identified); hence, phylogenetic analysis was performed for further confirmation. Phylogenetic analysis of whole mitogenome and CO1/cytb genes revealed and confirmed divergence of the Philippine Cynopterus cf. brachyotis from Asian Cynopterus (C. brachyotis and C. sphinx) with high branch support values. The results of this study altogether provide further evidence to confirm a previous proposition that Philippine Cynopterus cf. brachyotis is a unique species. Although part of the conflicting discussion on Philippine Cynopterus has been resolved by the results of this study, future studies to obtain nuclear marker/genome data, complementing the mitogenome data as well as morphological, behaviour and ecological data, are still needed so as to further establish the taxonomic resolution and the geographical extent of the Philippine Cynopterus cf. brachyotis.