Neuropeptides are evolutionarily ancient signaling molecules that act through diverse G protein-coupled receptor (GPCR) pathways to coordinate essential physiological and behavioral processes in insects. Despite their broad functional importance, the specific roles, receptor interactions, and regulatory mechanisms of individual neuropeptides remain underexplored, limiting our understanding of how these signaling systems operate across developmental and physiological contexts. In this study, we characterized transcript expression of five functionally diverse neuropeptides and one GPCR across Leptinotarsa decemlineata (Colorado potato beetle) tissues and life stages, and subsequently knocked down the glycoprotein hormone subunit alpha 2 (GPA2) and its cognate receptor, leucine-rich repeat-containing GPCR 1 (LGR1), for phenotypic investigation. Expression analyses revealed broad and tissue-specific variation among neuropeptides across developmental stages. No significant differences in mortality, larval weight, or defoliation were detected following dsRNA-mediated knockdown of GPA2 or LGR1 via injection or feeding. However, state-transition modeling identified significant effects of GPA2 and LGR1 knockdown on developmental timing, specifically during pupal-adult and larval-pupal-adult transitions, respectively, with time in stage emerging as a key predictor of overall survival. These findings suggest a potential hormonal role for the GPA2-LGR1 signaling system in regulating developmental progression and stabilizing life-stage transitions in L. decemlineata.

Genomic and transcriptomic analyses confirmed the presence of GPA2 (XP_023027110.1) and its receptor LGR1 (XP_023016362.1) in L. decemlineata. Tachykinin (XP_023020111.1) was identified with a complete and conserved sequence. Baratin (XP_023017557.1) and Ion Transport Peptide A (ITP-A; XP_023016721.1) were recovered as corrected precursor sequences, while Proctolin (GEEF01088300.1) was found to match the reference form without substitutions or receptor loss in Leptinotarsa (Veenstra 2019). Collectively, these neuropeptides represent conserved and verifiable signaling systems encompassing ligands and receptors with diverse and potentially functional roles in the L. decemlineata.

Dissections of L. decemlineata were performed on beetles from a laboratory-reared colony maintained in a growth chamber with 24 °C, 70 % humidity and a photoperiod of 16:8 L:D cycle. All lifestages were fed with non-insecticide treated, fertilized cv. ‘Silverton’ potatoes (Solanum tuberosum) reared continuously in a greenhouse. Adult and larval beetles were dissected to isolate the CNS, adipose (fat body), and midgut for tissue-specific expression analysis. These tissues were chosen because they represent major centers of neuroendocrine signaling and metabolic regulation in insects (Toprak 2020). Tissues were dissected in PBS under a stereoscopic microscope within 15 minutes in order to prevent significant RNA degradation, using micro-scissors, forceps, and insect pins. Tissues were stored in triplicates in 0.2 ml microcentrifuge tubes, snap-frozen in liquid nitrogen, and stored at -80 °C. Whole organism samples were also snap-frozen and stored at -80 °C, including eggs (frozen in groups of 10), first instar larvae (L1, stored in triplicates), third instar larvae (L3, stored in triplicates), pupae (sexed), and adults (sexed).

In order to investigate neuropeptide transcript expression, total RNA was extracted from whole body larvae and dissected tissues in biological triplicates using TRIzol reagent and a spin-column based purification with a PureLink RNA Mini Kit (Thermo Fisher). In the process of purification, RNA was treated with PureLink on column DNase I (Thermo Fisher). All RNA samples were stored at -80 °C, and quantified for concentration and purity using a Nanodrop. 

Single stranded cDNA for qPCR application was synthesized from total RNA using a SuperScript™ III First-Strand Synthesis kit (Thermo Fisher). A concentration of 800 ng was input into each reaction. Synthesized cDNA was stored at -20 °C until further downstream analysis. 

Gene expression of the six focal genes (GPA2, LGR1, Tachykinin, Baratin, ITP-A, and Proctolin) was quantified across developmental stages (first instar, third instar, pupae [sexed], and adults [sexed]; all feeding) using quantitative real-time PCR (qPCR). Gene expression data were relativized by the expression of three housekeeping genes (RP4 and Act-b; Primer list for housekeeping genes found in Supplementary Table 1 of the associated manuscript) and analyzed using the Pfaffl method (Pfaffl, 2001).

Primers for each neuropeptide (Supplementary Table 2) were designed using Integrated DNA Technologies (IDT). Prior to qPCR assays, conventional PCRs were performed for each primer pair to optimize annealing temperature and cycling times. qPCR reactions were run using GoTaq® qPCR Mastermix (Promega) with SYBR® Green I fluorescent dye in a final volume of 10 μL, consisting of 5 μL Mastermix, 0.2 μL of 10 mM forward primer, 0.2 μL of 10 mM reverse primer, 1 μL of cDNA template, and 3.6 μL nuclease-free water. The cycling program included an initial 10 min incubation at 95 °C to activate the polymerase, followed by 40 cycles of denaturation at 95 °C for 10 s, annealing at a primer-specific temperature for 15 s, and extension at 72 °C for 20 s. A melt curve was then generated from 60–95 °C in 1 °C increments every 5 s. All qPCR assays were performed in 0.2 μL tubes using a Qiagen Rotor-Gene-Q instrument.

Standard curves were included on each plate at dilutions of 1, 0.5, 0.2, 0.1, 0.01, 0.001, and 0.0001, along with no-template controls (NTCs). Thresholds were automatically calculated, and Ct values were collected for downstream analysis (Supplementary Tables 3 & 4). Runs were only included if amplification efficiency fell between 0.9 and 1.1, and NTCs did not cross threshold.

To validate target gene expression and confirm tissue specificity prior to qPCR analysis, we compared our candidate neuropeptides to publicly available transcriptomic data from the CPB Gene Expression Atlas (Wilhelm & Xu 2025). All neuropeptide and receptor sequences annotated in Leptinotarsa by Veenstra (2019) were extracted and cross-referenced with expression matrices spanning larval, pupal, and adult stages and corresponding tissues (central nervous system, midgut, fat body, hindgut, and reproductive tissues). Transcripts per million (TPM) values were log-transformed and visualized in a heatmap to compare relative expression profiles across developmental stages and tissue types. Concordant patterns between the atlas data and published beetle neuropeptide identification (Wilhelm & Xu 2025; Veenstra 2019) were used to confirm gene annotation accuracy and to validate tissue-specific expression trends observed in this study (Supplementary Figure 1).

Double-stranded RNA (dsRNA) was generated to conduct gene knockdown trials for GPA2 and LGR1 using the MEGAscript RNAi Kit (Thermo Fisher). We additionally generated a control dsRNA (a dsRNA foreign to the genome of L. decemlineata), using the control template from the MEGAscript RNAi Kit. T7 promoters were added to the 5’ end of the previously described primers and used to generate PCR products of ~200 bp of the neuropeptide amplicon sequence. These PCR products were checked for specificity by gel electrophoresis and purified using ethanol precipitation. The purified PCR products were then converted to dsRNA through in vitro transcription. The purification process involved DNase I and RNase treatments, as well as a spin column based purification. The dsRNA was then quantified by spectrophotometry using a NanoDrop 2000 (Thermo Scientific) and normalized to working concentrations for downstream injection and feeding assays.

Choice of dsRNA Delivery Methods

Direct injection of dsRNA into the hemolymph facilitates systemic knockdown and greater control over the amount of dsRNA applied to each subject, whereas the delivery of dsRNA through feeding minimizes injury-related mortality. In initial trials evaluating dsRNA delivery methods and concentrations, injection treatments resulted in high mortality attributed to puncture wounds. Therefore, we emphasize the results focused on feeding assays, although we provide detail on the injection assays (see Supplementary Methods and Results). 

Feeding of dsRNA

Adult L. decemlineata were collected from the Hancock Agricultural Research Station, Hancock, WI (HARS) and offspring were raised in the laboratory for the dsRNA feeding trial. Third (L3) and fourth (L4) instar larvae were fed dsRNA administered on 6 mm potato leaf disks. Treatment groups were established for both instars and included: FC1 (untreated control), FC3 (500 ng/μl control dsRNA), FGD1 (1000 ng/μl GPA2 dsRNA), FGD2 (3000 ng/μl GPA2 dsRNA), FLD1 (1000 ng/μl LGR1 dsRNA), and FLD2 (3000 ng/μl LGR1 dsRNA). Doses were raised compared to published injection protocols to account for nuclease activity in the gut that has been documented to decrease RNAi efficiency (Wynant et al. 2014). Each treatment was replicated with 20 larvae per instar, for n = 120 total larvae. A total volume of 1μl was applied via micropipette to each disk and air-dried before being fed to individual larvae. All larvae in the dsRNA feeding trial were allowed to feed for two hours on individual disks, and if the entire disk was consumed, larvae were transferred to 150 mm petri dishes with dampened filter paper in groups of 10 with fresh, unaffected potato leaves for the duration of the experiment. Gene knockdown for both candidates in both instars was confirmed via qPCR and gel electrophoresis. 

Phenotypic Assessments

The gene knockdown trials were conducted over 28 days to allow development from larval stages through pupation into adulthood. All L. decemlineata were maintained in a growth chamber with 24 °C, 70 % humidity and a photoperiod of 16 h light and 8 h dark per day. Phenotypic data was collected daily depending on life stage. Across the entire experiment, larvae were scored as alive, dead, or incapacitated, with incapacitated being defined as the inability to right oneself after being flipped dorsally. Lifestage and time to development were also documented with each molt. For the experimental trial, larvae were fed dsRNA as both third and fourth instars, respectively, marked as pupating once fully molted atop shallow soil, and adult when pupal casing was molted. Larval weight (g) and defoliation (%) of potato leaves in a petri dish were taken until all larvae in the replicate pupated or died. Once larvae began to show signs of pupation, they were moved into shallow dishes with soil to replicate a more biologically accurate environment. Once pupation was complete, emerged adults were allowed to feed and reproduce within treatment.

 In order to test the effect of the dsRNA neuropeptide treatments on typical larval development, phenotypic assessments were statistically analyzed separately in R v 4.5.1.

Larval Weight After Knockdown

Average larval weight per dish was analyzed with linear mixed models including treatment replicate (CPB_Group) as a random intercept/repeated factor. Models were fit by maximum likelihood to perform likelihood ratio testing (Faraway 2016). Significant differences between treatment groups were determined using Likelihood Ratio Tests (LRTs), where the full model was compared against reduced models excluding the interaction, treatment, or time. All terms were retained, as their removal significantly reduced model fit.

For overall model form, alternatives were compared using Bayesian Information Criterion (BIC) and residual diagnostics. A generalized additive model (GAM) provided the best fit and was used to visualize the data as average larval weight over time. 

Percent Defoliation After Knockdown

Percent defoliation was analyzed using beta regression with a logit link to account for proportional data. Observations were scaled to fall within 0-1 by applying a constant offset to raw percentages (y = (Defoliation + 0.001) / 100.002​). The model included Treatment, Assay Day, their interaction, and the number of surviving larvae (Alive) as fixed effects. Model precision (ϕ) was estimated with an identity link. Wald z tests were used for inference, and model fit was assessed with pseudo-R², residual diagnostics, and BIC comparisons to alternative linear and polynomial specifications. Predicted values on the response scale were visualized with treatment means ± SE displayed alongside individual dish trajectories.

Survival Rate After Knockdown

Larval survival was analyzed by Kaplan-Meier methods using individual-level records reconstructed from dish counts. Time was assay day (0-28). Treatments were coded as ordered factors, using the uninjected control as the reference. We estimated treatment-specific Kaplan-Meier curves and plotted survival with 95 % confidence bands and risk tables over all assay days. Differences among treatments were tested by the log-rank test (survdiff). Pairwise log-rank comparisons were performed with false-discovery rate (Benjamini-Hochberg) and Bonferroni adjustments for multiple testing.

Time to Development After Knockdown State-Transition Modeling of Beetle Development

For individual beetle development data, interval (start-stop) data conveyed state (lifestage) and the next observed destination state. States were encoded as: 1 = Alive, 2 = Incapacitated, 3 = Pupating, 4 = Adult, 5 = Dead (absorbing). Allowed/observed transitions included 1→2, 1→3, 1→5, 2→5, 3→4, 3→5, 4→5 (transitions with zero events were fit but not interpreted). All time was measured in days over the 28-day assay duration.

For each target transition i→j, there was a transition-specific data set that retained intervals spent in the origin state i. A binary event indicator status was set to 1 if the interval ended with the destination j, and 0 otherwise (including alternative states that were modelled).

Model details

Hazards starting from 1, which is Alive, used a Cox-proportional hazards model (Cox 1972) with a Markov assumption. This means that no prior information about previous states is modeled. Hazards that did not start from 1 relaxed this assumption by including the backward time in the current state, B(t) = stop- start (time since entry to i). This is called a semi-Markov assumption.

The general transition-specific Cox proportional hazards model is:

ijt | X, B(t) = 0,ijt exp{ij(trt)X + ij(B)B(t) + ij(B2)B(t)2},

Here, X is Treatment, which is a factor. To accommodate the difference between models starting from state 1 versus ones not starting from state 1, there are additional terms involving B(t). When starting from 1, both of the leading beta coefficients are set to 0. Only for the Pupating→Adult (3→4) transition, a quadratic term B2 was enabled in the analysis script. All other non-origin-1 transitions set the leading beta coefficient in front of the B(t) term to 0. 

Models were fit using the coxph function in R (Cox 1972; Therneau 2024; Terry & Gramsbach 2000).

Inference and Contrasts

Inference focused on treatment effects within each transition and the dependence of hazards on time in state B(t). Estimated marginal means and pairwise contrasts for treatment were fitted from Cox models using emmeans on the model’s linear predictor to control for multiple comparisons within a model. Each transition was analyzed in its own risk set (beetles “at risk” were those currently in origin state i).