Review of geological and biogeographic background
The north coast region
To a large extent, the northern part of New Guinea, including parts of the Central Range, are of Pacific affinity. This includes also the Bewani, Adelbert and Finisterre Mountains, parts of the Papuan Peninsula (see also below) as well as New Britain and parts of the Arfak Mountains, Yapen, Biak and Waigeo islands, among others (Fig. 1). Biologists have discussed the biogeographic significance of these elements of Pacific affinity, which has been referred to as the “Solomons Arc” (e.g., [28, 51]); summary in [33]. These authors interpreted this arc to represent an accreted island arc system that shows pronounced species level endemism, possibly derived from times when these terranes were islands adrift in the Pacific. This hypothesis implies an old origin of endemic lineages, as well as an initial stepping-stone dispersal from an Asian or a Pacific source area and then along an island arc with subsequent local speciation. This in turn requires that part of the island arc system was broadly exposed above sea-level for a considerable part of its history, e.g., since the mid-Miocene (15–11 Ma, [33] or earlier (> 30 Ma, [43]. We use this assumption as hypothesis H1 as suggested in the introduction.
A very poorly known area of Pacific affinity is the Gauttier Terrane, with the Foja Mountains as its prominent feature. It largely consists of basaltic lavas and breccia (and potentially ultramafic rocks) overlain by carbonates and volcanic ash [3, 77]. Few geological data exist from this region, but it is considered to be related to the Torricelli Terrane to the east, that is a terrane of Pacific affinity thought to be part of an island arc that accreted to New Guinea’s north coast after the early Miocene (< 23 Ma). Colonization of the Foja Mountains could thus be a comparably old biogeographic event, and this relates to our hypothesis H1A.
The ultramafics
Collision of other material of Pacific affinity with the northern edge of New Guinea likely occurred earlier and is represented by the various ultramafic/ophiolite rocks found mainly across central and western New Guinea [6, 51] (Fig. 1). These ophiolite and ultramafic rocks consist of oceanic lithosphere pushed southward and upward on top of the existing New Guinea landmass (i.e., due to obduction). They are considered to have been emplaced onto the northern margin of New Guinea at different times and locations between the late Cretaceous and the early Miocene (i.e., 100 Ma to ca. 20 Ma) [3, 14]. These seafloor rocks were subsequently ‘sandwiched’ between the fold and thrust belt of New Guinea’s Central Range, as well as the northern accreted terranes of Pacific affinity (“Solomons Arc” discussed above) during the late Miocene (11–5 Ma) to present day. The ophiolite sequences include the Irian Ophiolites (in Papua Province, Indonesia), the April Ultramafics (western Papua New Guinea, PNG), and the (east) Papuan Ophiolite (Papuan Peninsula, PNG) [2, 3] and these broadly correspond to the “Papuan Arc” terrane of [51]. We note that there are other areas of similar geology in the Cyclops Mountains [78] and the Mount Gamey area north of the Weyland Mountains [79], however, limited geological data in these regions means that little is known about their provenance and tectonic history.
Pioneering biogeographic work by [51, 58] suggested that the “Papuan Arc” was adrift in the Pacific and gave rise to genus level endemism. The ophiolites were suggested to have also played a role as stepping-stones for Asian fauna on their way to present day New Guinea. This hypothesis is not strongly supported by geological evidence, as: (1) ophiolites by definition are sections of the seafloor and the upper mantle that have been thrust on to continental crust—this means that ophiolites found on landmasses today, were covered by > 1000 m of seawater prior to their obduction, and could not harbor terrestrial life, and (2) the New Guinea ophiolites certainly stem from Pacific areas remote from an Asian source area. Note that the younger volcanic arc rocks of Pacific affinity would have likely been submerged volcanoes or small volcanic islands prior to their collision with New Guinea (see [41]). However, for our present phylogenetic analysis of one genus-level radiation, we suggested hypothesis H2, i.e. that some of the older clades should occur in, or nearby, areas of ophiolite and ultramafic rocks (e.g., Mesozoic ultramafic rocks, Fig. 2).
Central Range, Bird’s Head and Papuan Peninsula
The spine of New Guinea is the 1300 km long and up to 150 km wide central highland chain, otherwise known as the Central Range, consisting of the major geographic features, Maoke, Bismarck and Owen Stanley Range. It stretches from the easternmost end of the Bird’s Head area (ca. 135°E) to ca. 145°E, which is the westernmost end of the Owen Stanley range on the Papuan Peninsula. The latter is the Bird’s tail and eastward extension of the Central Range. More than 20 summits lie above 4000 m altitude, and summits above 3000 m are common. It includes a major fold-and-thrust belt in the Central Range which represents the deformed passive margin of the Australian continent. Prominent features include the Snow and Star Mountains as well as the Papuan Peninsula. Available geological data indicates that uplift propagated westward along the northern section of Papua New Guinea (i.e., within the ‘Solomons Arc’ region) at 8–5 Ma—this was associated with the strike-slip motion at the plate boundary at this time [18]. The Central Range was uplifted within the last 5 Ma, with uplift propagating from the north towards the south [4, 18]. There is some evidence that indicates this uplift may have propagated westward from Papua New Guinea, into Papua at ~ 2–3 Ma [5, 13], as well as eastward (e.g., [75]), including rapid rates of uplift after 3 Ma in the Papuan Peninsula (e.g., the Dayman Dome: [80] [81]).
There is no doubt that New Guinea also records evidence of earlier tectonic events that would have driven uplift in different parts of the island (e.g., a model proposed by [82] suggested that an underthrusting of the Australian continent beneath an Inner Melanesian arc resulted in an orogeny restricted to eastern New Guinea (forming the Papuan Peninsula Orogen, Fig. 1, ca. 30–35 Ma). This conceptual model is supported by low temperature thermochronology data that shows rocks of the Müller Range underwent a phase of early Eocene to Oligocene cooling and interpreted to provide potential evidence of the initiation of plate collision and associated uplift [83]. This event, together with a drop in global sea level during the late Oligocene-early Miocene (e.g., [84, 85]) resulted in large areas of New Guinea being exposed above sea-level (e.g., [7, 75, 86]). However, much of New Guinea was submerged again during the mid-late Miocene (ca. 5–15 Ma) due to subsidence outpacing a continued decline of global sea level. The time and extent to which land began to emerge during and after this time ultimately depends on the interpretation of the available geological data.
One piece of evidence that has been used to explain the emergence of land during the Miocene are the rocks classified as the Makats Formation, found in the Mamberamo region north of the Central Range. These sedimentary rocks were most likely deposited in relatively deep water, but contain detritus, including metamorphic rocks that must have been sourced from exposed land [86]. There is no firm evidence to indicate whether the source of the material was from the north (e.g., island arcs) or the south (e.g., uplifted ‘Australian’ crust) [86]. Despite this, [75] proposed that the metamorphic detritus within the Makats Formation were most likely sourced from a > 500 km distant landmass that first emerged in the westernmost part of the Central Range. The timing of this emergence was originally reported by [75] ca. 16–14 Ma corresponding to the maximum depositional age of the Makats Formation. We revised this reported age using the planktonic foraminifera that occur within the Makats Formation (i.e., first reported by [86] together with the currently accepted age ranges for the diagnostic planktonic foraminifera within the Makats Formation using the Mikrotax database (http://www.mikrotax.org/pforams—last accessed 7 April 2020) [87]. Based on this, the Makats Formation was most likely deposited between ca. 13.4 Ma and ca. 10.5 Ma (and most certainly no younger than ca. 6.8 Ma, based on the oldest ages reported for the overlying Mamberamo Formation). Smaller islands may also have been present from ca. 15 Ma. This is consistent with low-temperature thermochronology data and thermal history modelling (e.g., [4, 18, 83]. The Central Range region was considered to have grown laterally and in height, gaining maximum elevations of up to 2000 m by 8 Ma, and 4000 m by 6 Ma [75]. Towards the east (Papua New Guinea, PNG), the Central Range orogeny appears to be younger [82]. The Central Range continued to grow as the fold-and-thrust belt rapidly developed from 5 Ma to present. This period of time also saw the formation of mountains associated with crustal stretching and metamorphic core complex development (e.g., the Wandamen Peninsula in West Papua [13] and the Dayman Dome of the Papuan Peninsula, e.g., [80, 81]). New Guinea’s emergence from the sea was progressive as there were multiple phases of uplift and submergence in western New Guinea over the past 5 Ma [7]. So, while sections of the Central Range were being uplifted from 8 Ma, other regions such as the Papuan Peninsula, were not uplifted until the Early Pliocene, and the final stage of uplift across the island likely occurred between 1 Ma and present day [7].
Some studies suggest an older age for the emergence of the area spanning the Central Range of eastern PNG and the Papuan Peninsula including Herzog Mountains, e.g., referring to that region as “Woodlark Plate” [42], or, more focussed on the Papuan Peninsula including Herzog Mountains, as “East Papuan Composite Terrane, EPCT” [28, 33, 50]. The area interpreted as the EPCT is a concept that was proposed by [77]. It is said to consist of a series of crustal fragments that were torn off the northern margin of New Guinea in the Cretaceous or earlier. These rocks supposedly accreted with other crustal fragments from the Pacific between 52 and 23 Ma, all of which later collided with the northern margin of New Guinea during the Miocene [77]. The EPCT formation could also be linked to what [82] describe as the Oligocene Peninsular orogeny (ca. 35–30 Ma). The remnants of the EPCT today span the modern-day Papuan Peninsula including the Herzog Mountains. Based on this conceptual tectonic model, some biogeographers have assumed that the EPCT may have formed islands or a landmass north of New Guinea before 25 Ma, or at least representing one of the oldest terrestrial habitats in the proto Papuan archipelago (e.g., [23, 33, 42, 45, 48, 51]). In such conceptual models, the EPCT (and “Woodlark Plate” in [42] are proposed to harbor old biota and serve as source area for other parts of New Guinea.
However, the geological data on which the concept of the EPCT was formed simply indicates that a body of water developed between rocks of Australian Plate affinity along the northern margin of New Guinea [77]. It is unclear: (1) how wide this gap was; (2) whether all land connections were severed, and (3) whether this body of water was an ocean, or a rift valley that was subsequently inundated during higher sea levels. Available geological data does not indicate whether the EPCT was or was not exposed above sea-level before its final accretion to the northern margin of New Guinea during the Miocene. Considering these points, there is a high level of uncertainty as to whether the EPCT existed, and if it did, to what extent it may have been an emergent landmass within the Pacific Ocean prior to the middle to late Miocene. Rather than interpreting our data using a highly uncertain, conceptual tectonic model, we instead focus on whether there is any difference in which species are found within the different terrane affinities in the Papuan Peninsula, as well as the evidence that indicates the Papuan Peninsula was uplifted from the Miocene or later. To allow comparisons between biogeographic studies, in this paper, the rocks that others recognize as belonging to the EPCT correspond with the region marked by the East Papuan Ophiolite (Fig. 1). Similarly, the concept of the “Woodlark Plate” and its long-term biogeographic significance as discussed in [42] is highly problematic. This work used a map of the modern-day plate boundaries [88], where plate boundaries have largely been drawn on the basis of earthquake locations and GPS velocity data that suggest parts of the Earth’s crust are moving as a cohesive plate. Therefore, this does not outline entities for historical analyses that span millions of years.
Our hypothesis H3 relates to that, suggesting an early diversification on the present Central Range, possibly in an initial setting of a chain of islands, and subsequent colonization of surrounding areas such as the Bird’s Head and the Papuan Peninsula. This scenario was in part tested by [41], who suggested recent colonization of the Bird’s Head and Papuan Peninsula’s out of the Central Range, and [45] who found a colonization sequence from craton and later Papuan Peninsula (EPTC), to the Bird’s Head region and then areas with Pacific affinity.
We test this further under the following assumptions: (hypothesis H3A) the 1300 km long Central Range initially consisted of several islands, this implies localized radiations in different present-day highland blocks; (hypothesis H3B) there is a temporal sequence from west to the east; (hypothesis H3C) mountains of eastern PNG have existed as a separate island before 25 Ma and do therefore harbor fauna older than expected based on the above geological scenarios, and served as source area for other parts of New Guinea.
Assembly of the geological map
The geological terrane map that is shown in Fig. 2, and online at https://arcg.is/189zmz, is a simplification of numerous 1:250,000 scale geological maps of Indonesian New Guinea (currently Papua and West Papua Provinces, previously known as “Irian Jaya”) and Papua New Guinea. The Indonesian maps were developed by the Indonesian Geological Research and Development Centre (GRDC) (Pusat Penelitian dan Pengembangan Geologi (Indonesia)) and the Australian Bureau of Mineral Resources (BMR), now known as Geoscience Australia. The Papua New Guinea maps were developed by the Papua New Guinea Mineral Authority as well as the Australian Geological Survey Organisation (AGSO), now known as Geoscience Australia (for map sources and methods, see Additional file 1: Appendix 1).
Digital scans of the Indonesian New Guinea geological maps were orthorectified in ArcGIS v10.4 using the WGS84 datum. Each geological unit was mapped as a separate polygon and was assigned metadata according to the map rubric. Each polygon was classified according to one of seven geological terranes (Table 1). For Papua New Guinea, GIS maps were purchased from the PNG Mineral Authority (PNG_Geol250, 2002) and these were classified according to the same seven geological terranes as were used for Indonesian New Guinea. Some minor editing of the polygons was made to cut the digitized polygons to fit the current coastline. While every effort was made to quality check the data, readers should note that the digitization and reclassification involved individually modifying the attributes of > 20,000 polygons. Those polygons with a common terrane attribute and shared boundary were later merged to reduce the size of the datafile. Considering the size and complexity of the geological map, we expect that there will be minor errors. Readers should also note that the geological boundaries that we present have been digitized from hardcopy paper maps. There is some component of uncertainty associated with the original hardcopy maps. For instance, the regions of adjoining map sheets show considerable differences in the extent or continuity of particular rock types, particularly at the international border between Indonesia and PNG. Also, most of the maps were drawn before GPS was widely available, meaning there is an issue with the true location of the base maps that were used as well as the geologist’s ability to locate themselves on the map. Those who produced the PNG_Geol250 (2002) data estimated the geological boundaries in the digital dataset have an accuracy between 250 m (1 mm at 1:250,000) and 3.75 km (about 1.5 cm at 1:250,000 map scale), with the uncertainty being greatest in the highland regions and at the edge of adjoining map sheets. The Irian Jaya series maps likely have a similar level of spatial uncertainty. This does not account for uncertainty associated with distortion of the original paper maps before or during scanning. Having said this, readers should also note that the map presented in Fig. 2 is a much more detailed general terrane map of New Guinea compared to what is presented in earlier geological review papers (e.g., [2, 3]), and more importantly, the map is a much more accurate representation of the geology of New Guinea compared to the maps used in most existing biogeography papers.
Geospatial analysis
The latitude, longitude and altitudinal data obtained from a handheld GPS unit at each sample site were compiled in a database leading to 638 observations at 303 mapped localities. For instances where a reliable altitude could not be obtained in the field (e.g., due to dense tree cover), we assigned a value based on the sample site position and the corresponding cell from a digital elevation model (constructed from the ETOPO1 Bedrock Global Relief Model: [89]. Each database entry was classified with a number between one to five to represent the altitudinal range of each sampling site (< 500 m; 500–1000 m, 1000–1500 m, 1500–2000 m and > 2000 m), as well as a number between one and five to represent the geological terrane that each sampling site was encapsulated by (using the values recorded in Additional file 2: Appendix 2). The numerical value for the altitude and terrane was required to perform the Grouping Analysis tool (discussed below). Each database entry was double-checked as there were several instances where a sample site lay near a boundary of two terranes. In such instances, a decision was made to retain this classification, or to assign it to another terrane. This decision was most commonly employed when sample sites were located near the boundary between ‘Australian Plate affinity’ and ‘Pacific Plate affinity’ and were instead reassigned to the ‘Transition’ classification.
A one-to-many database join was used to assign the latitude, longitude, altitude and geological terrane to each identified species. Each database entry was also assigned a unique value within the database (Additional file 1: Appendix 2).
The Grouping Analysis tool, part of the ArcGIS Spatial Analysis toolbox extension was used to explore potential spatial relationships within the database. This technique utilizes unsupervised machine learning methods to determine natural groupings within a dataset. This technique is considered unsupervised because it does not require a set of pre-classified features to guide or train the algorithm that determines groupings within a spatial database. While our ultimate aim of employing this technique was to determine if there are spatial patterns within our database, we did not apply a spatial constraint to the algorithm, which means that the algorithm assumes there is no spatial correlation between any two sample locations. This means that if a spatial pattern is identified within the analysis, it is dependent on the altitude or underlying geology (or both). The tool requires the user to specify how many populations between 2 and 15 might exist within the dataset. After running the algorithm using different input parameters, we limited the final output to five populations. Readers should also note that it is possible that the numerical value assigned to represent the altitudinal range and geological terrane might also influence the results determined by this algorithm. Considering these limitations, we used the tool for data exploration purposes only, using it to test if any apparent relationship could be identified between species and altitude, species and underlying geology and species plus both altitude and geology.
Taxon sampling and molecular biology
We build upon the recent molecular framework of [41, 46], and added data from 41 additional species of Exocelina diving beetles compared to the most recent molecular treatment of the genus [46]. Most of the data sequenced for these new taxa was derived from older museum specimens. Complete genomic DNA was extracted with a Qiagen DNeasy Blood & Tissue kit (Hilden, Germany) using head and pronotum, or entire beetles. Using PCR protocols described in [41], we amplified and sequenced fragments of the following genes; mitochondrial cytochrome c oxidase I (cox1), cytochrome c oxidase II (cox2) and cytochrome b (cob), in addition to the nuclear histone 3 (H3), histone 4 (H4), 18S rRNA (18S), Carbomoylphosphate synthase (CAD) and Alpha-Spectrin (Asp). All new sequences will deposited in GenBank after publication of their formal names in [90, 91] so that new data can be uploaded with their final names. The matrix for this project has been deposited in Dryad: https://doi.org/10.5061/dryad.ns1rn8prj.
Alignment and phylogenetic inference
The resulting sequences were edited in Geneious R11 (Biomatters, USA) and checked for sequencing errors. Once assembled, the consensus sequences were aligned using MUSCLE [92] with existing datasets [41, 46] as well as a selection of outgroups chosen from the comprehensive phylogeny of [93] to facilitate the use of secondary calibrations in the BEAST divergence time analyses (see below). The resulting gene fragment alignments were checked for stop codons or indels and concatenated to produce a final matrix comprising 237 specimens (including 205 Exocelina specimens), and 4226 aligned nucleotides.
The concatenated matrix was analyzed in a maximum likelihood framework using IQ-TREE 1.6.6 [94]. The matrix was a priori subdivided into non-coding gene fragments and codon positions of coding gene fragments for a total of 22 partitions. The optimal partitioning scheme and corresponding models of nucleotide substitutions were searched simultaneously using ModelFinder [95] as implemented in IQTREE 1.6.6 among all available models and selected using the corrected Akaike Information Criterion (AICc) (see Additional file 1: Appendix 3 for the resulting partitioning scheme and models used in the IQ-TREE analyses). We performed 500 tree searches to avoid local optima. For each tree search, we performed branch support calculations with 1000 ultrafast bootstrap replicates (UFBoot, [96, 97] and 1000 SH-aLRT tests [98]. We used the hill-climbing nearest-neighbor interchange topology search strategy implemented in IQ-TREE to avoid severe model violations leading to biased ultrafast bootstrap estimations [97].
Divergence time estimation
There is no fossil known of the genus Exocelina, and the fossil record within Copelatinae is very scarce. Previous attempts at estimating absolute divergence times across Copelatinae have suggested different hypotheses pertaining to the temporal evolution of these beetles [41, 46, 65]. These discrepancies are largely linked to alternative calibration strategies. Recently, a new fossil-based dated phylogenetic framework for diving beetles has been developed based on the phylogeny of [66]. We rely on multiple secondary calibrations from this study to calibrate clocks and estimate divergence times. Specifically, we constrained the split between Dytiscinae and closely related subfamilies Copelatinae, Cybistrinae and Laccophilinae (i.e., the stem of Dytiscinae in [66] and the root of the tree in the current study) with a uniform distribution encompassing the 95% credibility interval recovered for this node in [66] (i.e., 122–141.4 Ma). We enforced the sister relationships between Laccophilinae and Cybistrinae in the BEAST analyses following [66]. We also constrained the crown of Cybistrinae (95% HPD = 37.4–81.5), crown of Laccophilinae (95% HPD = 44.7–98.3), crown of Copelatinae (95% HPD = 58–114.2), crown Cybistrinae + Laccophilinae (95% HPD = 71.6–120.6) and crown of Copelatinae + Cybistrinae + Laccophilinae (95% HPD = 88.8–135.5) with uniform prior distributions matching the 95% credibility intervals recovered in [66].
We used PartitionFinder 2 with the greedy algorithm, linked branch lengths and the set of models included in the program BEAST 1.10.4 [99], to select the optimal partitioning scheme and models of nucleotide substitution using the same 22 initial partitions as in the IQ-TREE analyses. The partitions and corresponding models of nucleotide substitution were selected with the Bayesian Information Criterion. We either assigned one clock to the mitochondrial partitions and another to the nuclear partitions (2 clocks in total), or a different clock to each partition recovered in PartitionFinder (12 clocks in total, see Additional file 1: Appendix 4 for the best partitioning scheme and models used in the BEAST analyses). We used relaxed molecular clocks with uncorrelated rates drawn from a lognormal distribution in BEAUti 1.10.4 [99]. The Tree Model was selected as Yule or birth death in different analyses. All other parameters were left to default. The analyses were conducted in BEAST 1.10.4 with 100 million generations, a parameter and tree sampling every 5,000 generations, and estimation of marginal likelihood using path-sampling and stepping-stone sampling with default parameters (chainLength = 1,000,000; pathSteps = 100; α = 0.3). The best scoring IQ-TREE topology out of 500 independent ML tree searches was constrained as a fixed input tree (with the unique modification being the enforced sister relationship between Laccophilinae and Cybistrinae) by manually editing the BEAUti.xml file. The four different analyses were compared based on their marginal likelihood estimates (MLE), and the one with the lowest MLE was used for further analyses.
Ancestral state reconstructions
We used the Bayesian Binary MCMC (BBM) method as implemented in RASP 4.2 [100] to estimate ancestral ranges across the New Guinean Exocelina radiation. To perform the reconstructions, we used the best BEAST maximum credibility clade tree based on MLE comparisons (Table 2). We coded each taxon with geology and geography (Additional file 1: Appendix 2). For altitude and geology, beetles were assigned to different categories using the geospatial model we compiled for this project. The altitude was coded using five discrete categories: (1) 0–500 m, (2) 501–1000 m, (3) 1001–1500 m, (4) 1501–2000 m, (5) above 2000 m. These intervals are subjective but provide the framework to visualize the three-dimensional distribution of the beetles. While we provide an outline of the elevational evolution of New Guinea in the introduction, note that it is currently not possible to suggest a reliable paleoaltimetric model for the island. Uplift was comparably fast with a rate of up to 10 mm/year−1 for some areas (e.g., the Wandamen Peninsula: [13] and up to 51 mm/year−1 for others (e.g., 17–51 mm/year−1 D’Entrecasteaux Island; [101]), at least from 3 Ma to the present. This might also apply for (parts of the) Central Range, which has seen rapid uplift in the past 5 Ma. But the rate of uplift might vary over time and region, and the supposed submergence of vast areas would add additional uncertainty. The geology was coded based on our geological terrane map of New Guinea. For the coding, the three principal terranes were in part subdivided as follows: Northern belt into: “accreted Pacific Plate affinity” and “ultramafic”; the transition belt into:”Post Collisional Volcanics” and “Transition”, and gondwanan material as “Uplifted Australian Plate affinity” (see introduction).
The geographic coding does not consider the geological history but looks at the existing island and its major regions (https://arcg.is/189zmz “New Guinea Regional Classification”). We here differentiate between six areas: the north coast mountain ranges including Wandamen (coded as A), the Bird’s Head including its satellite islands (B), the Bird’s Neck (Lengguru) (G), the Papuan Peninsula including the Herzog Mountains to its north east (C), as well as three blocks of the Central Range (from the Weyland/Paniai region in the west up to Baliem Valley at ca. 139° E (coded as D), from Baliem valley east to the Star Mountains including in Sandaun Province of PNG 142° E (coded as E), and finally the mountain block east of Sandaun to 145° E) (coded as F). This geographic delineation is simplified and meant to possibly reveal large scale biogeographic patterns and is not based on previously suggested New Guinea areas of endemism (e.g., [102]. It is also important to note that the delineation of the geographic areas is subjective to some degree; for example, in the west of area C, mainly Papuan Peninsula, we include the Herzog Mountains (that could alternatively be assigned to “A”, north coast ranges). The single species from the Bird’s Neck (Lengguru: Kaimana) region was coded in its own region (G). Alternatives would be assigning the Bird’s Neck to the Bird’s Head region (B) or the Central Range West (D) (Additional file 1: Appendix 11). Parts of the central range northern slopes, that are geologically also of Pacific affinity (Figs. 2, 4) were coded as geographically belonging to the Central Range.
All analyses were conducted in a Bayesian framework in RASP 4.2. using a Markov chain Monte Carlo method with 10 chains running for 1 million generation and with a sampling every 1,000 generations. We used the estimated F81 model for all runs. The rest of options were left to default in BBM.
Diversification rate dynamics estimation
We estimated diversification rate dynamics within Exocelina using the program BAMM 2.5.0 [103]. The analyses were performed with four reversible jump MCMC running for 1 million generations and sampled every 1000 generations. Parameter priors were estimated in R using the setBAMMpriors function (expectedNumberOfShifts = 1.0; lambdaInitPrior = 0.800; lambdaShiftPrior = 0.067; muInitPrior = 0.800). We used different priors (0.1, 1, 2 and 5) for the parameter controlling the compound Poisson process that determines the prior probability of a rate shift along branches of the chronogram. Our taxon sampling comprises 142 species of the New Guinean radiation (described and undescribed ones), yet we hypothesize that the extant species richness of the New Guinean radiation is closer to ca. 190 species. Therefore, the global sampling fraction was setup to 0.75 (New Guinean radiation only) in the different analyses. This is arguably more realistic than relying on the current described diversity (relying on the latter did not affect the results, data not shown). The BAMMoutput files were analyzed using BAMMtools 2.1.6 [104]. The posterior distribution of the BAMM analysis was used to estimate the best shift configuration and the 95% credible set of distinct diversification models.
We also tested the fit of various diversification dynamics scenarios to the entire New Guinean Exocelina diving beetle radiation. We relied upon constant‐rate, time‐dependent and diversity-dependent models of diversification as implemented in a maximum‐likelihood framework. The different models were fitted using the fit_bd function rom the R package RPANDA 1.8 [105] and the dd_ML function from the R-package DDD 4.3 [106], (see i.e., [107] for more details). Missing taxon sampling at the species level was also taken into account, using a global fraction of the expected species richness in the genus (i.e., 142/190 = ca. 0.75).
We tested the fit of the following models: (1) speciation rate constant through time with no extinction (BCST), (2) speciation and extinction rates constant through time (BCSTDCST), (3) speciation rate varying exponentially through time with no extinction (BtimeVarEXPO), (4) speciation rate varying linearly through time with no extinction (BtimeVarLIN), (5) speciation rate varying exponentially through time with constant extinction (BtimeVarDCSTEXPO), (6) speciation rate varying linearly through time with constant extinction (BtimeVarDCSTLIN), (7) extinction rate varying exponentially through time with constant speciation (BCSTDtimeVarEXPO), (8) extinction rate varying linearly through time with constant speciation (BCSTDtimeVarLIN), (9) speciation and extinction rates varying exponentially through time (BtimeVarDtimeVarEXPO), (10) speciation and extinction rates varying linearly through time (BtimeVarDtimeVarLIN), (11) speciation rate varying linearly with diversity without extinction (DDL), (12) speciation rate varying linearly with diversity with constant extinction (DDL + E), (13) speciation rate varying exponentially with diversity with constant extinction (DDX + E), (14) speciation and extinction rates varying linearly with diversity (DDL + EL).
In the time-dependent and diversity-dependent models, speciation and extinction rates (respectively λ and μ) could vary as a continuous function of time or diversity. This function was assumed to be either linear or exponential. The parameters α and β measure the sign and rapidity of time-variation for respectively speciation and extinction rates. Positive values of α or β can be interpreted as an indicator of speciation or extinction slowdown, while negative values indicate an acceleration of speciation or extinction. The parameter K measures the carrying-capacity in diversity-dependent models. These 14 models were compared with the AICc and ΔAIC to determine best-fit to the time-calibrated phylogeny.
