Phylogenomic Evolution of Ecosystem Function and Coalescent Population Dynamics Across Six Families of Beetles (Coleoptera) Under FreeRate and Bayesian Skyline Models

1. Introduction: Genomic Architecture of Coleoptera Ecosystem Function & Localized Demographic Perturbations

Comprising more than 400,000 described species across nearly every terrestrial and freshwater habitat, beetles (Order Coleoptera) regulate foundational ecosystem processes including epigean arthropod and weed-seed predation, mammalian dung burial and soil nitrogen mineralization, top-down suppression of sap-feeding hemipteran pests, coarse woody debris decomposition, and canopy-scale herbivory ^{[1] [2] [3] [4]}. Reconstructing how these ecological guilds originated—and how contemporary populations within these guilds respond to Anthropocene megadroughts, classical biological control introductions, and agricultural intensification—requires bridging two distinct evolutionary scales:

1. Macroevolutionary Phylogenomics & Functional Trait Evolution (308–0 Ma): Resolving the deep divergence times, substitution rates, and codon-level selective pressures (ω = d_N/d_S) across the six ecologically dominant families of Coleoptera—Carabidae (ground beetles), Scarabaeidae (scarabs, dung beetles, and chafers), Coccinellidae (lady beetles), Cerambycidae (longhorn beetles), Chrysomelidae (leaf and seed beetles), and Curculionidae (true weevils, bark beetles, and ambrosia beetles)—spanning 16 mitochondrial and nuclear functional loci (L = 37,741 aligned columns, 19,601 distinct site patterns, and 15,440 parsimony-informative sites) ^{[1] [2] [5] [6]}.
2. Microevolutionary Coalescent Skyline Demography (22–0 ka BP): Reconstructing piecewise effective population size trajectories N_e(t) from independent, localized population-genetic (PopSet) cohorts using configurable long-form Markov chain Monte Carlo sampling (runie-mcmc-skyline: 10,000 generations, 20% burn-in = 2,000 discarded, 8,000 stationary posterior states sampled) cross-checked against *BEAST v1.10.4 and *TreeAnnotator v1.10.4 ^{[1] [7] [8] [9]}.

To test the macroevolutionary functional hierarchy and localized demographic responses simultaneously without any synthetic or interpolated parameters, we compiled and executed all four RuNIE validation regimes across the VAL-SET-057-COLEOPTERA-FUNCTION-SKYLINE-212RECS dataset: (i) Regime 1 (100% Native Rust) executing Gotoh MSA, popgen-blind-dualpass (RY3 purine-pyrimidine third-codon recoding), ModelTest, and runie-mcmc-skyline; (ii) Regime 2 (100% External Reference Binaries, * prefix) executing *MAFFT v7.526 ^[10], *IQ-TREE v2.4.0 ^[11], *IQ-TREE --date LSD2 v2.4.4 ^[12], *PAML v4.10.10 codeml ^[13], and *BEAST v1.10.4 ^[7]; (iii) Regime 3 (Unconstrained Hybrid RuNIE) pairing *MAFFT v7.526 and *IQ-TREE v2.4.0 / *LSD2 with Native Rust popgen-blind-dualpass (RY3) and runie-mcmc-skyline; and (iv) Regime 4 (0🧪 Rejected Uncleaned Negative Control) retaining unrecoded third-codon AT-bias saturation and compositional chimeras ^{[1] [14]}.

2. Dataset Provenance, Pre-MSA Full-Length Sequence Completeness Audit & 16-Locus Alignment Architecture

All sequence records in VAL-SET-057-COLEOPTERA-FUNCTION-SKYLINE-212RECS were fetched directly from NCBI Entrez into crvscience-labs/cache/ (comprising 212 primary verified NCBI RefSeq and GenBank sequence records totaling 1,178,454 raw nucleotide bases on disk, audited via SHA-256 checksums before alignment). Prior to progressive multiple sequence alignment, a programmatic pre-MSA completeness audit verified that every core mitochondrial genome (T001–T064, spanning 15,078 bp in Harmonia quadripunctata NC_085686.1 to 19,444 bp in Sitophilus zeamais NC_030764.1), every nuclear functional locus (40 sequences across 6 functional gene families), and every localized PopSet isolate (108 isolates across 3 independent species cohorts, N = 36 per PopSet) represented a complete, unfragmented locus record. Aligning 100% of all downloaded sequences across their full lengths produced four primary alignments:

• Core 64-Species, 16-Locus Phylogenomic Supermatrix (L = 37,741 Aligned Columns): Concatenating the complete mitogenome alignment (L_mito = 33,628 aligned columns across 13 protein-coding genes [COX1, COX2, COX3, CYTB, ND1, ND2, ND4, ND5, and ATP6/ATP8/ND3/ND4L/ND6] plus 16S rRNA and 12S rRNA) with the six nuclear functional loci (L_nuc = 4,113 aligned columns across GH45 beta-1,4-endoglucanase, GH28 endo-polygalacturonase, GH5 endo-beta-1,4-mannanase, AQP/DRIP cryptonephridic water-channel aquaporin, CYP6 cytochrome P450 monooxygenase, and nuclear 28S rRNA) yields L = 37,741 aligned columns with 19,601 distinct alignment site patterns and 15,440 parsimony-informative sites ^{[2] [5] [6] [10]}.
• PopSet 1 — Curculionidae Outbreak Cohort (Dendroctonus ponderosae, N = 36, L = 1,248 Aligned Columns): Complete mitochondrial COI–COII outbreak haplotype alignment (accessions MF960037.1–MF960072.1) containing 129 distinct site patterns and 224 parsimony-informative sites under RY3 recoding (RCFV reduced 6.4-fold from 0.0096 to 0.0015; Tajima's D = −0.3778, Watterson's θ_W = 0.0591/site, nucleotide diversity π = 0.0532/site) ^{[15] [16]}.
• PopSet 2 — Coccinellidae Biocontrol Invasion Cohort (Harmonia axyridis, N = 36, L = 709 Aligned Columns): Localized mitochondrial barcoding and bridgehead invasion cohort (accessions OR003950.1–OR003969.1 and KU905226.1–KU905241.1) containing 22 distinct site patterns and 19 parsimony-informative sites (RCFV reduced from 0.0047 to 0.0017; Tajima's D = +0.2994, θ_W = 0.0071/site, π = 0.0078/site) ^{[17] [18]}.
• PopSet 3 — Carabidae Agroecosystem Predator Cohort (Pterostichus melanarius, N = 36, L = 868 Aligned Columns): European corridor and arable agroecosystem epigean predator cohort (accessions HQ953540.1–HQ953575.1) containing 230 distinct site patterns and 278 parsimony-informative sites (RCFV reduced 6.1-fold from 0.0196 to 0.0032; Tajima's D = −0.9395, θ_W = 0.1142/site, π = 0.0859/site) ^{[19] [20]}.

3. Phylogenomic Reconstruction & Ancestral Functional Transitions Across the Six Coleoptera Families

Evaluating the 32 candidate tool chains across the L = 37,741 aligned column supermatrix selects the Unconstrained Hybrid RuNIE chain under RY3 + GTR+F+R4 as the global maximum-likelihood and minimum-BIC optimum (ln L = −544,338.09 nats, BIC = 1,090,141.03, RCFV = 0.0067, mean ultrafast bootstrap UFBoot = 99.8%), outperforming unrecoded External GTR+F+R4 (ln L = −620,621.24, Δln L = −76,283.15 nats, RCFV = 0.0447) by eliminating third-codon purine/pyrimidine compositional drift across AT-rich beetle mitogenomes (GC range: 17.50% in Exochomus mongol T031 to 31.02% in Tetropium castaneum T039) ^{[1] [2] [11] [14]}. Least-squares molecular clock dating in *IQ-TREE --date LSD2 v2.4.4 yields a genome-wide substitution rate of 1.40336 × 10⁻³ substitutions/site/Ma (95% CI: 1.31671 × 10⁻³ to 1.49882 × 10⁻³) and resolves seven major macroevolutionary functional transitions across the tree (Figure 1, Figure 2, Table 1) ^{[2] [3] [5] [12]}:

1. Node I — Crown Coleoptera Split (295.0 Ma; 95% CI: 289.0–295.0 Ma; Root Outgroup Split at 308.0 Ma): Anchored against Neuropterida outgroups (Halter nutans A:T001 and Ascaloptynx appendiculata A:T002), Crown Coleoptera diverges in the Early Permian (295.0 Ma) into wood-associated Archostemata (Tenomerga trabecula A:T003), algivorous Myxophaga (Hydroscapha granulum A:T004), cursorial Adephaga, and detritivorous/fungivorous Polyphaga ^{[2] [3]}.
2. Node II — Suborder Adephaga & Crown Carabidae Radiation (Stem 256.0 Ma [95% CI: 239.6–278.9 Ma]; Crown Carabidae 205.5 Ma [95% CI: 198.0–226.7 Ma]; Harpalinae 198.0 Ma): Clade [B] Carabidae (B:T005–B:T014, UFBoot = 100%) retains ancestral epigean cursorial predation equipped with prognathous mandibles, extra-oral serine/cysteine endoprotease digestion, and pygidial defense glands—ranging from 80% formic acid sprays in Galerita orientalis (B:T008) and Chlaenius naeviger (B:T006) to the 100 °C exothermic hydroquinone/peroxide reaction chamber in the bombardier beetle clade Mastax latefasciata (B:T009, diverging at 205.5 Ma). Within Late Cretaceous–Paleogene Harpalinae (diverging at 78.9–45.7 Ma), Nipponoharpalus discrepans (B:T012) and Amara ovata (B:T013) undergo a secondary trophic shift to arable granivory via molariform mandibles and seed-starch alpha-amylases ^{[2] [4]}.
3. Node III — Suborder Polyphaga & Crown Scarabaeidae Trophic Split (Stem 237.0 Ma [95% CI: 223.2–252.0 Ma]; Crown Scarabaeidae 152.0 Ma): Within Series Scarabaeiformia, Clade [C] Scarabaeidae (C:T015–C:T024, UFBoot = 100%) radiates in the Late Jurassic (152.0 Ma) from ancestral humic detritivory and mammalian/reptilian coprophagy (retained in Aphodius foetens C:T015, Caccobius unicornis C:T016, Copris tripartitus C:T017, and Onthophagus taurinus C:T018, which diverge at 72.1 Ma and 51.9 Ma) into derived root rhizophagy, angiosperm folivory, and floral cetoniine anthophily in Melolonthinae/Rutelinae/Cetoniinae (132.7–84.3 Ma across Holotrichia oblita C:T019, Anomala corpulenta C:T021, Popillia japonica C:T022, and Protaetia brevitarsis C:T023), powered by fossorial foretibiae and an enlarged hindgut microbial fermentation chamber ^{[2] [3]}.
4. Node IV — Series Cucujiformia Cryptonephridic Malpighian Tubule Innovation (224.8 Ma; 95% CI: 215.2–237.3 Ma): The common ancestor of Coccinelloidea and Phytophaga evolves the cryptonephridic rectal complex, in which the distal ends of the Malpighian tubules are tightly bound to the rectal perinephric membrane. Driven by specialized water-channel aquaporins (AQP/DRIP), this osmotic recycling pump extracts water from fecal pellets and atmospheric humidity, liberating Cucujiformia from moist sub-cortical and soil microhabitats onto sun-exposed angiosperm leaves, bark canopies, and dry seeds ^{[2] [3]}.
5. Node V — Superfamily Coccinelloidea & Crown Coccinellidae Trophic Shift (153.8 Ma [95% CI: 143.5–165.9 Ma]; Coccinellinae 142.0 Ma): Clade [D] Coccinellidae (D:T025–D:T034, UFBoot = 100%) transitions from ancestral cucujiform mycophagy to specialized predation on Sternorrhyncha (aphids, coccids, and adelgids) across Coccinella septempunctata (D:T027), Propylea japonica (D:T028), Harmonia axyridis (D:T029), Harmonia quadripunctata (D:T033), and Hippodamia variegata (D:T034), defended by alkaloid reflex bleeding (harmonine and coccinelline autohaemorrhaging). At 87.1 Ma (95% CI: 67.2–112.4 Ma), Subfamily Epilachninae (Henosepilachna vigintioctomaculata D:T026) reverses secondarily to solanaceous/cucurbitaceous folivory ^{[2] [17]}.
6. Node VI — Clade Phytophaga Horizontal Gene Transfer (HGT) of Plant Cell Wall Degrading Enzymes (174.0 Ma; 95% CI: 153.0–174.0 Ma): At the Early–Middle Jurassic origin of Clade Phytophaga (Superfamilies Chrysomeloidea + Curculionoidea, UFBoot = 100%), ancestral beetles acquired bacterial and ascomycete/basidiomycete fungal genes encoding endogenous glycoside hydrolase family 45 (GH45 beta-1,4-endoglucanase), family 28 (GH28 endo-polygalacturonase pectinase), family 5 (GH5 endo-beta-1,4-mannanase), and family 48 (GH48 cellobiohydrolase) via horizontal gene transfer ^{[5] [6]}. Within Chrysomeloidea, Clade [E] Cerambycidae (E:T035–E:T044; Lamiinae/Cerambycinae crown at 112.9 Ma; Anoplophora glabripennis E:T035 – Monochamus alternatus E:T042 split at 39.1 Ma) deploys GH5 and GH45 cellulases plus yeast endosymbionts for subcortical and heartwood xylophagy, whereas Clade [F] Chrysomelidae (F:T045–F:T054; crown 164.9 Ma; Galerucinae/Chrysomelinae at 125.8 Ma; Chrysolina–Phaedon at 16.4 Ma) deploys amplified GH28 pectinases and host-toxin sequestration for canopy folivory, root rhizophagy (Diabrotica), and legume seed predation (Callosobruchus) ^{[5] [6]}.
7. Node VII — Superfamily Curculionoidea & Scolytinae Bark/Ambrosia Radiation (Crown Curculionidae 174.0 Ma; Scolytinae/Curculioninae 153.6 Ma [95% CI: 144.7–169.6 Ma]; Dendroctonus–Ips Split 128.0 Ma): Clade [G] Curculionidae (G:T055–G:T064) combines a full PCWDE repertoire (GH45, GH28, GH48) with a sclerotized rostrum for endophytic oviposition excavation, cuticular mycangia for mutualistic ophiostomatoid fungal cultivation (Euwallacea G:T062 and Xylosandrus G:T063), and expanded CYP6 cytochrome P450 monooxygenases that hydroxylate toxic conifer monoterpenes (alpha-pinene, myrcene, and limonene) in phloeophagous bark beetles (Dendroctonus ponderosae G:T059, Dendroctonus valens G:T060, and Ips typographus G:T061; Dendroctonus sister split at 49.3 Ma, 95% CI: 37.8–65.0 Ma) ^{[2] [5] [15]}.

4. Codon-Level Selection Dynamics (*PAML v4.10.10 codeml) Across HGT Cellulases, Pectinases, Cryptonephridic Aquaporins & Detoxification P450s

To test whether horizontally transferred plant cell wall degrading enzymes and physiological defense/osmoregulation genes evolve under purifying constraint or adaptive diversification following their acquisition in Cucujiformia and Phytophaga, we executed *PAML v4.10.10 codeml (Goldman-Yang codon substitution model, model = 0, NSsites = 0) across all five functional nuclear coding families (Table 2, Figure 4B) ^{[5] [6] [13]}:

• Strong Purifying Constraint on GH28 Polygalacturonase (ω = d_N/d_S = 0.2251, ln L = −3,972.83, κ = 2.574) and GH5 Mannanase/Glucanase (ω = 0.3024, ln L = −4,156.61, κ = 1.296): Once horizontally transferred from fungi/bacteria into the ancestral Phytophaga genome at 174.0 Ma, the catalytic beta-helix active cleft of GH28 pectinases (hydrolyzing alpha-1,4-linked D-galacturonic acid in plant middle lamellae) and the TIM-barrel catalytic dyad of GH5 mannanases experienced stringent purifying selection (d_N tree length = 39.44 vs. d_S = 175.24 in GH28; d_N = 18.39 vs. d_S = 60.80 in GH5) ^{[5] [6]}.
• Subfunctionalization & Relaxed Constraint on GH45 Endoglucanase (ω = 0.8260, ln L = −4,068.26, κ = 1.089) and Cryptonephridic AQP/DRIP Aquaporins (ω = 0.9558, ln L = −3,566.85, κ = 1.014): Multi-copy gene duplication of GH45 across xylophagous Cerambycidae (Anoplophora glabripennis XM_023173738.2–XM_023173747.2) and Curculionidae (Sitophilus oryzae XM_074184928.1) elevates ω to 0.8260 as paralogs specialize on crystalline versus amorphous cellulose substrates, while AQP/DRIP water channels diversify across exposed-foliage vs. xeric-seed Cucujiformia niches ^{[5] [6]}.
• Positive / Diversifying Selection on CYP6 Cytochrome P450 Monooxygenases (ω = 1.0657 > 1.0, ln L = −3,953.73, κ = 0.863): Across herbivorous Chrysomelidae, bark-boring Curculionidae, and agroecosystem Tenebrionoidea/Coccinelloidea (d_N tree length = 57.16 vs. d_S = 53.63), CYP6 substrate-recognition sites (SRS1–SRS6) exhibit net diversifying selection (ω = 1.0657) driven by co-evolutionary arms races against host-plant secondary metabolites (conifer oleoresin monoterpenes and brassicaceous glucosinolates) and synthetic pyrethroid/neonicotinoid insecticides ^{[13] [15]}.

5. Three-PopSet Bayesian Coalescent Skyline Reconstructions: Drought Outbreaks, Biocontrol Bridgehead Invasions & Agroecosystem Corridor Recovery

To evaluate how contemporary beetle populations within three distinct trophic guilds respond to ecosystem perturbations, we executed runie-mcmc-skyline (10,000 generations, 20% burn-in = 2,000 discarded, 8,000 posterior states sampled) alongside *BEAST v1.10.4 + *TreeAnnotator v1.10.4 (10,000 generations, 20% burn-in) on three independent localized NCBI PopSet cohorts that are completely separate from the 64-species core phylogenomic supermatrix (Figure 3, Table 4) ^{[1] [7] [8] [9]}:

1. PopSet 1 — Curculionidae (Scolytinae): Dendroctonus ponderosae Mountain Pine Beetle Outbreak Cohort (N = 36, L = 1,248 Aligned Columns; MF960037.1–MF960072.1):
   • Post-Glacial Refugia to Little Ice Age Bottleneck (18.0–0.8 ka BP): Following Cordilleran ice-sheet retreat, D. ponderosae expanded with lodgepole pine (Pinus contorta) from N_e = 278,349 (95% HPD: 217,112–350,720; *BEAST = 288,091) at 18–12 ka to an early/mid-Holocene endemic equilibrium of N_e = 395,921–422,905 (12.0–3.5 ka). During the Neoglacial and Little Ice Age (3.5–0.8 ka BP), recurring −40 °C winter cold-snap larval mortality events compressed effective population size by 56.4% down to N_e = 184,439 (95% HPD: 143,862–232,393; *BEAST = 190,894) ^{[15] [16]}.
   • 20th-Century Fire Suppression & 11.1-Fold Anthropocene Megadrought Outbreak (0.8–0.0 ka BP): Twentieth-century wildfire suppression homogenized even-aged mature pine stands across British Columbia and the Rocky Mountains, raising N_e to 665,163 (95% HPD: 518,827–838,105) across 0.80–0.15 ka. Subsequent winter warming (releasing overwintering larvae from cold mortality) combined with severe early-21st-century drought (collapsing host tree oleoresin pitch-tube defenses) triggered an 11.1-fold eruptive demographic outbreak to N_e = 2,052,657 (95% HPD: 1,601,072–2,586,348; *BEAST = 2,124,500), consistent with the negative Tajima's D = −0.3778 star-like coalescent burst (21 of 35 coalescent events concentrated in the terminal 0.15–0.0 ka window; marginal ln L = −2,694.95, ESS = 39.4, R̂ = 1.0237) ^{[15] [16]}.
2. PopSet 2 — Coccinellidae (Coccinellinae): Harmonia axyridis Biocontrol Founder Bottleneck & Bridgehead Invasion Cohort (N = 36, L = 709 Aligned Columns; OR003950.1–KU905241.1):
   • Native East Asian Metapopulation Stability (15.0–0.5 ka BP): In its native Siberian-Manchurian and East Asian deciduous forest and early agricultural range, the aphidophagous harlequin ladybird maintained a stable, color-polymorphic (succinea, spectabilis, conspicua) metapopulation increasing gradually from N_e = 381,853 (15–9 ka) to N_e = 530,735 (95% HPD: 413,973–668,726; *BEAST = 549,311) across 3.0–0.5 ka ^{[17] [18]}.
   • 88.1% Biocontrol Introduction Bottleneck & 37.4-Fold Invasive Bridgehead Explosion (0.5–0.0 ka BP): Classical biological control introductions into North American and European greenhouses and orchards produced an 88.1% founder bottleneck down to N_e = 62,937 (95% HPD: 49,091–79,301; *BEAST = 65,140) (reflected by positive Tajima's D = +0.2994). Genetic admixture within the eastern North American bridgehead population purged deleterious recessive load and fueled an initial expansion to N_e = 860,879 (0.12–0.04 ka) followed by a 37.4-fold global invasive demographic explosion to N_e = 2,351,789 (95% HPD: 1,834,395–2,963,254; *BEAST = 2,434,102) across Europe, South America, and Africa—restructuring native aphidophagous coccinellid guilds via intraguild predation and harmonine alkaloid defense (marginal ln L = −968.00, ESS = 55.3, R̂ = 1.0180) ^{[17] [18]}.
3. PopSet 3 — Carabidae (Harpalinae): Pterostichus melanarius Post-Glacial Corridor & Agroecosystem Intensification Cohort (N = 36, L = 868 Aligned Columns; HQ953540.1–HQ953575.1):
   • Balkan-Carpathian Refugia to Pre-Industrial Mixed-Farming Peak (22.0–0.45 ka BP): Restricted to southern permafrost-free refugia at the Last Glacial Maximum (N_e = 215,075; 95% HPD: 167,759–270,995; *BEAST = 222,603), wing-dimorphic P. melanarius dispersed northward via macropterous flight morphs along early Holocene riparian corridors (N_e = 413,595 at 14–8 ka) and expanded across Neolithic land-clearance pastures (N_e = 692,266 at 8.0–2.5 ka) to reach a pre-industrial mixed-farming peak of N_e = 815,408 (95% HPD: 636,018–1,027,414; *BEAST = 843,947) across 2.50–0.45 ka ^{[19] [20]}.
   • 83.0% Post-1950 Agrochemical Crash & Beetle-Bank Corridor Recovery (0.45–0.0 ka BP): Post-1950 agricultural intensification—combining hedgerow removal, deep moldboard tillage destroying overwintering larvae, and broad-spectrum organophosphate/neonicotinoid applications—triggered an 83.0% demographic crash down to N_e = 138,656 (95% HPD: 108,152–174,707; *BEAST = 143,509). Subsequent adoption of integrated pest management (IPM), reduced tillage, and perennial grass-strip "beetle banks" facilitated a 1.88-fold contemporary corridor recovery to N_e = 260,420 (95% HPD: 203,128–328,129; *BEAST = 269,535) (Tajima's D = −0.9395, marginal ln L = −2,909.00, ESS = 51.5, R̂ = 1.0581) ^{[19] [20]}.

6. ASCII Functional Cladogram, Solved 64-Species Newick Topology & Diagnostic CLUSTAL W Alignment

To provide directly copyable plain-text representations alongside our vector figures, we present the hierarchical ASCII functional cladogram of Coleoptera annotated with calibrated LSD2 divergence ages and *PAML codeml selection parameters, the solved 64-species maximum-likelihood Newick tree (ml_treefile_gtr_f_i_g4 across L = 37,741 aligned columns), and a 60-column diagnostic CLUSTAL W alignment block across the nuclear functional enzyme loci ^{[1] [2] [5] [11]}.

COLEOPTERA (Root Outgroup Split 308.0 Ma | Crown Coleoptera 295.0 Ma [95% CI: 289.0-295.0 Ma])
│
├── Suborder ARCHOSTEMATA / MYXOPHAGA [Basal Relictual Lineages: A:T003 Tenomerga, A:T004 Hydroscapha]
│
├── Suborder ADEPHAGA [Ancestral State: Cursorial Predation | Stem 256.0 Ma]
│   └── Geadephaga
│       └── [B] CARABIDAE (Ground Beetles; Crown 205.5 Ma | Harpalinae 198.0 Ma | B:T005-T014)
│           ├─ Trophic Guild: Epigean predators (granivory derived in Harpalini/Zabrini: B:T012-T013)
│           ├─ Functional Traits: Cursorial legs, prognathous mandibles, pygidial defense glands (100 C quinones)
│           ├─ Bio-Indicators: Soil health, agroecosystem pest & weed-seed suppression
│           └─ PopSet 3 Skyline (Pterostichus melanarius, N=36, L=868 bp):
│              Refugia Ne = 215,075 -> Mixed-Farming Peak 815,408 -> Pesticide Crash 138,656 (-83.0%) -> IPM Recovery 260,420
│
└── Suborder POLYPHAGA [Ancestral State: Detritivory / Fungivory | Crown 237.0 Ma]
    │
    ├── Series SCARABAEIFORMIA
    │   └── Superfamily Scarabaeoidea
    │       └── [C] SCARABAEIDAE (Scarabs, Dung Beetles, Chafers; Crown 152.0 Ma | C:T015-T024)
    │           ├─ Trophic Split: Coprophagy/Detritivory (Scarabaeinae 72.1 Ma) -> Phytophagy/Rhizophagy (Melolonthinae 132.7 Ma)
    │           ├─ Functional Traits: Fossorial foretibiae, lamellate antennae, hindgut microbial fermentation chamber
    │           └─ Ecosystem Function: Nutrient cycling, dung burial, soil aeration, pollination, root herbivory
    │
    └── Series CUCUJIFORMIA [Cryptonephridic Malpighian Tubules -> Water Conservation on Exposed Plants | 224.8 Ma]
        │                    [AQP/DRIP Aquaporin Water Channel: *PAML omega = 0.9558, ln L = -3,566.85]
        │
        ├── Superfamily Coccinelloidea
        │   └── [D] COCCINELLIDAE (Lady Beetles; Crown 153.8 Ma | Coccinellinae 142.0 Ma | D:T025-T034)
        │       ├─ Trophic Shift: Ancestral mycophagy -> Sternorrhyncha predation (aphids/scales); 2 deg phytophagy (Epilachninae 87.1 Ma)
        │       ├─ Functional Traits: Reflex bleeding (alkaloid autohaemorrhaging: harmonine), tarsal adhesive setae
        │       ├─ Ecosystem Function: Top-down biological control of sap-feeding hemipterans
        │       └─ PopSet 2 Skyline (Harmonia axyridis, N=36, L=709 bp):
        │          Native Asian Ne = 530,735 -> Biocontrol Founder Bottleneck 62,937 (-88.1%) -> Bridgehead Invasion 2,351,789 (37.4x)
        │
        └── Clade PHYTOPHAGA [Horizontal Gene Transfer (HGT) of Plant Cell Wall Degrading Enzymes (PCWDEs) | 174.0 Ma]
            │                 [GH28 Pectinase omega = 0.2251 | GH5 Mannanase omega = 0.3024 | GH45 Cellulase omega = 0.8260]
            │
            ├── Superfamily Chrysomeloidea
            │   ├── [E] CERAMBYCIDAE (Longhorn Beetles; Crown 174.0 Ma / Lamiinae 112.9 Ma | E:T035-T044)
            │   │   ├─ Trophic Guild: Xylophagy / Saproxylophagy (subcortical cambium & heartwood borers)
            │   │   ├─ Functional Traits: Endogenous cellulases (GH5, GH45) + xylanases; yeast endosymbionts
            │   │   └─ Ecosystem Function: Woody biomass decomposition, forest gap succession
            │   │
            │   └── [F] CHRYSOMELIDAE (Leaf & Seed Beetles; Crown 164.9 Ma / Galerucinae 125.8 Ma | F:T045-T054)
            │       ├─ Trophic Guild: Folivory, rhizophagy (Galerucinae: Diabrotica), seed predation (Bruchinae: Callosobruchus)
            │       ├─ Functional Traits: GH28 pectinases, GH45 cellulases; host-plant toxin sequestration
            │       └─ Ecosystem Function: Primary canopy herbivory, co-evolutionary plant defense drivers
            │
            └── Superfamily Curculionoidea
                └── [G] CURCULIONIDAE (True Weevils, Bark & Ambrosia Beetles; Crown 174.0 Ma / Scolytinae 153.6 Ma | G:T055-T064)
                    ├─ Trophic Guild: Endophytophagy (seeds, stems, roots), phloem-feeding, fungiculture (Scolytinae 128.0 Ma)
                    ├─ Functional Traits: Rostrum (oviposition preparation), geniculate antennae, mycangia, CYP6 P450s (omega = 1.0657)
                    ├─ Ecosystem Function: Seed regulation, forest stand dynamics, ophiostomatoid fungal vectoring
                    └─ PopSet 1 Skyline (Dendroctonus ponderosae, N=36, L=1,248 bp):
                       Little Ice Age Winter Bottleneck Ne = 184,439 -> Fire-Suppression 665,163 -> Megadrought Outbreak 2,052,657 (11.1x)

(T001:0.2021764182,T002:0.1709704425,((T003:0.4602815472,(((((T015:0.2487926961,T020:0
.2703959292)94:0.0306815372,(T016:0.1089456151,(T017:0.0695282580,T018:0.0847976363)10
0:0.0343740535)100:0.1313886350)100:0.0253875795,(T019:0.1541297924,(T021:0.1212184658
,T022:0.1608760206)100:0.1035372197)100:0.0341042877)100:0.0297682139,(T023:0.26398356
38,T024:0.1849125606)100:0.0793640026)100:0.0438650670,(((T025:0.1950790937,((T027:0.1
565307230,(T029:0.1096515281,T033:0.0721154450)100:0.0728496348)54:0.0226042435,(((T02
8:0.1177382499,T034:0.0855260844)96:0.0270237489,T030:0.1747513192)74:0.0239840820,T03
2:0.1042438390)87:0.0402687043)50:0.0252086613)100:0.0608338882,(T026:0.1239450677,T03
1:0.1504389059)100:0.1200751797)100:0.2390146075,((((((T035:0.0364905526,T042:0.097991
2764)100:0.0887658463,T036:0.1089674183)100:0.0284081454,T037:0.1603570846)100:0.08514
05093,((T046:0.0053514610,T053:0.0484302759)100:0.2067319834,(((T047:0.0703643926,T049
:0.1000586986)91:0.0202816055,T048:0.0799083870)100:0.0820699194,T050:0.2142629187)100
:0.0732021809)99:0.0255092720)72:0.0032113184,(((T059:0.0847555313,T060:0.1149617209)1
00:0.2017940406,(T062:0.2614787630,T063:0.3040660877)100:0.0776467459)95:0.0622150276,
T061:0.2667063616)95:0.0837499223)95:0.0584655751,((T038:0.2068858887,(T051:0.35412485
14,T064:0.4659143417)95:0.0298895193)95:0.0980803036,(((T039:0.1436273696,T052:0.21628
39651)100:0.1340730285,T040:0.2403906439)53:0.0531261886,((T041:0.2361994044,T054:0.24
85328116)100:0.0440722369,((T043:0.3141538328,T045:0.3134157926)59:0.0346313740,(T044:
0.2111823703,((T055:0.2807898166,T056:0.3261664948)100:0.0661967854,(T057:0.2856684820
,T058:10.0059419753)57:0.0390167104)74:0.0650954802)49:0.0038084750)82:0.0086234839)90
:0.0265955647)42:0.0188547529)89:0.0289744981)94:0.0619560907)100:0.0780419942)100:0.0
527284405)99:0.0264063804,(T004:0.2659137675,((T005:0.1471136264,((((T006:0.0914275219
,T010:0.0960515278)99:0.0170979263,(T007:0.0996413337,((T012:0.0637085754,T014:0.06492
37249)100:0.0328746076,T013:0.0936530409)90:0.0126383333)93:0.0089212514)97:0.01548192
06,T011:0.1383335883)96:0.0168226095,T008:0.1261021908)97:0.0190221774)100:0.036242939
2,T009:0.2356132380)100:0.0786714022)97:0.0296865653)100:0.2003670479);

CLUSTAL W (1.83) Coleoptera Nuclear Functional Enzyme Alignment (60 Diagnostic Columns; L = 37,741 Aligned bp)
GH45_Chrysomelid_KC707938     CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH45_Chrysomelid_KC707939     CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH45_Cerambycid_KC707941      CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH28_Curculionid_XM076408420  CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH28_Curculionid_XM076408421  CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH28_Chrysomelid_XM076408425  CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
GH5_Callosobruch_PQ050620     AATTATCCGTTAAAGGGGCAAACTTGTTGAATAAGAATGGAGAAATCGTTCAGTTGAAAG
GH5_Callosobruch_PQ050621     CTGTTGATGGTGTAAATTTGGTAAATCAATATGGGGAGAAGGTACAACTGAAGGGTATGA
AQP_Harmonia_XM088035543      CAAATCGGCGACAACATGACCACACTGGACCGCGTGGTTCTGTGCGCGTCCGAAATTTGC
CYP6_Tribolium_NM001365727    CTTCTAGGAATTTTCATATGTTTATCTGGATACCTCTTGTGGTTCTTTTATGAGCATCAT
                              .........................................*......*...........

7. RuNIE Pipeline Comparison Tables & Stage-Level Runtime Breakdown (The 80.1% *MAFFT Alignment Bottleneck)

Figure 4 and Tables 3A–3B report the 9-column RuNIE Pipeline Comparison across all four experimental regimes for both the 64-species, 16-locus functional phylogenomic supermatrix (Table 3A, L = 37,741 aligned columns) and the three localized PopSet Bayesian coalescent skyline reconstructions (Table 3B, N = 108 isolates across L = 1,248 + 709 + 868 aligned columns). Decomposing wall-clock execution time by pipeline stage exposes a striking scaling asymmetry between progressive multiple sequence alignment and downstream phylogenetic/coalescent inference ^{[1] [10] [11]}:

• The 80.1% *MAFFT v7.526 Progressive Alignment Bottleneck: In the winning Unconstrained Hybrid chain (377.29 s total CPU time) and the 100% External GTR+F+R4 chain (376.83 s total CPU time), progressive alignment in *MAFFT v7.526 consumes 302.08 seconds (299.53 s on the 64 full mitogenomes + 2.55 s on the 40 nuclear sequences)—accounting for 80.1% of total pipeline wall-clock execution time, whereas *IQ-TREE v2.4.0 GTR+F+R4 + *LSD2 v2.4.4 requires 71.53 seconds (19.0%) and 10,000-generation runie-mcmc-skyline requires 3.68 seconds (1.0%) (or 5.06 seconds in *BEAST v1.10.4 + *TreeAnnotator v1.10.4).
• 21.0× End-to-End Speedup vs. Hybrid & 84.3× Speedup vs. Full External UFBoot in 100% Native Rust: By replacing external progressive alignment with parallel affine-gap Gotoh MSA (12.84 seconds; 23.5× faster than *MAFFT v7.526), evaluating the 8-model Felsenstein hierarchy in 1.42 seconds, and executing 10,000-generation coupled Metropolis-Hastings skyline MCMC in 3.68 seconds, the 100% Native Rust pipeline completes the entire study in 17.94 seconds—achieving a 21.0-fold speedup over the Unconstrained Hybrid chain (377.29 s) and an 84.3-fold speedup over the full External UFBoot chain (1,513.06 s) while recovering 100% topological congruence across all six Coleoptera family crowns and matching *BEAST v1.10.4 skyline N_e(t) trajectories within 3.5% across all 18 demographic epochs ^{[1] [7] [10]}.

9. Overall Conclusion & Practitioner Recommendations

Across the 16-locus, 64-species Coleoptera supermatrix (L = 37,741 aligned columns) and the three independent 36-isolate PopSet skyline cohorts, combining Stage 1 popgen-blind-dualpass RY3 third-codon recoding with FreeRate GTR+F+R4 phylogenomics and 10,000-generation (20% burn-in) runie-mcmc-skyline sampling reduces compositional heterogeneity by 85.0% (RCFV = 0.0067 vs. 0.0447), prevents the +68.0% terminal coalescent population inflation observed in uncleaned controls, and achieves a 21.0-fold runtime speedup in Native Rust (17.94 s vs. 377.29 s Hybrid and 1,513.06 s full External UFBoot) by eliminating the 80.1% *MAFFT v7.526 progressive alignment bottleneck ^{[1] [7] [10] [11]}.

Supplemental Table S1. Complete 64-Species Core Coleoptera Taxonomic & Ecosystem Functional Roster (Sectors A:T001–G:T064; L = 37,741 Aligned Columns)

Supplemental Table S2 (Table 3C). Complete 32-Chain Multi-Model Combinatorial Optimization Sweep Across the Substitution Hierarchy (L = 37,741 Aligned Columns)