DNA extraction was conducted from whole specimens (excepting the large-bodied N. scirpi for which a leg was used) using non-destructive procedures and Omega Mag-Bind® Blood and Tissue DNA Kit (Omega Bio-tek) in the KingFisher robotic system (Thermo Fisher Scientific inc.). PCR amplification was done for the 5’ end COI gene (standard barcode region for Metazoa; Hebert et al. 2003) using degenerate Folmer barcode primers (FoldF: ‘TCNACNAAYCAYAARRAYATYGG; FoldR: ‘TANACYTCNGGRTGNCCRAARAAYCA’) (Folmer et al. 1994; Yu et al. 2012). PCR conditions were: 10 min at 95°C in 10 min, followed by 40 cycles of 30 s at 95°C, 30 s at 48°C, and 3 min at 72°C; 10 min at 72°C, and holding at 10°C. PCR products were cleaned using exonuclease and rapid alkaline phosphatase, and were Sanger-sequenced with ABI technology in Macrogen, Spain. This procedure was applied to 16 Raymondionymini specimens and Notaris scirpi (Table 1). We used Geneious Prime 2023 to visualize and edit chromatograms and to generate an HKY distance matrix and UPGMA tree to explore similarity between specimens. The agreement between the obtained barcode sequences and the morphological identification of specimens corroborated the presence of 8 species of Raymondionymini belonging to 6 genera.

One representative per species was selected for mitogenome sequencing and assembly following the mitochondrial metagenomics approach (Andújar et al. 2015; Crampton-Platt et al. 2015; Dettai et al. 2012), where complete mitochondrial genomes are assembled from shotgun sequencing of specimen DNA mixtures. The dsDNA concentration of raw DNA extracted from each specimen was measured using a Qubit 2.0 Fluorometer (Life Technologies Corp., Carlsbad, CA), and five TruSeq nano DNA libraries were constructed respectively from five equimolar DNA pools, each of these including one or several target specimens of this study plus a number of other Coleoptera from distant lineages (Table 1). Sequencing was performed with the Illumina MiSeq platform (Illumina Inc., San Diego, CA) (2 × 300 bp; 700–900 bp insert size), aiming a coverage of 1% of the MiSeq run per specimen. Illumina output was processed with Trimmomatic 0.30 (Bolger et al. 2014) for Illumina adapter removal. Reads were subsequently filtered using Blast 2.2.27 (Altschul et al. 1990) against a reference database including 2344 mitochondrial genomes longer than 5,000 bp retrieved from NCBI nucleotide database (accessed 5^th November 2020). Retrieved mitochondrial reads were then assembled using RAY 2.3.1 (Boisvert et al. 2012) (–K 61; –minimum-seed-length 100 –minimum-contig-length 1000), SPADES 3.14 (Prjibelski et al. 2020) (–k 21,33,55,77,99,127), and IDBA 1.1.3 (Peng et al. 2012) (–maxk 300 to –mink 50). The resulting contigs from the three assemblers were re-assembled in Geneious using the de novo assembly function and showed wide overlap, minimizing potential problems associated with the formation of chimeric mitogenome sequences. Obtained mitogenomes were annotated using gene predictions with MITOS (Bernt et al. 2013) with additional manual editing performed in Geneious. All mitogenomes were structured following the putatively ancestral gene order for the Coleoptera. Finally, mitogenomes were unambiguously assigned a taxonomic identity by comparison against the COI barcode sequences obtained from the same specimens with PCR-Sanger sequencing (see above). For two of these specimens (Table 1), mitogenome sequencing was performed twice in two independent libraries, allowing to corroborate reliability of obtained sequences.

A first and preliminary analysis was designed to construct a guide tree including our nine newly generated mitogenomes plus available Curculionoidea mitogenomes within the NCBI nucleotide database. This guide tree was to serve two purposes. Firstly, we wanted to preliminarily replicate the basal weevil branching events reported in earlier studies and summarized in Shin et al. (2018). Secondly, the obtained tree will guide the selection of the near and distant outgroups for a subsequent and statistically more exhaustive analysis using a lesser number of terminals. For these purposes, all 423 mitogenomes longer than 5000 bp classified as Curculionoidea available from GenBank on 14^th July 2021 were downloaded. Of these, 59 mitogenomes classified as unspecified Curculionoidea were excluded. From the remaining 382 mitogenomes, 94 lacked gene annotations and where de novo annotated in MITOS as indicated above. The final dataset included 391 gene annotated mitogenomes. Single gene datasets for each of the 13 protein coding genes (PCGs) and the two ribosomal genes were extracted using Geneious and individually trimmed and aligned using the FFT-G-INS-i algorithm of MAFFT (Katoh et al. 2002). Individual gene alignments were concatenated, yielding (i) a dataset of 15 genes and 13,491 bp (Preliminary Dataset 1; PD1); (ii) a dataset with exclusively the 13 PCGs and a length of 10,842 bp (PD2); and (iii) a dataset with amino acids sequences obtained from the 13 PCGs (invertebrate mitochondrial code) with a length of 3,856 AAs (PD3). These three datasets were used for Bayesian inference with PhyloBayes (Lartillot and Philippe 2004) running 2 chains under a GTR-CAT model for a minimum of 5,000, 6,000 and 3,500 generations respectively for PD1, PD2, and PD3. A consensus tree was obtained for each dataset combining trees from both chains after discarding the first 2000 generations as a burn-in fraction. In this and all subsequent analyses we used FigTree (Rambaut 2012) to visualize the obtained topologies.

We designed our restricted phylogenetic analysis based on congruence between trees obtained for the preliminary dataset (391 terminals) and the well-resolved weevil topology of Shin et al. (2018). The Raymondionymini ingroup finally included eight newly generated mitogenomes plus that of Coiffaitiella (MK692586, Andújar et al. 2019) (Table 1). The nearest outgroup was formed by four mitogenomes of Brachycerinae available from GenBank, plus the newly generated mitogenome of Notaris scirpi. The distant true weevil outgroup consisted of the representatives of the following clades/subfamilies (Shin et al. 2018): Dryophthorinae (5 mitogenomes), Platypodinae (3), Bagoinae (1) and “higher” weevils consisting of the CEGH clade (9) and CCCMS clade (19). To adequately place Curculionidae within the Curculionoidea phylogenetic framework, we added representatives of Brentidae (3 mitogenomes; the supposed sister group of Curculionidae), Attelabidae (3), Anthribidae (3), and Nemonychidae (1). We rooted all topologies on the clade formed by Anthribidae plus Nemonychidae consistently with the current hypotheses on the first split within weevils (Oberprieler et al. 2007). In total, 61 mitogenomes were analysed, nine of them newly generated. The 13 protein coding genes (PCGs) and the two ribosomal were individually extracted, aligned, and trimmed as before in Geneious, with the exception of ribosomal genes that were aligned using the online version of MAFFT with the Q-INS-i algorithm (Katoh and Toh 2008).

Individual gene alignments from the reduced dataset were concatenated, yielding (i) a dataset of 15 genes and 12,603 bp (Reduced Dataset 1; RD1); (ii) a dataset with exclusively the 13 PCGs and a length of 10,191 bp (RD2); and (iii) a dataset with amino acids sequences obtained from the 13 PCGs (invertebrate mitochondrial code) with a length of 3,396 AAs (RD3). These three datasets were used for maximum likelihood (ML) and Bayesian phylogenetic analyses. ML trees were obtained using RAxML v.8 (Stamatakis 2014) and IQTree (Nguyen et al. 2015), in both cases using gene partitions (DNA and AA datasets RD1, RD2 and RD3) and gene and codon partitions (DNA datasets RD1 and RD2). RAxML analyses were conducted on the CIPRES Science Gateway (Miller et al. 2010), applying an independent GTRGAMMA (DNA datasets) or PROTCATGTR model (AA dataset) to each data partition. The best scoring ML tree was selected among 1,000 searches on the original alignment with different randomized parsimony starting trees. Support values were obtained with 1,000 bootstrap replicates (Felsenstein 1985). IQTree analyses were run on the IQ-TREE web server at http://iqtree.cibiv.univie.ac.at (Trifinopoulos et al. 2016) using the best fitting substitution model for each gene partition as estimated with ModelFinder (Kalyaanamoorthy et al. 2017). Nodal support was obtained by 1,000 ultrafast bootstrap (UFBoot) replicates (Minh et al. 2013). For each dataset and partitions scheme, IQTree analyses were repeated twice. PhyloBayes (Lartillot and Philippe 2004, Lartillot et al. 2013) analyses were done on the CIPRES Science Gateway, running 2 independent chains under a GTR-CAT model. For each dataset (RD1, RD2 and RD3) analyses were duplicated, allowing to run on CIPRES for 48 and 72 hours respectively (and using between 64 and 96 cores). A consensus tree was obtained for each dataset combining trees from both chains after discarding the first 500 generations as a burn-in fraction. All together we conducted 21 phylogenetic analyses on the reduced datasets, as summarised in Table 2.

Our results are remarkably consistent with the gradually emerging phylogenetic framework of weevil (see Introduction). Our topologies display the composition and arrangement of the main weevil clades (Figs 4, S2–S25), as well as their statistical support (Table 2), nearly identical to those obtained for the analyses of hundreds of protein coding nuclear genes by Shin et al. (2018) and (Li et al. 2023). These consistent results include the monophyly of Curculionidae and its sister relationship with Brentidae, Brachycerinae as a grade of early splitting lineages within Curculionidae, the sister relationships of Dryophthorinae and Platypodinae, and a well supported group including Bagoinae plus the well supported clades CEGH and CCCMS as found by Shin et al. (2018). Also consistent between both studies is the lower support for the internal relationships within CEGH and CCCMS clades and a topology suggesting the polyphyletic nature of several large subfamilies within Curculionidae, such as Molytinae and Curculioninae. The herein recovered non-monophyly of Scolytinae, even though highly unlikely given morphological and biological data (summarized in Johnson et al 2017), has also been detected by Shin et al. (2018) and Mugu et al. (2018). These two groups of scolytines correspond to the basal split within this subfamily, between Scolytini and the rest (Pistone et al. 2018) and currently we consider this result as an artefact of our analysis.

Four of the five non-raymondionymine Brachycerinae genera here studied (i.e., Ocladius, Brachycerus, Lissorhoptrus and Echinocnemus) were among the seven Brachycerinae analyzed by Shin et al. (2018) and (Li et al. 2023). Ocladius and Brachycerus are consistently found as early splitting lineages in the case of Shin et al. (2018). Each of these lineages include an additional genus here not sampled (Ocladius + Schizomicrus and Brachycerus + Synthocus), whereas in the case of Li et al. (2023)Schizomicrus appears as an additional early splitting lineage. The additional Brachycerinae early splitting lineage consistently found corresponds to the supported clustering of Lissorhoptrus and Echinocnemus, with disagreement regarding their position. In 17 of our 21 analyses we found Lissorhoptrus + Echinocnemus as a sister clade to the rest of non-raymondionymine Curculionidae, whereas results of Shin et al. (2018) supported the sister relationship of a clade formed by Lissorhoptrus + Echinocnemus + Tanysphyrus with the clade grouping Platypodinae + Dryophthorinae. Further analyses with additional Brachycerinae representatives will be required to improve our knowledge on the early evolution of the polyphyletic Brachycerinae, where special attention should be placed to avoid long branch attractions effects among multiple early splitting lineages. In this sense, it will be especially relevant incorporating additional molecular markers from the nuclear genome to our dataset, in order to discard phylogenetic wrong conclusions derived from the use of single marker dataset and the potential misleading effect of compositional heterogeneity and long branch attractions (Sheffield et al. 2009; Song et al. 2010). Despite of these potential limitations, complete mitochondrial genomes have been shown to be highly informative and robust for relatively deep phylogenetic nodes as those targeted in the current study (e.g., Cameron et al. 2007; Liu et al. 2018; Talavera and Vila 2011; Timmermans et al. 2015, and including CurculionidaeGillett et al. 2014), with long branch attraction problems minimised and high consistency with nuclear markers reported when protein sequences are used and/or the site heterogeneous mixture model (CAT; Lartillot and Philippe 2004; Lartillot et al. 2013) are applied, as performed in this study (Liu et al. 2018; Timmermans et al. 2015).

Within the lineage formed by the Mediterranean raymondionymines, we have found that the genus Ubychia forms the sister clade to the rest of the subfamily in all our analyses. This result is consistent with (i) Ubychia being the easternmost representative of the subfamily and the only genus inhabiting the Caucasus and Elburs mountains (Hlaváč and Nakládal 2018), and (ii) with its notably distinct morphology, by having elytra smooth, rather than striate and deeply punctured (Fig. 3).

Mediterranean raymondionymine weevils, here represented by seven of the nine known genera, form a strongly supported lineage sister to the remaining Curculionidae, in agreement with the phylogenetic position reported by Andújar et al. (2019) for the genus Coiffaitiella. This statistically significant basal dichotomy within the family Curculionidae, together with the documented polyphyletic nature of the subfamily Brachycerinae, points to the need of a reclassification from their current taxonomic rank as a tribe within the subfamily Brachycerinae. Two alternatives are allowable: to consider raymondionymines either as a subfamily of the family Curculionidae, or elevate them to the status of a full beetle family as previously suggested by Alonso-Zarazaga and Lyal (1999). We advocate for the former option based on both morphological and physiological traits. Curculionidae is definable by at least five morphological synapomorphies (Oberprieler et al. 2007), all five are observable in raymondinymines: (i) adult geniculate antennae, (ii) compact adult antennal club, (iii) 3–4 dorsal folds in the larval abdominal segments, (iv) a prothoracic position of the thoracic spiracle in larvae, and (v) the frontal sutures of the larval head blocked by a frontoepicranial bracon. Despite the limited knowledge of raymondinymine larvae, these three characters are distinguished on the larva of Raymondionymus perrisi (Grenier), the only member of the subfamily with documented immature stages (Remillet 1968a). Although limited, available knowledge on feeding habits also support the inclusion of Raymondionyminae within Curculionidae. As presently defined, true weevils (Curculionidae) is a monstrous clade of some 51,000 species sister to the family Brentidae containing about 4,000 species (Oberprieler et al. 2007). Unlike other weevils feeding predominantly on fungi and conifers, the Brentidae plus Curculionidae clade is cladistically defined by an evolutionary novel colonizing of angiosperms and their use as the food source (Oberprieler et al. 2007). Among 55,000+ species of this clade, only a few feed on non-angiosperms. All such cases are likely subsequent evolutionary novelties. The only record of Raymondionyminae feeding habits (Remillet 1968b) indicates that at least Raymondionymus perrisi feeds on roots of angiosperm trees and shrubs, which is consistent with the herein proposed placement of this species inside Curculionidae.

The absence within our analyses of representatives from the six non-Mediterranean genera (19 species) prevents any conclusion about the global monophyly and the limits of the group. The relationship between the non-Mediterranean genera and the Mediterranean lineage has been previously questioned (Grebennikov 2010; Grebennikov and Anderson 2021), as the non-Mediterranean members were added to the group mostly based on the easily converging and often misleading similarity of small, eyeless, wingless, and depigmented deep soil dwellers. In this sense, using a morphological dataset, Grebennikov (2010) suggested that the genus Alaocybites, distributed in California and Russian Far East, is an unlikely member of the tribe. Similarly, Grebennikov and Anderson (2021) using a three DNA-marker dataset found no support for the relationship between an American unnamed genus standing close to Bordoniola and a well supported clade formed by the Mediterranean genera Alaocyba and Raymondiellus. European representatives were not included in the study by Shin et al. (2018), where the only raymondionymine weevil analysed corresponded to the genus Schizomicrus from California, which is not represented in our study. Still, the sister relationship between Schizomicrus and Ocladius in Shin et al. (2018), and the lack of such relationship between the European raymondionymines and Ocladius in our study, may be interpreted as an indirect evidence for the non-monophyly of the here resurrected subfamily Raymondionyminae.

Consequently, the phylogenetic position of all non-Mediterranean genera previously placed within Raymondionymini, including Alaocybites (two species endemic to California and two species endemic to the Russian Far East), Bordoniola (seven species in Ecuador and Venezuela), Gilbertiola (two species in California), Homosomus (three species in Madagascar), Neoubychia (monotypic in Mexico), and Schizomicrus (monotypic in California), can not be established according to current phylogenetic or morphological evidence and are left as incertae sedis within Raymondionyminaestat. rev., requiring further phylogenetic corroboration. Given the results of the present study, we hypothesize that additional currently undetected species-poor early offshoots of the true weevil radiation might await taxonomic recognition as subfamilies of Curculionidae. Most likely these lineages are hidden in what we consider the evolutionary twilight zone of true weevils. The latter is formed by the remaining members of the herein circumscribed subfamily Brachycerinae and the six non-Mediterranean genera of the subfamily Raymondionyminae. Bringing these intuitively classified organisms under a phylogenetic spotlight will significantly improve our knowledge on the early evolution and allow to fine-tune the systematics of the charismatic and megadiverse clade of true weevils.