Literature DB >> 24465949

Species radiation of carabid beetles (broscini: mecodema) in new zealand.

Julia Goldberg1, Michael Knapp2, Rowan M Emberson3, J Ian Townsend4, Steven A Trewick5.   

Abstract

New Zealand biodiversity has often been viewed as Gondwanan in origin and age, but it is increasingly apparent from molecular studies that diversification, and in many cases origination of lineages, postdate the break-up of Gondwanaland. Relatively few studies of New Zealand animal species radiations have as yet been reported, and here we consider the species-rich genus of carabid beetles, Mecodema. Constrained stratigraphic information (emergence of the Chatham Islands) and a substitution rate for Coleoptera were separately used to calibrate Bayesian relaxed molecular clock date estimates for diversification of Mecodema. The inferred timings indicate radiation of these beetles no earlier than the mid-Miocene with most divergences being younger, dating to the Plio-Pleistocene. A shallow age for the radiation along with a complex spatial distribution of these taxa involving many instances of sympatry implicates recent ecological speciation rather than a simplistic allopatric model. This emphasises the youthful and dynamic nature of New Zealand evolution that will be further elucidated with detailed ecological and population genetic analyses.

Entities:  

Mesh:

Substances:

Year:  2014        PMID: 24465949      PMCID: PMC3900486          DOI: 10.1371/journal.pone.0086185

Source DB:  PubMed          Journal:  PLoS One        ISSN: 1932-6203            Impact factor:   3.240


Introduction

Biologists have long been perplexed with New Zealand's biotic composition which cannot be classed as typically oceanic or continental [1], [2]. The landmass is continental in stratigraphic composition and origin, founded on a tectonic fragment separated from Gondwanaland ∼80 million years ago (Ma) - referred to as Zealandia [3], [4], [5]. Its enigmatic biota includes so-called relicts (tuatara Sphenodon), a supposed Gondwanan element (e.g. leiopelmatid frogs, Onychophora, kauri Agathis) and recent colonists [6], [7], [8], [9], [10]. This eclectic biotic composition has singled New Zealand out as being a particularly intriguing example of island biology [11]. An initial emphasis by phylogeographers on the origin of New Zealand lineages (stem groups) tended to undervalue the greater evolutionary significance of crown groups in the assembly of the New Zealand biota. Increasingly, molecular studies are revealing the story of diversification in a wide variety of New Zealand animals and plants (see for example [10], [12], [13], [14], [15], [16], [17], [18], [19]). There appears to be little congruence in phylogeographic patterns expressed by New Zealand taxa with no overwhelming pattern among lineage formation, landscape history and distribution of taxa (see reviews in [8], [9], [20]). In retrospect this may not be surprising, for the early expectation that species level phylogeography would reveal scaled-down patterns of the type seen in the northern hemisphere have rarely been realised (but see [21], [22], [23]). On the contrary, even within species, genetic diversity is often very high [24] and spatial patterns often complex, indicating that even on the relatively small scale of New Zealand population diversity relates to older events and/or involves persistence of larger populations than those commonly identified for European taxa [18], [20]. In many respects northern hemisphere biology has been dominated by Pleistocene climate extremes causing fluctuating population size, but New Zealand biology appears to have an older signature and identification of the respective roles of past geophysical events is in its infancy [but see for example tiger beetles- [19], weta- [24], cicadas- [25]). We know that after separation from Gondwanaland starting 80 Ma, the continent of Zealandia gradually submerged beneath the sea, and that modern New Zealand is primarily the product of tectonic activity initiated ∼25 Ma [3], [4], [5]. How much land persisted in the region prior to this phase is debated, but extreme reduction is certain and this period is likely to have caused a biological bottleneck at the very least [26]. Thus, whatever the source of lineages in New Zealand (oversea dispersal or Zealandian survivors), crown group diversity is unlikely to be older than Miocene. Tectonic activity since that time has resulted in substantial remodelling of the archipelago, including crustal uplift and orogenesis since Pliocene time (∼5 Ma; [27], [28]). Biotic turnover since early Miocene time is also evident from the fossil record and this might reflect a shift from tropical to temperate climate (e.g. [29], [30], [31]). There is naturally an increasing interest in how diversification is distributed through time and space. In this study we examine the phylogenetic relationships and timing of radiation of Mecodema (Blanchard, 1843) carabid beetles (tribe Broscini). This endemic genus of large, flightless beetles constitutes a prominent species radiation in New Zealand and presents a good opportunity to explore species level diversification. We utilise the fact that the genus is represented on the Chatham Islands, which are located approximately 850 km east of the South Island, New Zealand in the Pacific Ocean. Geological evidence for the formation of this archipelago within the last 4 Myr is compelling [3], [32], [33], [34] and corroborated by genetic data for many taxa (e.g. insects – [18], [35], plants – [36], [37], [38], parakeets – [39], pigeons – [40], [41], cicadas – [42], invertebrates and plants – [43]). In this study the earliest possible establishment of an island biota (4 Myr) on the Chathams [32] is used as a maximum possible calibration for estimating the timing of diversification in Mecodema. Furthermore a substitution rate for Coleoptera [44] is employed to further explore timing of lineage diversification of this beetle genus.

Materials and Methods

Sampling

The genus Mecodema (Blanchard, 1853) belongs to the tribe Broscini (Carabidae). Broscini has a worldwide distribution but has its main diversity in the southern hemisphere (subfamily Nothobroscinae) [45] and consists of at least 27 genera, comprising about 80 species (see http://www.landcareresearch.co.nz/research/biosystematics/invertebrates/carabid/carabidlist) [46], [47]. Six endemic genera of Broscini are recognized in New Zealand, but Mecodema is especially species-rich. Adult Mecodema beetles are relatively slow-moving, nocturnal, flightless (with fused elytra), generally active throughout the year, and usually scarce [48]. As with other Carabidae, adults and larvae of the New Zealand taxa are predatory. Mecodema is a diverse genus with species distributed throughout the New Zealand mainland from alpine to coastal habitats. In contrast, there is a single species (Mecodema alternans) on the Chatham Islands. The same species occurs in southeast New Zealand near Dunedin. Although M. alternans may be better treated as a species complex [47], no morphological characters have yet been described that distinguish Chatham Island populations from those in mainland New Zealand ([49] & I. Townsend pers. obs.). In total our sampling comprised 113 specimens, with 88 Mecodema representing 35 described species, and 4 undescribed species of the 66 recognized Mecodema species (after [45], [46] and http://www.landcareresearch.co.nz/research/biosystematics/invertebrates/carabid/carabidlist), see Table 1 for details and authorities. Putative outgroup New Zealand Broscini in our sample included Oregus (Putzeys, 1868), Diglymma (Sharp, 1886), Brullea antarctica (Laporte de Castelnau, 1867), Metaglymma (Bates, 1867), and Bountya insularis (Townsend, 1971); plus one representative of Broscini from Australia, Chylnus ater (Putzeys, 1868) (Table 2). As many of the species in Mecodema are scarce and therefore difficult to collect we also made use of material from museum collections, and supplemented outgroup sampling with available GenBank sequences (Oregus septentrionalis AF466847 & AF466848; Oregus crypticus AF54423; Oregus inaequalis AF466850 [50]; Calathus aztec GU254333 [51]; Broscosoma relictum AF012502; Promecoderus sp. AF012499; Creobius eydouxi AF012498 [52]). The resulting sample represents the geographic and ecological range of Mecodema in New Zealand (Table 2). Mecodema species show examples of allopatric, parapatric and sympatric distribution. Species but not always species groups are limited to particular areas and only few species (e.g. M. crenicolle) are documented as being present in both main islands of New Zealand [47].
Table 1

List of Broscini species used in this study with information regarding authority, distribution (S.I. = South Island New Zealand, N.I. = North Island New Zealand; Ch.Is. = Chatham Islands) and number of individuals employed.

SpeciesDistribution# used in studyAuthority
M. alternans S.I. & Ch.Is.12Laporte de Castelnau, 1867
M. alternans hudsoni The Snares1Broun, 1909
M. crenicolle N.I. & S.I.9Laporte de Castelnau, 1867
M. crenaticolle N.I.5Redtenbacher, 1868
M. curvidens N.I.1(Broun, 1915)
M. fulgidum (cf fulgidum)S.I.5Broun, 1881
M. howittii S.I.2Laporte de Castelnau, 1867
M. longicolle N.I.1Broun, 1923
M. lucidum S.I.2Laporte de Castelnau, 1867
M. occiputale N.I.3Broun, 1923
M. cf oconnori N.I.4Broun, 1912
M. oregoides S.I.2(Broun, 1894)
M. rugiceps S.I.1Sharp, 1886
M. sculpturatum S.I.4Blanchard, 1843
M. huttense (cf huttense)S.I.3Broun, 1915
M. simplex N.I.4Laporte de Castelnau, 1867
M. spinifer N.I.4Broun, 1880
M. strictum S.I.1Britton, 1949
M. sulcatum N.I. & S.I.1(Sharp, 1886)
M. validum N.I.1Broun, 1923
M. rectolineatum S.I.1Laporte de Castelnau, 1867
M. politanum S.I.1Broun, 1917
M. impressum S.I.1Laporte de Castelnau, 1867
M. constrictum S.I.3Broun, 1881
M. costellum lewisi S.I.1Broun, 1908
M. costellum obesum S.I.1Townsend, 1965
M. allani S.I.1Fairburn, 1945
M. laterale S.I.1Broun, 1917
M. minax S.I.2Britton, 1949
M. elongatum S.I.1Laporte de Castelnau, 1867
M. metallicum S.I.1Sharp, 1886
M. ducale S.I.1Sharp, 1886
M. morio S.I.1(Laporte de Castelnau, 1867)
M. infimate S.I.1Lewis, 1902
M. punctatum S.I.1(Laporte de Castelnau, 1867)
Meta. moniliferum S.I.2Bates, 1867
Meta. aberrans S.I.5Putzeys, 1868
Meta. tibiale S.I.1(Laporte de Castelnau, 1867)
Brullea antarctica N.I. & S.I.1Laporte de Castelnau, 1867
Bountya insularis Bounty Is.1Townsend, 1971
O. aereus S.I.6(White, 1846)
O. inaequalis S.I.2(Laporte de Castelnau, 1867)
O. crypticus S.I.1Pawson, 2003
O. septentrionalis S.I.2Pawson, 2003
D. clivinoides S.I.4(Laporte de Castelnau, 1867)
D. obtusum S.I.2(Broun, 1886)
D. seclusum S.I.1(Johns, 2007)
Table 2

List of investigated Broscini samples with sample numbers corresponding to numbers in Figures, including species name, genes sequenced per specimen and sampling locations.

Sample IDSpeciesGenesLocation
MB 01* M. alternans aChatham Is., South East Is.
MB 02* M. alternans bChatham Is., South East Is.
MB 70* M. alternans aChatham Is., Mangere Is.
MB 71* M. alternans bChatham Is., South East Is.
MB 86* M. alternans bChatham Is., Mangere Is.
MB 87* M. alternans bChatham Is., South East Is.
MB 86.1* M. alternans dChatham Is., Mangere Is
MB 88* M. alternans dChatham Is., South East Is.
MB 88.1* M. alternans dChatham Is., South East Is.
MB 190* M. alternans dS.I., Otago, Takahopa River Mouth
MB 14* M. alternans aS.I., Dunedin, Taieri Mouth
MB 16* M. alternans aS.I., Dunedin, Sandfly Bay
MB185* M. alternans hudsoni dThe Snares [LUNZ2804]
MB 61* M. crenaticolle dN.I., Taranaki, Lake Rotokare
MB 62* M. crenaticolle cN.I., Taranaki, Lake Rotokare
MB 72* M. crenaticolle dN.I., Tongariro, Ohakune
MB 12* M. crenicolle dS.I., Nelson, Shenandoah
MB 44* M. crenicolle dS.I., Abel Tasman, Awaroa
MB 82* M. crenicolle dS.I., Abel Tasman, Awaroa
MB 79* M. crenicolle bS.I., Marlborough Sounds, Pelorus Bridge
MB 103* M. crenicolle bS.I., Nelson Lakes, Wairau River
MB 117* M. crenicolle dS.I., Lake Rotoroa, Braeburn Walk
MB 163* M. crenaticolle cN.I., Taranaki, Kaiteke Ra., Lucy's Gully
MB 186* M. crenaticolle dN.I., Waikato, Skyline Cave
MB 54* M. morio dS.I., Catlins, Purakaunui Falls
MB 105* M. infimate cS.I., Codfish Island
MB 110* M. fulgidum aS.I., Clarence Valley, Mt. Percival
MB 91* M. cf fulgidum bS.I., Seaward Kaikoura Ra., Mt. Lyford
MB 123* M. fulgidum dS.I., Hamner Springs, Forest Park
MB 81* M. fulgidum dS.I., Cobb Valley, Lake Sylvester Tr.
MB 07* M. fulgidum. dS.I., Seaward Kaikoura Ra., Mt. Fyffe
MB 17* M. punctatum dS.I., Rock & Pillar Range
MB 98* M. howittii aS.I., Canterbury, Banks Peninsula
MB 99* M. howittii bS.I., Canterbury, Banks Peninsula
MB 50* M. spinifer bN.I., Hawkes Bay, Mohi Bush
MB 18* M. spinifer aN.I., Auckland, Waitakeres, Arataki
MB 10* M. spinifer dN.I., Auckland, Waitakeres Ra.
MB 50.1* M. spinifer dN.I., Hawkes Bay, Mohi Bush
MB 100a* M. costellum lewisi cS.I., Road nr. Mt. White Station
MB 101* M. costellum obesum dS.I., Nelson, Takaka Hill, Canaan
MB 195* M. allani dS.I., Nelson Lakes, Matakitaki V. [LUNZ2372]
MB 197* M. laterale dS. I., Fiordland, Routeburn Track [LUNZ2383]
MB 178 M. minax dS.I., Catlins, Mt. Pye Summit [P23]
MB 180* M. minax cS.I., Catlins, Mt. Pye Summit [P22]
MB 160* M. elongatum cS.I., Otago, Kinloch
MB 143* M. metallicum cS.I., Buller, Fox River
MB 176* M. rectolineatum dS.I., Remarkables Range, Wye Creek [P15]
MB 196* M. politanum dS.I., The Remarkables, Rastus Burn [LUNZ2640]
MB 128* M. impressum dS.I., Queenstown, Kinloch
MB 134* M. constrictum dS.I., Craigieburn Forest P., Education Ct.
MB 35* M. constrictum cS.I., Fog Peak, Porter's Pass
MB 27* M. constrictum dS.I., Canterbury, Craigieburn Rec. area
MB 121* M. ducale cS.I., Lewis Pass, Lake Daniels Walk
MB 20* M. cf oconnori dN.I., Levin, 30B The Avenue
MB 73* M. cf oconnori cN.I., Te Urewera, Ngamoko Trig Tr.
MB75* M. cf oconnori dN.I., Dannevirke, Norsewood Res.
MB 21* M. cf oconnori aN.I., Wellington, Levin, Ohou
MB 76* M. simplex dN.I., Tararua Forest Park, Putara
MB 77* M. simplex cN.I., Tararua Ra., Mt Holdsworth
MB 25* M. simplex bN.I., Manawatu, Pahiatua Track
MB 64* M. simplex aN.I., Manawatu, Palmerston North
MB 63* M. longicolle aN.I., Ruahine Ra., Pohangina Valley
MB 69* M. validum aN.I., Tongariro NP, Whakapapanui Track
MB 90* M. oregoides aS.I., Christchurch, Ahuriri Scenic Res.
MB 147.1* M. oregoides cS.I., Christchurch, Ahuriri Scenic Res.
MB 19* M. lucidum aS.I., Otago, Carrick Range
MB 111* M. lucidum dS.I., Pisa Range
MB 03* M. rugiceps aS.I., Fiordland, Lake Harris
MB 45* M. sculpturatum aS.I., Dunedin, Ross Reserve
MB 125* M. sculpturatum dS.I., Catlins Forest Park, River Walk
MB 04* M. sculpturatum dS.I., Dunedin, Leith Saddle
MB 06* M. sculpturatum dS.I., Dunedin, Mosgeil, Silver St.
MB 09* M. huttense cS.I., Canterbury, Peel Forest
MB 46* M. cf huttemse dS.I., Canterbury, Peel Forest
MB 108* M. cf huttense bS.I., Canterbury, Peel Forest
MB 96* M. strictum aS.I., Nelson, Takaka Hill, Canaan
MB 95* M. sulcatum aS.I., Kaikoura, North of Ohau Point
MB 66* M. curvidens aN.I., BOP, Rotorua
MB 68* M. occiputale aN.I., Mangatoi, Otanewainuku Forest
MB 65* M. occiputale dN.I., BOP, Ohope Scenic Res.
MB 67* M. occiputale dN.I., Mangatoi, Otanewainuku Forest,
MB 11* M. n.sp.aS.I., Central Otago, Old Man Range
MB 37* M. n.sp.bS.I., Arthurs Pass NP, Dome
MB 51.1* M. n.sp.aN.I., Hawkes Bay, Havelock North
MB 124* M. n.sp.cS.I., Lewis Pass, Lewis Tops
MB 30* M. crenicolle dS.I., Abel Tasman, Awaroa
MB 38* M. crenicolle dS.I., Takaka, The Grove
MB 49* M. crenicolle dS.I., Takaka, Rameka Track
MB 106 Meta. moniliferum aS.I., Canterbury, Christchurch, Quail Is.
MB 107 Meta. moniliferum bS.I., Canterbury, Christchurch, Quail Is.
MB 198 Meta. aberrans dS.I., Canterbury, Lake Tekapo [LUNZ3089]
MB 187* Meta. aberrans dS.I., Otago, Cromwell, Bendigo
MB 181 Meta. aberrans dS.I., Otago, Old Man Range, Omeo Gully [P16]
MB 182 Meta. aberrans dS.I., Otago, Old Man Range, Omeo Gully [P17]
MB 184 Meta. aberrans dS.I., Otago, Old Man Range, Omeo Gully [P10]
MB 194 Meta. tibiale dS.I., Central Otago, Upper Clutha [LUNZ3088]
MB 85 Brullea antarctica cN.I., Manawatu, Himatangi Beach
MB 199 Bountya insularis dBounty Is., Proclamation Is.
MB 192 Chylnus ater dAustralia
MB 13 O. aereus bS.I., Otago, Danseys Pass
MB 41 O. aereus aS.I., Dunedin, Morrison St.
MB 28 O. aereus dS.I., Dunedin, 46 Morrison St.
MB 29 O. aereus dS.I., Dunedin, Sandfly Bay
MB 47 O. aereus dS.I., Dunedin, Silver Peaks
MB53 O. aereus dS.I., L. Onslow, Lammarlaws
MB 5 O. inaequalis dS.I., Dunedin, Miller Rd.
MB 48 D. clivinoides bS.I., Seaward Kaikoura Ra., Tinline Va.
MB 31 D. clivinoides aS.I., NW Nelson, Heaphy Track
MB 161 D. obtusum dS.I., Fiordland, Kepler Track
MB 162 D. obtusum cS.I., Otago, Catlins Coast, Tautuku
MB 159 D. clivinoides cS.I., Otago, Kinloch
MB 8 D. clivinoides dS.I., Nelson, Cobb Valley
MB 175 D. seclusum dS.I., Fiordland, Spey River Valley

(M. = Mecodema; Meta = Metaglymma; O. = Oregus; D. = Diglymma; S.I. = South Island New Zealand; N.I. = North Island New Zealand; BOP = Bay of Plenty; genes analysed: a = COI, COII, 16S, 18S; b = COI, COII, 16S; c = COI, 18S; d = COI; LUNZ = Lincoln University Entomological Research Museum; P = S.M. Pawson collection);

COI sequences from samples of the Mecodema ingroup marked with * have previously been deposited in Genbank (#JN409817–JN409904) [18].

(M. = Mecodema; Meta = Metaglymma; O. = Oregus; D. = Diglymma; S.I. = South Island New Zealand; N.I. = North Island New Zealand; BOP = Bay of Plenty; genes analysed: a = COI, COII, 16S, 18S; b = COI, COII, 16S; c = COI, 18S; d = COI; LUNZ = Lincoln University Entomological Research Museum; P = S.M. Pawson collection); COI sequences from samples of the Mecodema ingroup marked with * have previously been deposited in Genbank (#JN409817–JN409904) [18]. Fresh specimens were obtained in accordance with Department of Conservation collection permits and preserved in 95% ethanol after hand collection in the field. These specimens are in interim storage with unique voucher numbers as part of the Phoenix Collection, Massey University, Palmerston North. Additional taxa were loaned from relevant collections (labelled accordingly in Table 2). Species identification relied on morphological characters [46].

DNA extraction, amplification and sequencing

DNA was extracted from a single leg of each collected specimen using a salting-out extraction protocol [53] or a CTAB and phenol/chloroform extraction [54]. DNA extractions from pinned museum specimens were undertaken in a dedicated ancient DNA laboratory, remote from modern DNA facilities, using a CTAB and phenol/chloroform extraction [54] and following stringent protocols for handling ancient DNA [55]. We employed three mitochondrial genes and one nuclear for determining the species level relationships in Mecodema. Partial cytochrome oxidase I (COI) fragment was amplified for 24 specimens in this study and the remaining 89 COI sequences were drawn from GenBank (JN JN409817–JN409904), see Table 2. Additionally, partial cytochrome oxidase II (COII) and 16S were amplified for a subset of 39 specimens. The COI gene region (789 bp) was amplified using primers C1-J-2195 and L2-N-3014, COII (589 bp) using primers TL2-J-3037 and C2-N-3661, and 16S (769 bp) using N1-J-12585 and LR-N-13398 [56]. For museum specimens, which were expected to yield fragmented, low concentration DNA, Mecodema-specific COI primers were designed using the program Oligo4 (Molecular Biology Insights, Inc., Cascade, CO) to generate a series of short (∼130–200 bp) overlapping fragments (Table 3). Additionally a 967 bp fragment of nuclear rRNA (18S) for a subset of 44 Broscini was obtained using primers 18S-S22 and 18S-A1984 [57].
Table 3

List of Mecodema specific oligonucleotide primers for mtDNA cytochrome oxidase I (COI), designed for this study.

Primer nameSequence
MB_COI38 3′ GCTGATGTAAAATATGCTCG
MB_COI39 5′ GAGCATATTTTACATCAGCAAC
MB_COI03 3′ CAATGAATAAATCCTCCAA
MB_COI05 5′ GGAGGATTTATTCATTGAT
MB_COI91 3′ GTWGATCCAATTGATGAAAC
MB_COI88 5′ GTARTTTCATCAATTGGATC
MB_COI12 3′ CTTAAATATGATCATGTRG
Polymerase chain reactions were performed in 10 µl volumes and the amplified products then checked on a 1% agarose gel and purified using SAP/EXO1 digest (USB Corporation). Purified PCR products were sequenced using standard protocols for the ABI Prism BigDye Terminator Ready Reaction Kit (Applied Biosystems, Mulgrave, Australia) and run on an ABI Prism 377 automated sequencer (Applied Biosystems). Sequence identity was confirmed by comparison with published data, checked for nucleotide ambiguities in Sequencher 4.2 (Gene Codes Corporation, Ann Arbor, MI, www.genecodes.com) and aligned using Se-Al v2.0a11 [58]. The sequences have been deposited with accession numbers KF913050–913193 at GenBank (16S: KF913050–913088; 18S: KF913089–913130; COII: KF913131–913169; COI: KF913170–913193).

Phylogenetic analysis

To test whether Diglymma and Oregus species are the natural sister group to Mecodema we first analysed data from two genes separately (COI and 18S) as it was not possible to gain sequences for outgroup taxa outside of New Zealand for all the employed species and sequence availability in GenBank within Broscini was also very poor. Although Oregus and Diglymma represent two of the New Zealand Nothobroscina genera considered closest to Mecodema [45], it was crucial to verify this relationship within Broscini as two other potential outgroup taxa, Metaglymma and Brullea, exist. MrBayes 3.1.2 [59] was used to implement Bayesian analysis with the datasets applying a GTR model with gamma-distributed rate variation across sites and a proportion of invariable sites. Analyses with MrBayes used four independent Markov Chain Monte Carlo (MCMC) runs for ten million generations with a burn-in of 10% and a tree sampling frequency of 1000. Results were checked for convergence. Resulting posterior probabilities on the nodes were recorded. To examine the species phylogeny of the Mecodema group in New Zealand we employed all four genes (three mitochondrial and one nuclear) with a subset of 50 taxa (44 ingroup and 6 outgroup samples). The outgroup sampling was chosen after consideration of the results from the prior outgroup analyses. All taxa with data missing for no more than one of four genes were included in the phylogenetic analysis (Table 2). Partition-homogeneity tests (PHT [60]) were implemented in PAUP*4.0b10 [61] with 500 replicates for the combination of the gene regions to detect significant heterogeneity among the data sets. MrBayes 3.1.2 [59] was then used to implement Bayesian analysis with the concatenated dataset, applying a GTR model with gamma-distributed rate variation across sites and a proportion of invariable sites. The same model was applied to the partitions with rates and nucleotide frequencies for each gene unlinked. Analyses with MrBayes used four independent Markov Chain Monte Carlo (MCMC) runs for two million generations with a burn-in of 25% and a tree sampling frequency of 1000. Resulting posterior probabilities on the nodes were recorded. The same data were subjected to Maximum Likelihood analysis with bootstrap resampling incorporating a GTR model with gamma-distributed rate variation. ML analysis used RaxML [62] implemented via the CIPRES portal [63]. The data were partitioned by gene (COI, COII, 16S, 18S) and bootstrap resampling was halted by RaxML

Molecular dating

As fossil remains of Mecodema that could provide information for calibrating a molecular clock have not been found, yet, we had to rely on geological information and a substitution rate calculated for Coleoptera COI [44] for calibration. In order to gauge the timing and extent of species radiation in Mecodema within New Zealand COI sequences were obtained for additional specimens in addition to previous analyses (Table 2). In some cases, this drew upon museum specimens to further assess the stability of our inferences about the distribution of diversity and timing of radiation in this beetle group. In total 113 unique COI sequences were used for molecular dating in Mecodema (Table 2). To obtain estimates for the maximum age of lineage formation within the genus Mecodema we used this dataset of COI with two different calibration strategies. First we employed a stratigraphic calibration to estimate divergence times using the split between Chatham Island M. alternans and its closest relative on mainland New Zealand. We assumed a normal distribution for the age around a calibration value of 4 Myr, derived from the maximum age for the Chatham Islands land surface [3], [33], assuming that colonization was most likely sooner after emergence of the islands than later. Alternatively, to capture the minimum likely diversification dates we also calibrated the COI dataset with the substitution rate estimated by Pons et al. [44] for Coleoptera COI. This included a normally distributed prior on the substitution rate of 0.08606 subst/site/myrs/l, and a 95% HPD interval from 0.0253–0.147 subst/site/myrs/l as an approximation to the posterior distribution provided by Pons et al. [44]. This rate obtained from analysis of numerous beetle taxa is amongst the highest estimated for any animal gene, and other rates obtained for particular beetle lineages are much slower (e.g. 0.0211 subst/site/myrs/l [64], 0.02 subst/site/myrs/l [65]). Even these rates are nearly twice the widely employed 0.0115 subst/site/myrs/l estimate of Brower [66]. Age estimation for both datasets and both calibration strategies were conducted under the assumption of a strict molecular clock as well as assuming a relaxed molecular clock with a lognormal distribution of rates along the phylogenies [67]. The fit of both priors was compared using Bayes Factors. For all datasets and calibration strategies the relaxed lognormal clock fitted the data decisively better than the strict clock [68]. We used the software BEAUTI 1.4.8 [69] and BEAST 1.7.5 [70] for molecular dating with the given calibrations. All analyses were conducted with a Birth-Death tree prior and a random starting tree under the GTR+I+Γ model of nucleotide substitution. The MCMC was run for 50–100 million generations, sampling every 5000th–10,000th step after a discarded burn-in of 5–10 million steps. Each analysis was run at least two times. The program Tracer 1.4 [71] was used to summarize posterior distributions of all parameters in question, to verify convergence of the MCMC and to estimate Effective Sample Sizes (ESS). If the effective sample size was less than 200, a third MCMC run was conducted for the respective analysis. After convergence of the MCMC was confirmed, the posterior distributions of all parameters were estimated from the combined posterior distributions of all runs conducted for each analysis. The program FigTree 1.4.0 [69] was used to visualize the reconstructed phylogenies.

Results

Three widely used mitochondrial gene regions were employed to gauge the scale of genetic diversity among the Mecodema specimens. Overall we observed relatively low genetic distances among species of Mecodema with a maximum ML-distance of 0.0179 in COII (COI: 0.0161, 16S: 0.053). Lower values for 16S compared to COI and COII reflect the comparatively low proportion of variable sites in this gene (16.9%). Separate analyses of COI and 18S DNA sequences from 41 specimens of Mecodema and 6 outgroup taxa (Table 2) resulted in similar topologies even though sequence variation in 18S was low. These analyses confirmed Oregus and Diglymma as the sister group to Mecodema, and revealed the placement of Metaglymma and Brullea within the Mecodema radiation (Figs. 1A & B).
Figure 1

Bayesian phylogenetic trees of New Zealand Broscini.

Bayesian trees of a subset of Broscini samples to confirm ingroup and outgroup relationship; A) COI phylogney (10 million generation, run four times, burn in 10%), three taxa outside New Zealand Broscini (Bountya insularis, Chylnus ater, Calathus aztec) were chosen as outgroup (not shown); B) 18S phylogeny (10 million generation, run four times, burn in 10%), three taxa outside New Zealand Broscini (Broscosoma relictum, Promecoderus sp., Creobius eydouxi) comprised the outgroup (not shown).

Bayesian phylogenetic trees of New Zealand Broscini.

Bayesian trees of a subset of Broscini samples to confirm ingroup and outgroup relationship; A) COI phylogney (10 million generation, run four times, burn in 10%), three taxa outside New Zealand Broscini (Bountya insularis, Chylnus ater, Calathus aztec) were chosen as outgroup (not shown); B) 18S phylogeny (10 million generation, run four times, burn in 10%), three taxa outside New Zealand Broscini (Broscosoma relictum, Promecoderus sp., Creobius eydouxi) comprised the outgroup (not shown). The alignment of data from four gene regions comprising 50 specimens sampled across New Zealand (including 6 outgroup specimens –Table 2, Fig. 2) was 3114 bp long in total. All three mitochondrial genes displayed the average insect A-T content of about 75%. The partition homogeneity test (PHT) revealed no significant heterogeneity of lineage partitioning among the data sets (p = 0.866), suggesting their concatenation was appropriate. The GTR+I+Γ model of nucleotide substitution was identified as the best fitting model by the hLRT and the AIC as implemented in Modeltest 3.5 [72].
Figure 2

Maximum likelihood phylogeny for Mecodema.

A) Spatial distribution of samples used in this analysis. Symbols correspond to those in Fig. 2B and code for different clades. B) Analysis of concatenated dataset including mitochondrial COI, COII and 16S plus nuclear 18S. Values at nodes indicate ML bootstrap support returned by analysis using RaxML. Specimen numbers at tips are given as in Table 2.

Maximum likelihood phylogeny for Mecodema.

A) Spatial distribution of samples used in this analysis. Symbols correspond to those in Fig. 2B and code for different clades. B) Analysis of concatenated dataset including mitochondrial COI, COII and 16S plus nuclear 18S. Values at nodes indicate ML bootstrap support returned by analysis using RaxML. Specimen numbers at tips are given as in Table 2. Bayesian and ML analysis of the concatenated dataset supported a single topology with the same groupings and branching order and well-supported nodes (Fig. 2). Mecodema was confirmed as paraphyletic with respect to Metaglymma and Brullea as these fall inside the Mecodema complex throughout all datasets and analyses in this study. We note that in all cases, Metaglymma and Brullea are placed within the M. curvidens/origoides group. This phylogenetic position contradicts the current taxonomic classification and needs to be addressed further in the future. Clades revealed in this analysis comprise species that are, in many cases, widely distributed across New Zealand, and therefore spatial structure is not evident (Fig. 1A). The complete alignment of the COI gene region comprising 113 specimens (including 20 outgroup specimens within Broscini and another one outside of this tribe) was 789 bp long and had an overall A-T bias of 73%. Diglymma and Oregus species that represent two of the New Zealand Nothobroscina genera taxonomically closest to Mecodema [45] grouped, as expected, outside of the Mecodema clade (Fig. 3). The other potential sister taxa to Mecodema within Broscini, Metaglymma and Brullea, were not supported as being sister to Mecodema, but were nested within the Mecodema radiation (Fig. 3). The relationship among Mecodema species was consistent with other analyses conducted in this study, with clades similar to species groupings previously proposed [47]. Despite a good level of resolution in the gene trees, there was little evidence for spatial correlation of clades with the current terrain of New Zealand. In keeping with our four gene analysis, it was clear that no groupings of species in specific North/South Island clades are apparent, nor are specific lineages correlated to landscapes (Figs. 2 & 3).
Figure 3

The timing of Mecodema diversification.

A) COI Bayesian tree generated in BEAST. Numbers on nodes show age estimates based on stratigraphic calibration of 4 Myr. Outgroup taxa were the same as used in the outgroup test (Fig. 1). Grey branches indicate lineages present in North Island New Zealand, and black branches indicate South Island lineages. Tip symbols correspond to clade and location identifiers in Fig. 2. Coloured symbols match symbols in Fig. 3B, highlighting disjunct lineages in southern North Island and northern South Island. Asterisks indicate age estimates between N.I. and S.I. lineages (* = 1.60 Myr, ** = 1.69 Myr, *** = 2.03 Myr). B) Reconstruction of the paleogeographic environment in lower North Island, New Zealand ca. 3 million years ago, green areas indicating likely land above sea level during this time (modified from [28]). Black outlines indicate present day New Zealand land area with coloured symbols corresponding to those in Fig. 3A, showing the present sampling locations of some sister taxa in North Island and South Island. (Scale bar = 100 km.)

The timing of Mecodema diversification.

A) COI Bayesian tree generated in BEAST. Numbers on nodes show age estimates based on stratigraphic calibration of 4 Myr. Outgroup taxa were the same as used in the outgroup test (Fig. 1). Grey branches indicate lineages present in North Island New Zealand, and black branches indicate South Island lineages. Tip symbols correspond to clade and location identifiers in Fig. 2. Coloured symbols match symbols in Fig. 3B, highlighting disjunct lineages in southern North Island and northern South Island. Asterisks indicate age estimates between N.I. and S.I. lineages (* = 1.60 Myr, ** = 1.69 Myr, *** = 2.03 Myr). B) Reconstruction of the paleogeographic environment in lower North Island, New Zealand ca. 3 million years ago, green areas indicating likely land above sea level during this time (modified from [28]). Black outlines indicate present day New Zealand land area with coloured symbols corresponding to those in Fig. 3A, showing the present sampling locations of some sister taxa in North Island and South Island. (Scale bar = 100 km.) Accommodating a relaxed molecular clock approach we calculated divergence times by first using BEAST with a stratigraphic calibration and set priors for the NODE of 4 Myr (normal distribution), followed by use of a substitution rate estimated for COI in Coleoptera [44]. The two approaches resulted in different age estimates within the Mecodema radiation. For the given dataset of 113 taxa the ROOT was estimated at a maximum age of 13.42 Myr with the stratigraphic calibration of 4 Myr. However, using the rate of molecular evolution proposed from independent data for COI in Coleoptera [44], we found that date estimates were more than an order of magnitude smaller with 0.83 Myr for the ROOT and 0.27 Myr for the NODE (Table 4). Time estimates for principal nodes of the stratigraphic calculation are also shown in Fig. 3A.
Table 4

BEAST time estimates based on stratigraphic and COI substitution rate of Coleoptera [44] calibration.

NodeEstimates Myr (95% HPD) with stratigraphic calibrationEstimates Myr (95% HPD) with rate calibration
A 4 0.27 (0.08–0.46)
B1.24 (0.5–2.09)0.07 (0.02–0.15)
C8.38 (5.29–12.06)0.52 (0.22–0.9)
D6.39 (3.97–9.41)0.4 (0.18–0.72)
E9.8 (5.86–13.94)0.6 (0.27–1.07)
F10.4 (5.86–15.93)0.64 (0.26–1.12)
G 13.42 (7.84–19.51) 0.83 (0.37–1.45)

Letters correspond to letters in Fig. 3A, bold letters denote NODE (A) and ROOT (G).

Letters correspond to letters in Fig. 3A, bold letters denote NODE (A) and ROOT (G). Most of the Mecodema radiation appeared to be geologically young, with estimates of many clade origins less than 5 Myr. There was little spatial correlation within the COI dataset even at a broad scale such as between the two islands. Instead we noted multiple instances of North Island/South Island (NI/SI) splits within the phylogram (Fig. 3A, grey/black branches with coloured symbols correlating to map in Fig. 3B), at least one in each of the taxonomic groups in Mecodema. All but one of these splits (between the simplex/validum group to the rest of the inner clade) were younger than 5 Myr. The results estimated based on the substitution rate of Coleoptera [44] were even younger and placed the origin of all analysed Mecodema clades within the last 1 Myr (Tab. 4).

Discussion

In this study we explored the pattern and depth of species diversity in the beetle genus Mecodema, which was sampled broadly across known species subgroups [47] and geographic range in mainland New Zealand, the Snares Islands and the Chatham Islands. Trees with well-supported structure were returned from the combined data set, although some shallow species-level relationships remain unresolved. Most clades are stable across analyses and are largely consistent with the taxonomic groupings proposed by [46] and [47]. Incomplete resolution of some species-level relationships in Mecodema reflects the low levels of DNA sequence divergence in some cases, for example, at the COI locus, a maximum of 3.6% among more than 8 species in one sub-clade is observed. This observation in itself suggests recent speciation, which cannot readily be dismissed as taxonomic over-splitting as the fine-scale taxonomic subdivision of Mecodema is based on sound morphological and ecological criteria. Interestingly Mecodema diversity does not comprise only allopatric populations, which might be expected under a model of recent spatial isolation; many species exist in sympatry with congenerics. Furthermore, Mecodema is commonly also found with other carabids (e.g. Megadromus (Motschulsky, 1866)), implying additional niche competition [48] resulting in fine scale ecological delimitation. Thus, high species diversity in Mecodema might well be the product of intense selection yielding adaptive radiation. A striking example is that of Metaglymma (Bates, 1867), which by virtue of its distinct morphology was classified in a separate genus, but is probably better treated as an ecologically specialized Mecodema. Metaglymma may, along with the coastal monotypic Brullea, be included in Mecodema following further morphological and genetic investigation. The close relationship between these three genera, compared to the two other mainland New Zealand outgroup genera Diglymma and Oregus, is consistent with the degree of morphological differences among them [45]. Additional sampling of taxa and populations is necessary to interpret phylogeographic and taxonomic patterns in detail, but the current level of sampling is sufficient for the purposes of gauging the timing of species radiation. Molecular dating with appropriate calibrations provides an empirical approach to estimate timing of past speciation events and phylogeography [73]. Determining the age of endemic biota has been largely dependent on the use of molecular clock calibrations, and the dating of speciation in New Zealand's plants and animals is especially problematic due to the generally poor fossil record for many lineages. Even where fossils are present, their use requires good time constraint and confidence in their ancestral status for respective extant taxa [73]. Substitution rates of mitochondrial genes differ greatly among genes and lineages and therefore the use of a general invertebrate divergence rates would be inappropriate for this study. Bayesian relaxed clock methods on the other hand allow rates to vary among lineages although accuracy might still be influenced greatly by the setting of priors. In this study the inferred mitochondrial divergence times based on the fast COI rate obtained for Coleoptera [44] are consistently substantially younger than would be expected for a highly diversified genus even in the relatively youthful landscape of New Zealand. Despite a perception that New Zealand is an ancient land mass with an ancient biota, diversity is increasingly shown to be the product of comparatively recent speciation congruent with Plio-Pleistocene tectonic activity since the mid-Miocene [7], [8], [74]. Speciation of plants (e.g., [75]) and animals including vertebrates (e.g., [76], [77]) often correlates with relatively recent but profound environmental changes and habitat diversification. We conclude that radiation of the Mecodema crown group is unlikely to have started before the mid to late Miocene, with most lineage formation most likely in the Pliocene and Pleistocene. This relatively shallow radiation is therefore consistent with the timing of radiations inferred for other New Zealand invertebrates (cicada [43], weta [78], [79], cockroach [80], isopod [81]), vertebrates (galaxiid fish [82], [83], geckos [84]) and plants (buttercup [85], Pachycladon, [38], [86]). On the basis of a conservative rate estimate for molecular evolution most diversification seems to have happened since the Pliocene (>5 Myr). The inference from this is that net diversification (comprising species origination and species extinction) within this genus has been strongly influenced by recent geophysical processes, such as mountain building in the Pliocene, habitat shifting and land formation in the Pleistocene. On a broad scale, and considering the relationships within Mecodema and calibrations used, at least eight exchanges between the two main islands of New Zealand can be inferred (Fig. 3). All of the splits between North Island and South Island sister taxa include species now distributed in the mid or lower North Island of New Zealand; areas that were mostly below sea level until recent geological time [28]. Considering the availability of land during this time frame, biotic exchange between the two main islands was most probably from south to north. South Island has been of similar size as it is today, at least since early Pliocene, whereas southern North Island formed only in the last 1 or 2 Myr ([28] and references therein). Prior to this, the proto-North Island was probably smaller than today and at least in the Pliocene was further from northern South Island than it is now. More complex scenarios involving multiple range shifts, speciation events and extinction are possible, but cannot be inferred from the present data. The splits between sister taxa currently occupying either North or South Island range in estimated age using stratigraphy from 6.39 Myr to 1.60 Myr (Fig. 3A). This suggests numerous exchanges between these regions since the Pliocene involving, most probably, colonisation after over-sea dispersal, rather than in situ radiation from a single ancestor in each island. The dominance of South Island lineages within Mecodema, Metaglymma and the Diglymma/Oregus sister clade are consistent with the notion that colonisation has for the most part been northwards. The current sampling is not complete in terms of living species or populations. Additional sampling is most likely to reduce the inferred age of formation of species lineages (by dividing existing branches), which could increase the number of inferred dispersal events. Alternatively, detailed population sampling could indicate taxonomic revision is required, and this might lead to slightly older inferred ages of named species. Nevertheless, a number of observations can be made: Geographically proximate species are rarely sister taxa (e.g. M. costellum, M. crenicolle, M. strictum in north east South Island). Some localised species might have arisen via allopatry on past islands, for example M. n. sp. in eastern central North Island (light green triangle in Fig. 3). Many populations and species have their range (sampled and known) in areas that were under water during late Pliocene and early Pleistocene (e.g. M. crenicolle/crenaticolle, M. oconnori, M. simplex). The existence of these indicates speciation since land became available and we note a mismatch between estimated lineage age using Chatham stratigraphy and estimates of land age. For example the split within North Island and South Island M. crenicolle/crenaticolle complex is estimated at 2.87 Myr. This predates evidence for land in southern North Island by more than 1 Myr, and predates by at least 1.5 Myr the most substantial period of land connection between North Island and South Island in the late Pleistocene. This type of mismatch between dates derived from genetics and those from geology might result from lineages evolving before the extant species they yielded arrived where they occur today as species ranges do change over time. In North Island, New Zealand, molecular studies showing this type of range expansion include tree weta [78] and stick insects [22]. However, evidence of recent range expansion is most easily found when populations of the taxon concerned still exist at source locations as well as colonised locations. In Mecodema, range shift after lineage origination would require lineage splitting, persistence of two lineages at a source location, expansion of at least one lineage and extinction of that lineage at the source. An alternative and simpler explanation is that the molecular rates inferred for these data by Chatham stratigraphic calibration are underestimated, which has thus yielded over estimation of lineage age. Such an inference is supported by the much younger times yielded by calibration with the Coleoptera COI rate [44]. It is not possible to determine using the present data when, after emergence of the Chatham Islands, the Mecodema alternans lineage arrived there. In all likelihood the age estimates from stratigraphic calibration can be considered maximal. Using this approach, we can infer a rate of molecular evolution for Mecodema COI of 0.0059 subst/site/myrs/l [87], which is substantially lower than comparable estimates for Coleoptera that range from 0.02 [65] to 0.086 subst/site/myrs/l [44]. Together with a lack of taxonomic distinction between mainland and Chatham populations this is consistent with a scenario of colonisation late in Chatham Island history. If so, then the entire Mecodema radiation is supported as even younger than we are currently able to demonstrate, with a larger proportion of speciation since the late Pliocene. The fact that this genus and other large flightless insects are present on the Chatham Islands [8], [18] and must have got there over sea supports the plausibility of trans-Tasman colonisation of the ancestral lineage to New Zealand. While the lineage (stem group) in New Zealand might date back to the separation of Zealandia from Gondwanaland having persisted despite marine transgression before the Miocene, it might equally have arrived in New Zealand in the Miocene as additional land became available. Regardless of the timing and mode of origination of Mecodema in New Zealand, the more important evolutionary feature of this beetle, with respect to the assembly of the native biota, is diversification of the crown group. Evidently, a few million years have been sufficient for the evolution of a complex ecosystem comprising not simply allopatric subunits but an array of sympatric species with distinct ecologies; for this to happen a long geological history was not required. While it has been predicted that intense phylogeographic structuring and speciation dating to the Plio-Pleistocene might be observed more frequently in naturally subdivided alpine conditions than in lowland forests [88], Mecodema speciation appears to provide an example where diversification has proceeded across space and into diverse habitats, from coast to above the tree line. Future work on the detailed ecology of these species will be instrumental in demonstrating the population genetic and ecological mechanisms of diversification (e.g. [89]). Increasingly, the fields of species phylogenetics and population phylogeography are merging as it becomes easier to generate appropriate DNA data, and the focus in taxonomy is shifted towards an evolutionary paradigm (e.g. [90]. Teasing apart the interaction of abiotic and genetic processes on population subdivision remains challenging but Mecodema is one taxon group that will provide helpful insight, and it is already evident that Mecodema is an impressive example of recent species radiation in the New Zealand fauna. In recent years, synthesis of phylogenetic, ecological and taxonomic evidence has indicated that the biology of New Zealand is primarily the story of recent adaptation and speciation [8], [9], [19]. Understanding properly how our view of “recent” geological time relates to the intergenerational population genetics of large invertebrate populations will provide the basis of exciting research.
  39 in total

1.  Complex pattern of coalescence and fast evolution of a mitochondrial rRNA pseudogene in a recent radiation of tiger beetles.

Authors:  Joan Pons; Alfried P Vogler
Journal:  Mol Biol Evol       Date:  2005-01-12       Impact factor: 16.240

2.  The resurrection of oceanic dispersal in historical biogeography.

Authors:  Alan de Queiroz
Journal:  Trends Ecol Evol       Date:  2004-11-25       Impact factor: 17.712

3.  RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models.

Authors:  Alexandros Stamatakis
Journal:  Bioinformatics       Date:  2006-08-23       Impact factor: 6.937

Review 4.  New Zealand phylogeography: evolution on a small continent.

Authors:  Graham P Wallis; Steven A Trewick
Journal:  Mol Ecol       Date:  2009-08-07       Impact factor: 6.185

5.  Mega-island or micro-continent? New Zealand and its fauna.

Authors:  C H Daugherty; G W Gibbs; R A Hitchmough
Journal:  Trends Ecol Evol       Date:  1993-12       Impact factor: 17.712

6.  The New Zealand biota: Historical background and new research.

Authors:  R A Cooper; P R Millener
Journal:  Trends Ecol Evol       Date:  1993-12       Impact factor: 17.712

7.  Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae).

Authors:  P Sunnucks; D F Hales
Journal:  Mol Biol Evol       Date:  1996-03       Impact factor: 16.240

8.  Phylogeographical pattern correlates with pliocene mountain building in the alpine scree weta (Orthoptera, anostostomatidae).

Authors:  S A Trewick; G P Wallis; M Morgan-Richards
Journal:  Mol Ecol       Date:  2000-06       Impact factor: 6.185

9.  Ancient DNA provides new insights into the evolutionary history of New Zealand's extinct giant eagle.

Authors:  Michael Bunce; Marta Szulkin; Heather R L Lerner; Ian Barnes; Beth Shapiro; Alan Cooper; Richard N Holdaway
Journal:  PLoS Biol       Date:  2005-01-04       Impact factor: 8.029

10.  Exploring Phylogeographic Congruence in a Continental Island System.

Authors:  Julia Goldberg; Steven A Trewick
Journal:  Insects       Date:  2011-08-03       Impact factor: 2.769

View more
  2 in total

1.  Extreme convergence in egg-laying strategy across insect orders.

Authors:  Julia Goldberg; Joachim Bresseel; Jerome Constant; Bruno Kneubühler; Fanny Leubner; Peter Michalik; Sven Bradler
Journal:  Sci Rep       Date:  2015-01-16       Impact factor: 4.379

2.  Diversification of Angraecum (Orchidaceae, Vandeae) in Madagascar: Revised Phylogeny Reveals Species Accumulation through Time Rather than Rapid Radiation.

Authors:  Herinandrianina N Andriananjamanantsoa; Shannon Engberg; Edward E Louis; Luc Brouillet
Journal:  PLoS One       Date:  2016-09-26       Impact factor: 3.240

  2 in total

北京卡尤迪生物科技股份有限公司 © 2022-2023.