An investigation into the taxonomy of Abavorana luctuosa (Peters, 1871) (Anura, Ranidae) and the resurrection of Rana decorata Mocquard, 1890 from Borneo

The taxonomic status of the ranid frog Abavorana luctuosa (Peters, 1871) was investigated using a combination of molecular and morphological data. The analyses revealed that A. luctuosa sensu lato is composed of two species in Borneo. One of these species agrees with the description of Rana decorata Mocquard, 1890 which is resurrected in the combination Abavorana decorata comb. nov. (Mocquard, 1890). Abavorana decorata is recovered as the sister lineage to the remainder of Abavorana and differs by a 16.0– 17.0 % uncorrected pairwise sequence divergence from its congeners A. nazgul and A. luctuosa, respectively. It is distinguishable morphologically from A. luctuosa and A. nazgul by its ventral pattern (bold, black and white reticulations on its venter along with bold banding on the underside of hind limbs vs. generally immaculate and spotted in the latter two species), and a prominent white streak beneath the eye and/or tympanum extending to the corner of the jaw. Abavorana decorata further differs from A. luctuosa by having a significantly wider head and snout, larger interorbital and tympanum diameters, longer femur in both sexes, and various combinations of other mensural characters. Both species are sympatric in Borneo and this discovery adds to a growing number of widespread Sundaic species shown to be species complexes with distinct forms in Borneo.


Introduction
The Mahogany frog (Abavorana luctuosa) is a colourful riparian species heard at night calling from tangled thickets along the edges of small streams and rivers as well as standing bodies of water such as ponds, pig wallows and pools in intermittent streams on the Thai-Malay Peninsula and Borneo. Abavorana luctuosa was described by Peters in 1871 from an unspecified location in Sarawak, East Malaysian Borneo (Capocaccia 1957). It was once thought to be a single wide-ranging species across Sundaland until a second species A. nazgul, was recently described from Gunung Jerai in Peninsular Malaysia . In their analyses, Quah et al. (2017) also recovered populations of A. luctuosa from Sabah as polyphyletic and only distantly related to one another.
Building upon these preliminary findings, we hypothesize that these populations are distinct evolutionary lineages and test this hypothesis with increased sample sizes from Borneo, Thailand, Peninsular Malaysia and detailed morphological, colour pattern, and genetic analyses based on the mitochondrial genes 16S, ND1 and flanking tRNAs (tRNA-Leu, tRNA-Ile, and tRNA-Gln). Given the results from the preponderance of data analyses herein, we rescind Boulenger's synonymization of Rana decorata Mocquard, 1890 with Abavorana luctuosa and resurrect the species in the combination Abavorana decorata comb. nov. (Mocquard, 1890).

Phylogenetic analysis
The new sequences were added to the molecular phylogenetic dataset of Quah et al. (2017) along with additional GenBank sequences from Shimada et al. (2011) (Table 1). Ingroup samples consisted of 30 individuals representing three nominal species of Abavorana (Table 1). Outgroups species used to root the tree were Pulchrana sundabarat Chan, Abraham, Grismer & Brown, and Meristogenys stenocephalus Shimada, Matsui, Yambun & Sudin, based in part on the relationships of Quah et al. (2017). New sequence data generated for this study are deposited on GenBank (Table 1).
Maximum likelihood (ML) and Bayesian Inference (BI) were used to generate phylogenetic trees. The data were partitioned into four partitions (ND1 by codon positions [1-3] and 16S+tRNAs) and best-fit models of evolution were determined in the IQ-TREE webserver version (Nguyen et al. 2015;Trifinopoulos et al. 2016) using the Bayesian information criterion (BIC ;Schwarz 1978) implemented in ModelFinder (Kalyaanamoorthy et al. 2017). The model K2P+G4 was the best-fit model of evolution for ND1 codon position 1, HKY+I for codon position 2, TN+I for codon position 3 and TIM2e+G4 for 16S+ tRNAs. The ML analysis was performed using the IQ-TREE webserver (Trifinopoulos et al. 2016) with 1000 bootstrap pseudoreplicates using the ultrafast bootstrap (UFBoot) option (Minh et al. 2013;Hoang et al. 2017). The Bayesian analysis was performed in BEAST v2.5.1 (Bouckaert et al. 2014) using bModelTest (Bouckaert & Drummond 2017) to numerically integrate over uncertainty in models of substitution, while estimating the phylogeny using Markov chain Monte Carlo (MCMC). The models selected by bModelTest for ND1 codon position 1 was 123141 (posterior support 35.22%), codon position 2, 121131 (posterior support 6.22%), codon position 3, 121131 (posterior support 29.94%) and 123454 (posterior support 29.01%) for the noncoding partition 16S+tRNAs. We ran the MCMC chain for 100 million generations while sampling every 10,000 generations. Stationarity was checked with Tracer v1.7 (Rambaut et al. 2018) and the effective sample sizes for all parameters were greater than 200. A maximum clade credibility tree was generat-ed by summarizing the posterior distribution of trees after discarding the first 25%, using mean node heights with TreeAnnotator v2.5.1 (Bouckaert et al. 2014). We considered Bayesian posterior probabilities (BPP) of 0.95 and above and ultrafast bootstrap support values (UFBoot) of 95 and above as strong nodal support (Huelsenbeck et al. 2001;Minh et al. 2013). Uncorrected pairwise sequence divergences (p-distance) were calculated in MEGA v6.06 using default settings which used pairwise deletion for the treatment of gaps and missing data in the dataset (Tamura et al. 2013) (Table 2).

Morphological analysis
Morphological data from 62 adult specimens (34 males; 28 females) of Abavorana luctuosa sensu lato from 21 populations throughout the Thai-Malay Peninsula and Borneo (Appendix 1), and eight specimens of A. nazgul (seven males; one female) from the type locality in Penin-sular Malaysia, were taken following the methodology of Quah et al. (2017) (Table 3). Notes on colour pattern were taken from digital images of living and euthanized specimens prior to preservation. Sex was determined by the presence of the humeral glands in adult male Abavorana as they lack vocal sacs and nuptial pads (Inger 1966). The following characters were measured by ESHQ with a Control Company digimatic caliper to the nearest 0.01 mm on the right side of the body for symmetrical characters: snout-vent length (SVL), from the tip of snout to the vent; head length (HL), from the posterior margin of mandible to the tip of snout; head width (HW), measured at the level of the jaw articulation points; snout length (SL), from the anterior margins of eyes to the tip of snout; snout width (SW), distance between the anterior margins of eyes; interorbital diameter (IOD), distance across the top of the head between the medial margins of the upper eyelids at their closest points; internarial distance (IND), distance between the medial margins of the nostrils; eye diameter (ED), distance between the anterior and poste-   rior margins of the eyeball; tympanum diameter (TD), the horizontal width of the tympanum at its widest point; brachium length (BL), distance from the axilla to most distal point of inflection on flexed elbow; forearm length (FAL), distance from outer margin of flexed elbow to the base of the inner metacarpal tubercle; manus length (ML), distance from the proximal edge of the outer metacarpal tubercle to the tip of the third finger; femur length (FL), distance from the vent to the outer margin of the flexed knee; crus length (CL), distance from the outer face of flexed knee to the outer margin of tarsal inflection; tarsal length (TL), the outer margin of flexed tarsus to the base of the inner metatarsal tubercle; pes length (PL), measured from the proximal edge of the inner metatarsal tubercle to the tip of fourth toe; humeral gland length (HG), the horizontal length of the humeral gland. Discrete co-  Table 3 continued.

A. luctuosa
Borneo lour pattern data, along with univariate and multivariate analyses of morphological data (see statistical analyses below) were then used to search for characters and morphospatial patterns bearing statistically significant differences that were consistent with the previous designations of the species-level hypotheses, thus providing independent diagnoses to complement the molecular analyses.

Statistical analyses
All analyses were performed using the platform R version 3.2.1 (R Core Team 2015). Owing to a lack of genetic material from the type locality of Rana decorata at Mount Kinabalu, Sabah, an a priori decision as to the membership to which species to assign specimens prior to additional downstream analyses was based on the colour pattern of their venters which has been shown to be a distinguishing characteristic for the genus ( Fig. 9; Quah et al. 2017). Specimens with bold black and white reticulations along flanks, and banding on the underside of their legs are identified as the decorata pattern morph, while those with generally immaculate venters are identified as the luctuosa pattern morph. To reduce allometric bias as a result of ontogeny, the characters of each species were scaled using the following equation: X adj = log(X)-β[log(SVL)log(SVL mean )], where X adj = adjusted value; X = measured value; β = unstandardized regression coefficient for each population; SVL = measured snout-vent length; SVL mean = overall average SVL of all samples (Thorpe, 1983;Lleonart et al. 2000). Each species was scaled separately and the data concatenated so as not to conflate interpopulational variation. The measurements for the sexes were analysed separately to eliminate bias caused by sexual dimorphism. All raw measurements are presented in Table 3. A MANOVA analysis was performed to test for significantly different (p < 0.05) multivariate structure in the mensural data and subsequently viewed in a principal component analysis (PCA). PCA and discriminant analysis of principal components (DAPC) were performed using the ADEGENET package in R (Jombart et al. 2010) to determine if the morphospatial position of Bornean Abavorana luctuosa populations with boldly marked venters and other populations of A. luctuosa with generally immaculate venters from Borneo and the Thai-Malay Peninsula coincided with the putative species boundaries delimited by the molecular phylogenetic analyses. All mensural data were incorporated in the analyses except HG which is absent in females. Principal components of the DAPC with eigenvalues accounting for 90-99% of the variation in the data were retained for the DAPC according to the criterion of Jombart et al. (2010). For this analysis, the first seven PCs out of 17 (males) and 16 (females) were retained which accounted for 94.9% (males) and 97.6% (females) of the variation.
Independent of the MANOVA, an analysis of variance (ANOVA) was conducted on mensural characters. Since the species and population sample sizes were unequal, a Levene's test for homogeneity of variances among the mensural characters was conducted prior to an analysis of variance (ANOVA). A Welch's ANOVA was conducted on characters with unequal variances (p ≥ 0.05) and a standard ANOVA was conducted on characters with equal variances (p ≤ 0.05) to test for the presence of statistically significant mean differences (p ≤ 0.05) between all combinations of species or population pairs in the data set. Characters with equal variances and bearing statistical differences were subjected to a TukeyHSD test to ascertain which population pairs differed significantly from each other for those characters (Table 4). Characters with unequal variances bearing statistical differences were subjected to a Games-Howell test to ascertain which population pairs differed significantly from each other for those characters (Table 4). Summary statistics were generated for the mensural data (Table 5).

Phylogeny and genetic divergence
The BI and ML analyses recovered trees with identical topologies. Specimen FMNH 231052 with bold black and white mottling on the ventral side of its flanks and banding on the underside of its legs from Danum Valley, Sabah ( Fig. 6A & B) was recovered as the strongly supported (BPP=1.0, Ufboot=89) sister lineage of the remaining Abavorana (Fig. 1). The remaining A. luctuosa sensu lato from the Thai-Malay Peninsula and Borneo bearing generally immaculate venters and underside of their legs (Figs 1, 4) formed a strongly supported (BPP=1.0, Ufboot=96) clade, and were recovered as the sister species to A. nazgul (BPP=1.0, Ufboot=99). FMNH 231052 differs by an uncorrected pairwise sequence divergence of 16.0-17.0 % from A. nazgul and other A. luctuosa populations, respectively (Table 2). With the inclusion of additional Bornean specimens in this study, two reciprocally monophyletic lineages (BPP=1.0, Ufboot=96) of A. luctuosa with generally immaculate venters and underside of hind limbs were recovered-one from Peninsular Malaysia and another from Borneo (Fig. 1). Nevertheless, these populations differ from each other by only a 2.0-3.0 % uncorrected pairwise sequence divergence (Table 2).

Systematics and nomenclature
The analyses demonstrate unequivocally that Abavorana luctuosa on Borneo are composed of two genetically and morphologically distinct forms. Abavorana luctuosa described by Peters (1871) on the basis of a single female specimen from Sarawak (Fig. 3) has only very faint speckling on its venter and underside of the limbs. Therefore, specimens from southern Thailand, Peninsular Malaysia, and Borneo that possess this phenotype shall be referred to as A. luctuosa (Figs 4, 7E-H, 8A-F, 9A-B). The second lineage is readily identifiable by the reticulated pattern on their venters and boldly banded underside  of their hind limbs, and is endemic to Borneo. Nine specimens examined from the states of Sabah (Mount Kinabalu, Kota Belud, Maliau Basin, Tambunan and Danum Valley) and Sarawak (Bario) bear this ventral colour patter and match the description of Rana decorata Mocquard, 1890 that was described from Mount Kinabalu, Sabah, East Malaysian Borneo. They resemble the colour pattern of the type series (MNHNP 1889.226-28) (Figs 5, 6, 7A-D, 9C-D). Rana decorata was synonymised with A. luctuosa by Boulenger (1891) without clear justification yet followed in subsequent works (e.g., Boulenger 1920;Van Kampen 1923;Quah et al. 2017 and references therein). Thus, given the preponderance of data, R. decorata Mocquard, 1890, in the combination Abavorana decorata comb. nov. (Mocquard, 1890) is resurrected for the populations bearing this phenotype.
Diagnosis. Body robust, medium-sized; head moderate; snout short, rounded, canthus rostralis smoothly rounded; interorbital space broader than the upper eyelid; tympanum distinct, not quite two-thirds the size of the eye with no pale colouration on the margins of the tympanum; no vocal sacs in males; vomerine teeth in two small oblique patches on a level with the posterior edge of the choanae; length of 1st finger greater than 2nd finger; disc width to finger width ratios of finger 3 and toe 4 is 1-1.5; dorsolateral fold indistinct or absent; the humeral gland in males is prominent, raised and centrally positioned on the ventral surface of the upper arm; a weak or absent rictal ridge; outer metatarsal tubercle weak or absent; subarticular tubercles moderate; skin of dorsum smooth or finely shagreened; throat, abdomen, and flanks smooth; posterior section of venter and back of the thigh rugose; dorsum reddish-orange to chocolate-brown, encircled by a white or cream coloured dorsolateral line that encircles the snout, canthus rostralis, outer edge of the upper eyelids, and dorsum along the dorsolateral fold to the vent; lower flanks dark-brown or black below the dorsolateral line grading into a paler venter; dorsal colouration of the limbs same as the flanks with whitish or light-grey speckles or stripes. Abavorana luctuosa can be easily differentiated from its congeners on the basis of its ventral colour pattern which is usually immaculate or with only very faint, sparse, light speckling (Figs 3, 4, 9A (Table 5).

Abavorana decorata
Diagnosis. In addition to its phylogenetic placement (Fig.  1), Rana decorata is reassigned as a member of the genus Abavorana based on the combination of having a robust, medium-sized body; no vocal sacs in males; length of 1st finger greater than 2nd finger; disc width to finger width ratios of finger 3 and toe 4 is 1-1.5; dorsolateral fold indistinct or absent; the colour of the dorsolateral line being white or yellow; the humeral gland in males is prominent, raised and centrally positioned on the ventral surface of the upper arm; a weak or absent rictal ridge; outer metatarsal tubercle weak or absent; skin of dorsum smooth or finely shagreened; throat, abdomen and flanks smooth; posterior section of venter and back of the thigh rugose; no pale colouration on the margins of the tympanum; flanks dark-brown or black below the dorsal fold grading into a pale venter (Inger 1966;Oliver et al. 2015;Quah et al. 2017). Dorsum reddish-orange to rust-brown, encircled by a white or cream coloured dorsolateral line that encircles the snout, canthus rostralis, outer edge of the upper eyelids, and dorsum along the dorsolateral fold to the vent; lower flanks dark-brown or black below the dorsolateral line grading into a paler venter; dorsal colouration of the limbs light-grey or brown with speckling and prominent dark-brown or black bands. Abavorana decorata can be differentiated from its congeners on the basis of its ventral colour pattern which is reticulated in black and white especially on the lower flanks, underside of the limbs (especially hind-limbs) are boldly banded in black and white, and a prominent white streak under the eye and/or tympanum to the corner of the jaw ( (Table 5). Fig 10). Endemic (Fig. 7C) and Bario (Fig. 6A-B); and West Kalimantan, Indonesia: Melawi Regency (Fig. 7D). The species is expected to be wider ranging on the island. Distribution. Abavorana nazgul is thus far confirmed from the upper elevations of Gunung Jerai, Kedah, Peninsular Malaysia . Specimens resembling this species have been photographed at Thale Ban National Park, Satun Province in southern Thailand at approximately 500-600m in elevation around wild boar wallows but are only tentatively identified as A. cf. nazgul pending positive evidence (Fig. 8H).

Discussion
Prior to this study and that of Quah et al. (2017), Inger (1966) and Grandison (1972) noted that Peninsular Malaysian Abavorana luctuosa were smaller than Bornean populations and lacked stripes on the undersides of their limbs but did not realise they were comparing two different species, A. luctuosa and A. decorata. Abavorana decorata is recovered as the sister species to the remainder of the genus and this study suggests that the ancestor of the genus possibly originated in Borneo before dispersing across Sundaland where one lineage gave rise to A. nazgul in the uplands of northwestern Peninsular Malaysia and possibly southern Thailand, and the second lineage gave rise to the widespread A. luctuosa. These findings lend credence to the out-of-Borneo hypothesis that has been postulated for a number of other herpetofaunal taxa and other biotic groups (de Bruyn et al. 2014;Grismer et al. 2016;Quah et al. 2020b). Alternatively, the present-day endemism of A. decorata on Borneo might be the result of extinction in other parts of Southeast Asia, the mainland in particular. Abavorana decorata may have had a far more extensive distribution across Sundaland previously. Nevertheless, both hypotheses can be explored in greater detail in the future with biogeographic range evolution analysis based on larger molecular datasets and broader taxon sampling for both Abavorana and other related genera across Sundaland. Borneo has long been a focal point of herpetological research in Southeast Asia (Das 2006bInger et al. 2017;Malkmus et al. 2002;Stuebing et al. 2014 and references within), and new discoveries continue to be made especially with the aid of molecular data and careful taxonomic reappraisals (for example, Davis et al. 2019;Grismer et al. 2018;Hertwig et al. 2012;Karin et al. 2018;Matsui et al. 2020;Mediyansyah et al. 2019;Pui et al. 2017;Quah et al. 2019aQuah et al. , 2020a. Our findings continue to lend support to the growing evidence that widely distributed species tend to be complexes and that their true diversity is not fully realised in the region (such as Arifin et al. 2018;Chan et al. 2020a,b;Grismer et al. 2018;Mediyansyah et al. 2019;Munir et al. 2019;Quah et al. 2019aQuah et al. , 2020a. Our data also reveal some significant differences in the measurements of Abavorana luctuosa populations from the Thai-Malay Peninsula and Borneo which, along with the genetic data and their allopatry, indicates that these populations are evolving on separate evolutionary trajectories, which makes them candidates for their recognition as distinct species. However, we continue to consider them as conspecific due to their low genetic divergence (2.0-3.0 %) based on ND1, 16S and tRNAs of mitochondrial DNA. We are adopting a conservative approach for now pending a more extensive genetic examination as other studies using genomic datasets have revealed past hybridization events between Bornean and Thai-Malay Peninsula populations of other anuran species complexes (Chan et al. 2020c).