Genomic signatures accompanying the dietary shift to phytophagy in polyphagan beetles

A representative sampling of the two major coleopteran suborders

Reliable estimation of gene gain and loss events requires a robust evolutionary framework, i.e., a phylogeny that includes the species studied, as well as the characterization of gene families across complete gene sets from these same species. To study adaptation to phytophagy, we sampled from both Adephaga (mostly predaceous) and Polyphaga (with diverse trophic habits, including a very large number of phytophagous species). A balanced sampling of each suborder was achieved comprising 12 transcriptomes and 6 genomes, with Benchmarking Universal Single-Copy Ortholog (BUSCO) completeness estimates [31, 32] ranging from 71.9 to 97% (Fig. 1, Table 1). The molecular species phylogeny estimated using protein sequences of 405 BUSCO genes found to be complete in all species and in the strepsipteran outgroup, Stylops melittae, was used to build the time-calibrated species phylogeny (Fig. 1). Subsequent analyses of gene gain and loss rates and of signatures of adaptive gene family expansion employed this ultrametric species tree, which importantly shows no significant difference in node ages within each suborder. Protein-coding sequence predictions ranged from 9844 to 24,671 genes per beetle species. These sequences matched 14,908 Arthropoda orthologous groups (OGs) containing at least 1 species of Coleoptera in the OrthoDB v8 catalog [40]. This represented a minimum of 6742 and a maximum of 11,149 OGs for Carabus frigidus and Leptinotarsa decemlineata, respectively. OGs containing genes from only 1 of the 2 sampled suborders were excluded, resulting in a total of 9720 OGs for the analysis that have evolutionary histories traceable to the last common ancestor of beetles. Functional annotations of the sequences within these OGs were used to identify and assign several of them to enzyme families relevant to the tested hypothesis. These candidate OGs comprised 4 UGTs, 22 P450s, 19 CEs, 6 GSTs, 4 SERs, 7 CYSs, 28 ABCs, and 1 GH, for a total of 91 candidate OGs from 8 families of genes (i.e., functional categories) (Table 2).

Species	Short form	Suborder	Family	Type of assembly	Accession	Bioproject	Source	BUSCO score (1658 genes)	Predicted proteins	OrthoDB groups (OG)
Cicindela hybrida	CHYBR	Adephaga	Carabidae	Transcriptome	GDMH01000000	PRJNA286505	1KITE, this study	C: 90.4 [S: 84.9%, D: 5.5%], F: 3.4%, M: 6.2%	13,916	8111
Calosoma frigidum	CFRIG	Adephaga	Carabidae	Transcriptome	GDLF01000000	PRJNA286499	1KITE, this study	C: 78.3% [S: 73.5%, D: 4.8%], F: 9.8%, M: 11.9%	9844	6742
Elaphrus aureus	EAURE	Adephaga	Carabidae	Transcriptome	GDPI01000000	PRJNA286520	1KITE, this study	C: 92.4% [S: 85.9%, D: 6.5%], F:2.2%, M: 5.4%	12,808	8020
Noterus clavicornis	NCLAV	Adephaga	Noteridae	Transcriptome	GDNA01000000	PRJNA286561	1KITE, Vasilikopoulos et al. [34]	C: 90.9% [S: 84.1%, D: 6.8%], F: 4.4%, M: 4.7%	12,981	7918
Haliplus fluviatilis	HFLUV	Adephaga	Haliplidae	Transcriptome	GDMW01000000	PRJNA286525	1KITE, Vasilikopoulos et al. [34]	C: 91.8% [S: 76.3%, D: 15.5%], F: 4.3%, M: 3.9%	19,408	8528
Cybister lateralimarginalis	CLATE	Adephaga	Dytiscidae	Transcriptome	GDLH01000000	PRJNA286512	1KITE, Vasilikopoulos et al. [34]	C: 88.2% [S: 83.1%, D: 5.1%], F: 3.7%, M:8.1%	13,916	7256
Sinaspidytes wrasei	SWRAS	Adephaga	Aspidytidae	Transcriptome	GDNH01000000	PRJNA286492	1KITE, Vasilikopoulos et al. [34]	C: 87.0% [S: 76.7%, D: 10.3%], F: 4.0%, M: 9.0%	13,392	7721
Dineutus sp.	DINEU	Adephaga	Gyrinidae	Transcriptome	GDNB01000000	PRJNA286516	1KITE, Vasilikopoulos et al. [34]	C: 71.9% [S: 51.6%, D: 20.3%], F: 13.6%, M: 14.5%	14,644	7089
Gyrinus marinus	GMARI	Adephaga	Gyrinidae	Transcriptome	GAUY01000000	PRJNA219564	1KITE, Misof et al. [35]	C: 81.5% [S: 79.0%, D: 2.5%], F: 8.1%, M: 10.4%	13,867	7663
Aleochara curtula	ACURT	Polyphaga	Staphylinidae	Transcriptome	GATW01000000	PRJNA219522	1KITE, Misof et al. [35]	C: 91.4% [S: 87.5%, D: 3.9%], F: 4.0%, M: 4.6%	20,280	8513
Anoplophora glabripennis	AGLAB	Polyphaga	Cerambycidae	Genome	GCF_000390285	PRJNA167479	I5k, McKenna et al. [21]	C: 96.9% [S: 95.8%, D: 1.1%], F: 2.7%, M: 0.4%	22,035	10,959
Agrilus planipennis	APLAN	Polyphaga	Buprestidae	Genome	GCF_000699045	PRJNA230921	I5k, unpublished	C: 92.5% [S: 91.2%, D: 1.3%], F: 4.5%, M: 3.0%	15,497	9089
Dendroctonus ponderosae	DPOND	Polyphaga	Curculionidae	Genome	GCF_000355655	PRJNA162621	Keeling et al. [36]	C: 91.2% [S: 86.0%, D: 5.2%], F: 4.1%, M: 4.7%	13,457	8518
Leptinotarsa decemlineata	LDECE	Polyphaga	Chrysomelidae	Genome	GCF_000500325	PRJNA171749	I5k, Schoville et al. [37]	C: 88.9% [S: 87.5%, D: 1.4%], F: 9.9%, M: 1.2%	24,671	11,149
Laparocerus tessellatus	LTESS	Polyphaga	Curculionidae	Transcriptome	10.5281/zenodo.1336288	N/A	This study, Seppey et al. [38]	C: 93.8% [S: 91.9%, D: 1.9%], F: 1.4%, M: 4.8%	18,448	8616
Meloe violaceus	MVIOL	Polyphaga	Meloidae	Transcriptome	GATA01000000	PRNJA219578	1KITE, Misof et al. [35]	C: 90.3% [S: 85.6%, D: 4.7%], F: 5.9%, M: 3.8%	14,295	8480
Onthophagus taurus	OTAUR	Polyphaga	Scarabaeidae	Genome	GCF_000648695	PRJNA167478	I5k, unpublished	C: 96.2% [S: 93.9%, D: 2.3%], F: 2.5%, M: 1.3%	17,483	9315
Tribolium castaneum	TCAST	Polyphaga	Tenebrionidae	Genome	GCF_000002335	PRJNA12540	Richards et al. [39]	C: 97.0% [S: 96.5%, D: 0.5%], F: 1.6%, M:1.4%	16,645	9429
Stylops melittae	SMELI	Outgroup	Stylopidae	Transcriptome	GAZM02000000	PRNJA219603	1KITE, Misof et al. [35]	C: 76.5% [S: 55.0%, D: 21.5%], F: 7.1%, M: 16.4%	13,026	6104

Gene family category	InterProScan (Pfam or InterPro identifiers) or Gene Ontology	UnifRef keyword	Number of OGs
UDP-glycosyltransferases (UGTs)	PF00201	name:“cluster UDP glucuronosyltransferase” OR name:“cluster UDP glycosyltransferase”	4
Cytochrome P450 oxidases (P450s)	PF00067	name:“cluster Cytochrome P450”	22
Carboxylesterases (CEs)	PF02230, PF00135	name:“cluster carboxylesterase” OR name:“carboxylic ester hydrolase”	19
Glutathione S-transferases (GSTs)	PF00043, PF02798	name:“cluster Glutathione S-transferase”	6
Serine proteases (SERs)	PF00450, PF12146, PF05577, GO:0008236	name:“cluster Serine protease” OR name:“cluster Serine peptidase”	4
Cysteine proteases (CYSs)	PF00112	name:“cluster cysteine protease” OR name:“cluster cystein protease” OR name:“cluster Papain”	7
ABC transporters (ABCs)	PPF00005, PF00664	name:“cluster ABC”	28
Glycoside hydrolases (GHs)	IPR000334, IPR000743, IPR001360, IPR001547	name:“cluster Glycoside hydrolase”	1
Total			91

Polyphaga exhibit more frequent gains across a larger set of OGs

Analysis of per-species gene counts of the complete set of 9720 OGs was performed with the Computational Analysis of gene Family Evolution (CAFE v3) [41] tool. CAFE analyses changes in gene family sizes using a stochastic birth and death process to model gene gain and loss across a species phylogeny. It estimates gain and loss rates from the extant gene count data, taking errors into account to allow for accurate inferences even with incomplete datasets, and identifies gene families with significantly accelerated rates of gain and loss. The mode considering distinct gene gain (λ = 0.0019 gain/gene/million years) and gene loss (μ = 0.0018 loss/gene/million years) was preferred over a single value for λ and μ, having a significantly greater maximum likelihood score (see the “Methods” section). The λ (gain) and μ (loss) values predicted when CAFE was run on each suborder separately were λ = 0.0020 and μ = 0.0027 for Adephaga, versus λ = 0.0023 and μ = 0.0021 for Polyphaga, showing a tendency for Adephaga to lose genes and for Polyphaga to gain genes. Among the 9720 OGs were 21 with reported expansions originating at the Adephaga root and 126 at the Polyphaga root (see Fig. 1 to locate the nodes). Conversely, 240 OGs showed gene losses for Adephaga and 354 for Polyphaga. Two expansions and 21 losses affected the candidate OGs for Adephaga, and 9 expansions and 6 losses for Polyphaga. Other polyphagan nodes leading to phytophagous-rich clades (i.e., Chrysomeloidea and Curculionidae) also exhibited more candidate OGs expanding than contracting (Fig. 1). All counts of gene gains and losses per node are presented in Additional file 1: Figures S1 and S2. Additionally, CAFE assigned individual OG p values of < 0.01 to a subset of 910 (9.3%) OGs, which, according to De Bie et al. 2006 [42], indicates gene families likely to have experienced accelerated rates of gain and loss. These are interesting to investigate further as they may represent large OGs of potentially unequal size between the suborders. Among these were 26 of the 91 candidate OGs (28.6%), a significantly larger proportion (two-sample test for equality of proportions, chi-square test, p value < 0.0001) compared with just 9.3% of non-candidate OGs.

Signatures of adaptive expansion are more prevalent in Polyphaga

All 910 OGs with significant variations in their gene content were tested for signatures of adaptive expansion in each suborder, by comparing Brownian motion (BM, neutral) to Ornstein-Uhlenbeck (OU, selective pressure) evolutionary models [43]. As mentioned previously, these included 26 OGs that belong to one of the functional categories listed in Table 2 (“candidate” OGs). The models consider per-species gene count as a trait that can evolve towards a value, which may or may not differ between the two suborders and may or may not be guided by selective pressure; we call this the “optimum” value in models integrating selection. The BM models assume no selection where differences between the suborders result from stochastic processes whose rates are estimated. The OU models assume that reaching an optimal gene family size in each suborder is driven by selective pressure. No adaptive LSE is represented by the null hypotheses of BM models with a single rate for the whole tree (BM1) or with different rates for each suborder (BMS) or the OU model with selection towards the same optimum for both suborders (OU1). Adaptive LSE is represented by the alternative hypotheses of OU models with selection towards two optima having the same variance (OUM) or with selection towards two optima having two variances (OUMV). In total, 21 OGs displayed a higher optimum for Adephaga (0.2% of the initial 9720 OGs) and 88 for Polyphaga (0.9%). Eight of these 88 OGs with higher optima in Polyphaga (Table 3; Additional file 1: Table S1; and gene trees in Fig. 2 and Additional file 1: Figures S3-S9) are candidate OGs belonging to one of the candidate gene families of Table 2, while none of the 21 OGs with higher optima in Adephaga belong to any of the candidate gene families. The proportion of OGs with expansions and higher optima in the background (all “candidate” and remaining “control” OGs) was significantly larger for Polyphaga compared to Adephaga (two-sample test for equality of proportions, chi-squared, 88/9720 vs. 21/9720, p value < 1e−09), indicating that Polyphaga have experienced globally more LSE under selection on protein-coding genes. Additionally, a test for enrichment (see the “Methods” section) of OGs with LSE under selection from the candidate families (Table 2) compared to the background was significant for Polyphaga (8/91 vs. 88/9720, p value < 1e−09). The same test applied individually on each candidate gene family within the candidate dataset demonstrated that categories enriched for LSE under selection in Polyphaga were GSTs (3/6 positive tests, FDR-corrected p value < 1e−09) and CEs (3/19, FDR-corrected p value < 0.005), as shown in detail in Table 4. Furthermore, functional genomics data from the polyphagan Asian longhorned beetle [21] supported a biological role in plant feeding activities for the candidate OGs that tested positive for adaptive expansions, which were enriched with genes upregulated in larvae fed on sugar maple trees vs. a nutrient-rich artificial diet (36/114 vs. 1391/12,461, p value < 1e−10).

Category	ODB8 ID	BM1 AICc H0.1	BMS AICc H0.2	OU1 AICc H0.3	OUM AICc H1.1	OUMV AICc H1.2	Mean Adephaga	Mean Polyphaga
P450	EOG805VG7	148.37	153.21	143.35	148.10	137.95	34.13	34.47
CE	EOG87DCWX	143.23	143.77	143.63	138.91	141.15	6.55	18.78
CE	EOG8KD911	87.08	91.23	82.90	89.10	79.74	0.89	2.86
CE	EOG876NDC	80.64	85.67	80.08	87.42	77.23	1.72	3.48
GST	EOG87WR3Z	86.24	87.76	76.05	74.40	72.12	1.71	3.16
GST	EOG81RS7Z	108.77	114.44	109.19	113.76	103.79	6.85	11.69
GST	EOG85F05D	117.62	115.88	111.53	107.62	106.32	5.69	9.16
CYS	EOG8JDKNM	91.85	91.62	89.74	85.66	88.25	1.80	3.78

Category	Positive/total OGs	Category enrichment FDR
P450	1/22	0.36268
CE	3/19	0.00252
GST	3/6	0.00016
CYS	1/7	0.16627
UGT	0/4	1.00000
SER	0/4	1.00000
ABC	0/28	1.00000
GH	0/1	1.00000
Category	Positive/total OGs	Candidate vs. background enrichment p value
Background (Polyphaga)	88/9720	0
Candidates (Polyphaga)	8/91
Background (Adephaga)	21/9720	1
Candidates (Adephaga)	0/91

Dataset heterogeneity

For the comparison of gene repertoires between the two groups to be unbiased, the gene content of all analyzed species should be of similar accuracy and completeness. The number of predicted proteins for the genomic resources for each beetle species (Table 1, mean 15,977 and standard deviation 3748) was within the range expected for insects (see [44]). The average total gene count for Adephaga species (all transcriptomes) was about 4200 fewer than for Polyphaga, which include 2 genomes with more than 22,000 genes. This difference in average gene counts is reduced to just 1384 when considering only genes assigned to the 9720 OGs selected for the analysis. Our conservative orthology filtering therefore ensured that the comparisons focused on gene families with reliably traceable evolutionary histories that span both groups of beetles. Secondly, assessments of completeness showed that the majority of the datasets contained more than 90% of complete BUSCOs (Fig. 1, Table 1). While the dynamically evolving families that are the focus of this study are clearly not universal single-copy orthologs, the high levels of BUSCO completeness support the assumption that the datasets represent good coverage of the species’ gene content. Furthermore, the transcriptomes were sequenced from adults collected from their natural habitats, so members of the gene families that are the focus of our study were likely being actively transcribed at the time of sampling. Re-analyses of our data that exclude the two adephagan beetle species with fewer than 80% complete BUSCOs reduced the power of the model tests but nevertheless still identified the three GST OGs that favor a model with higher optima for Polyphaga (see Additional file 1: Table S2). Three of the adephagan transcriptomes showed more than 10% of duplicated BUSCOs, which could have arisen from suboptimal filtering of the transcriptomes, i.e., failure to remove alternative transcripts of the same gene. While such potentially inflated gene counts for these adephagans might prevent the identification of some true expansions in Polyphaga, they do not invalidate those that were identified. Conversely, transcriptome assemblies might collapse very similar paralogs into a single transcript and thereby underestimate true gene counts. Our gene trees (Fig. 2 and Additional file 1: Figures S3-S9) nevertheless show several examples of closely related paralogs from the transcriptomes, indicating that they can be successfully recovered. Finally, half of the OGs representing positive results showed a higher mean value for polyphagan transcriptomes than genomes, including the three GST OGs (Additional file 1: Table S1), and explicitly testing for effects due to using both genome and transcriptome data for the species of Polyphaga, by performing a modified OUwie analysis with data type as the regime under selection, identified only one CE (EOG8KD911) for which the favored model linked gene family expansion to species with genomes (see Additional file 1: Table S3).

Candidate OG identification

The annotation strategy was designed to link OGs to candidate gene families based on manually selected keywords used to filter sequence search results, as well as Pfam and InterPro identifiers (Table 2), with the aim of excluding false positives (see the “Methods” section). This conservative strategy may not have fully captured all possible candidate OGs, which would therefore have remained in the background set of OGs that were used as controls. For example, we identified nine GST OGs (six were retained as candidates after filtering) while ten subclasses have been identified in arthropods [45]. While the strict (conservative) strategy we employed to identify candidate OGs may have resulted in an underestimate of the extent of the observed effects, this does not invalidate those that were identified. In addition, filtering the OGs to retain only those with genes from both Adephaga and Polyphaga excluded from the analyses any genes that were specific to either suborder. These might include genes with key roles in phytophagy, e.g., enzymes acquired by horizontal gene transfer identified from the A. planipennis, A. glabripennis, and D. ponderosae genomes [21]. While acknowledging their importance, here, we explicitly tested for adaptive LSE in one lineage vs. the other so gene evolutionary histories were required to span the two suborders and thus be traceable to their last common ancestor.

The more speciose Polyphaga exhibit more dynamic gene repertoire evolution

The loss (μ) and gain (λ) values reported by CAFE on all 9720 OGs are consistent with assessments of other insect clades [46, 47]. Although the overall gain rate is slightly higher than the loss rate, the number of OGs losing genes reported by CAFE at each individual node is generally larger than the number of OGs with gains. This is reconciled by considering that across Coleoptera, many OGs lost a few genes while few families gained many genes. As most OGs display a low number of genes per species, i.e., they are evolving under “single-copy control” [48], losing more than one ortholog per species is understandably rare, while there is no theoretical limit for an OG to gain new members. Comparing the two clades, Polyphaga has a higher rate of gene gain and six times more OGs with gains, and while the Adephaga rate of gene loss is higher, Polyphaga have 1.5 times more OGs that have experienced gene losses. Importantly, these rates were estimated using a time-calibrated ultrametric species phylogeny with no significant difference in node ages between the two suborders. Hence, the gene repertoires of Polyphaga exhibit a more dynamic evolutionary history with more gains (rate) in more OGs (counts) and fewer losses (rate) spread out over more OGs (counts). It is possible that this greater dynamism may be generally linked to the greater species richness of Polyphaga, with no specific role for phytophagy underpinning this trend. However, among the candidate OGs for detoxification and digestion, there are also more gains in Polyphaga and, in contrast to the background, fewer losses. Thus, both gain and maintenance are higher for candidates in Polyphaga, which is consistent with a key role for phytophagy in driving dynamic gene repertoire evolution, and particularly LSEs.

Support for adaptive expansions of gene families involved in detoxification in polyphagan beetles

In addition to observing more expansions among candidate OGs in the suborder Polyphaga, the positive results from the OUwie analysis support the hypothesis that selective pressures drive detoxification enzymes towards larger gene family sizes. This is especially pronounced for GSTs, for which half of the OGs tested positive and for which a significant enrichment compared to the background was found. The importance of GSTs in dietary shifts to phytophagy has been noted in mustard-feeding flies, where duplicated GSTs involved in the mercapturic acid pathway showed signatures of positive selection [49]. Our results therefore suggest that comparable phenomena have been acting at the level of polyphagan beetles. The CEs also show a statistically significant enrichment for positive results compared to the background, further supporting the diet detoxification hypothesis. The other positive results include a P450 OG and a CYS OG, neither of which led to a category enrichment compared to the background. The P450 OG is by far the largest among the positive results (Additional file 1: Figure S3), highlighting the importance of P450s in beetle (and generally insect) physiology with diverse roles beyond detoxification, e.g., hormone biosynthesis [50]. However, while considered as significantly expanded and under selection by our model, the actual mean values in the suborders are not dramatically different. Importantly, the enrichment of positive results among candidates still holds if this OG is excluded (see Additional file 1: Supplementary Results). The involvement of P450s in many other processes may explain why a broader difference between the suborders was not identified. Apart from one positive result among the cysteine proteases (no significant category enrichment), our study did not highlight additional expansions in other digestive enzymes or in transporters within a suborder. The lack of evidence for expansion in Polyphaga with respect to ABC transporters, which is the candidate functional category encompassing the largest number of OGs, may indicate that the ancestral diversity of transporters was sufficient for maintaining the excretion of toxins, despite variations in the substrates imposing a selective pressure on early stages of the detoxification pathway. Alternatively, if such pressure were acting on later stages of the pathway, i.e., transporters, its strength could have been too low for the detection power of our methods and data, unlike for GSTs or CEs.

Data sources

This study included 6 genomes and 13 transcriptomes representing a balanced sampling of polyphagan and adephagan beetles, along with 1 representative of the sister group to Coleoptera, Strepsiptera, to root the species phylogeny. Annotated gene sets from 4 genomes were sourced from the i5k pilot project datasets [53] (Anoplophora glabripennis v0.5.3 [21], Leptinotarsa decemlineata v0.5.3 [37], Onthophagus taurus v0.5.3, Agrilus planipennis v0.5.3) and 2 were independently published Dendroctonus ponderosae Ensembl Metazoa v1.0 [36] and Tribolium castaneum Ensembl Metazoa v3.22 [39]. One Polyphaga transcriptome, Laparocerus tessellatus (Additional file 1), was sequenced for this project [38], and the others were provided by the 1KITE project (Additional file 1, [34, 35, 54]). A detailed list is presented in Table 1.

Coding sequence predictions, transcriptome, and genome quality assessments

Coding sequences and peptide sequences were predicted from all transcriptomes using TransDecoder (v2.0.1 https://transdecoder.github.io [last accessed May 8, 2019]) along with a custom python script to retain the best-scoring entry among overlapping predictions. The coding sequences and peptide sequences from the genomes were retrieved from their official annotated gene sets. All genome and transcriptome gene sets were assessed using BUSCO (v2.0, python 3.4.1, dataset insecta_odb9/2016-10-21, mode proteins) [32]. This tool identifies near-universal single-copy orthologs by using hidden Markov model profiles from amino acid alignments. CD-HIT-EST v4.6.1 [55] was run on the protein sequences with a 97.5% identity threshold to ensure that all species datasets were filtered to select a single isoform per gene.

Orthology delineation

The OrthoDB [40] hierarchical orthology delineation procedure was employed to predict orthologous protein groups (OGs). Briefly, protein sequence alignments are assessed to identify all best reciprocal hits (BRHs) between genes from each pair of species, which are then clustered into OGs following a graph-based approach that starts with BRH triangulation. The annotated proteins from the genomes of A. planipennis, O. taurus, and all transcriptomes were mapped to OrthoDB v8 at the Arthropoda level (with 87 species including 4 of the beetles with sequenced genomes). Mapping uses the same BRH-based clustering procedure but only allows genes from mapped species to join existing OGs. These OGs were then filtered to identify the 9720 OGs with representatives from both Polyphaga and Adephaga to focus the study on OGs with evolutionary histories traceable to the last common ancestor of all the beetles, i.e., 5188 OGs with genes from only 1 of the 2 suborders were removed.

Time-calibrated species phylogeny

To build an ultrametric phylogeny required for the CAFE analyses, the maximum likelihood molecular species phylogeny was first estimated based on the concatenated superalignment of orthologous amino acid sequences from each of the datasets. Protein sequences of single-copy BUSCO genes and the best-scoring duplicated genes present in all species were individually aligned for each set of BUSCO-identified orthologs using MAFFT with the --auto parameter [56], and each result was manually reviewed to exclude poor-quality alignments. Four hundred and five alignments were retained out of 436 and concatenated into a superalignment, partitioned according to the best model for each set of orthologs using aminosan 1.0.2015.01.23 [57]. RAxML v8.1.2 (-f a -m PROTGAMMA -N 1000) [58] was used to compute the maximum likelihood tree. The monophyly of Geadephaga and Hydradephaga was constrained to match the generally accepted resolution of Adephaga (as [6]). The chronos function of the R package ape (v3.4 on R 3.2.1, relaxed model) [59] was used to obtain an ultrametric tree, and the tip to root length was adjusted to match the approximately 250 million-year evolutionary history of crown group Coleoptera [6].

Functional annotation and definition of candidate genes

InterProScan was run on all species' protein sets (-appl Pfam --goterms, 5.16.55) [60] to identify protein families. Additionally, blastp 2.3.0 [61, 62] was run against uniref50 (version June 22, 2016; [63]) with an e value cutoff of 1e−20. An OG was included in the set of candidate OGs when it had a match to both the uniref50 clusters and Pfam families [64] or gene ontologies [65] as detailed in Table 2.

CAFE analysis

The number of genes in OGs for each species was counted. All candidates and remaining (control) OGs were pooled together and processed with CAFE 3.1 [41], to infer gene family evolution in terms of gene gains and losses. First, the python script provided by CAFE was used to estimate the error in our dataset. The CAFE software was then run using the mode in which the gain and loss rates are estimated together (λ) and a second mode in which they are estimated separately (gains = λ, losses = μ). The more complex model was retained as it reached a significantly better score (− 199,989 for a single estimated parameter and − 199,981 for two distinct estimated parameters, 2× delta log-likelihood = 16, chi-squared distribution, df = 1). For the entire analysis, the CAFE overall p value threshold was kept at its default value (0.01). To run CAFE on each suborder separately, the newick file was pruned to retain only the required species using newick utils 1.1.0 [66].

Evolutionary models

To evaluate adaptive OG expansion, the likelihood of the count data was tested by optimizing parameters considering two methods provided by the OUwie R package v1.51 [43, 67]. First, a Brownian motion (BM) approach was used, which assumes no selection and thus differences result from a stochastic process whose rate is estimated. Second, Ornstein-Uhlenbeck (OU) models were used. They take into account an optimal family size that is obtained by selective pressure. Two groups were defined in the phylogeny, namely Polyphaga and Adephaga, to which the two different regimes to consider were assigned, plus a third regime to the root. This represents a simplified scenario allowing for the comparison of gene contents between one group and the other rather than attempting to estimate “levels” of phytophagy or zoophagy across the phylogeny. The models BM1 (Brownian motion with a single rate for the whole tree), BMS (Brownian motion with different rates for each group), and OU1 (selection towards the same optimum for both groups) were optimized as null hypotheses (H0) and compared to OUM (selection towards two optima, same variance) and OUMV (selection towards two optima, two variances) models as alternative hypotheses (H1). The Akaike information criterion corrected for small sample size (AICc) [68] was used to compare the models, and an AICc > 2 between the best H0 and the best H1 model was considered as significant to prefer the H1 model.

Statistical enrichment

All results for candidates and controls were pooled together to obtain a background distribution of positive and negative results. Positive results are those OGs that passed the OUwie analysis, and negative results are all of the 9720 OGs that did not obtain a significant overall CAFE p value or did not pass the OUwie analysis. Then, 100,000 random draws (using the R function sample, without replacement) having the sample size of the candidate category to test for enrichment were taken from the background, and the significant outcomes for Polyphaga and Adephaga were counted. A p value was calculated for each group as follows: the number of random draws reaching the amount of significant outcomes found for the candidate category, or more, divided by 100,000. Additionally, the multiple tests conducted on each individual candidate category were corrected for false discovery rate (FDR) using the R p.adjust function (method BH, Benjamini Hochberg).

Gene trees

The alignments for the gene trees for the eight OGs that tested positive for adaptive expansions were produced using MAFFT with the --auto parameter. The gene trees were computed with RAxML v8.1.2 (-f a -m PROTGAMMALGF -N 100) and plotted with EvolView [69, 70].

Further methods details are presented in Additional file 1: Supplementary Methods, and a chart summarizing the main steps of the analysis is available in Additional file 1: Figure S10.