V. vulgaris nests were collected from around the Dunedin/Mosgiel region after advertising free excavation services on social media in collaboration with Nichol’s Garden Centre. Excavated nests were stored in plastic buckets in a 4°C cold room overnight to induce stasis and inhibit function. Individual adults, pupae, larvae, and adults of both sexes were then selected and stored at −80°C for RNA extraction.

Developmental stages included in this study were third-instar larvae, pupae, and three-day-old adults. Developmental stages were identified by specific features: a grub-like appearance for larvae and visualised as the middle size of the larvae present; pupae were selected from capped cells and had clear eye definition; and adults were selected upon emergence from their cells and monitored for 72 hours before freezer storage.

For the sexes, adult workers and adult drones were used. Gynes/queens were excluded from this treatment as they are genetically identical to workers and were not as accessible as workers (e.g., 15 nests with queens would have had to be excavated, as opposed to a single nest for sufficient workers).

The stability of different reference genes in V. vulgaris was assessed across (i) developmental stages and (ii) sexes. Wasps used for life stages were selected the day after nest collection. Larvae and pupae were collected from nest combs, and adults were collected after emerging from their silk caps and stored in foraging boxes for three days then placed in 1.5 mL Eppendorf tubes and stored at −80°C until RNA extraction. Three samples were used for each life stage/experiment, with each experiment conducted in triplicate. Under the sexes condition, three sets of five adult workers and drones were selected separately.

The Qiagen RNeasy Mini Kit was used for RNA extraction. Whole samples were removed from -80°C storage, transferred to a 1.5 ml Eppendorf tube, and placed on dry ice. 300 µl of Trizol was added, and the sample was homogenized using a 1000 µl pipette tip for a minimum of three minutes. More Trizol was added to make up the final volume of 1 ml, along with 200 µl of chloroform. The tube was shaken vigorously by hand for 30 seconds and then incubated at room temperature for three minutes. Samples were centrifuged at maximum speed (10,000 xg) for 15 minutes. 450 µl of the upper aqueous upper phase was transferred to a new sterile 1.5 ml Eppendorf tube with the rest discarded. 1x volume of 100%, RNA-free ethanol was added and mixed by pipetting. Samples were transferred to a RNeasy Mini spin column placed in a 2 ml collection tube, and centrifuged at 12,000 xg for one minute. The flow-through was discarded, and samples were loaded with 350 µl of buffer RW1 and centrifuged at 13,000 xg for one minute. 10 µl DNase 1 stock solution was added to 70 µl of RDD buffer, mixed by gently inverting the tube, and centrifuged briefly. This was transferred directly to the RNeasy column membrane, and incubated at room temperature for 15 minutes. The column was washed with 350 µl RW1 buffer and centrifuged for one minute at 13,000 xg with flow-through discarded. 500 µl buffer RPE was added and centrifuged as above. An additional 500 µl buffer RPE was added and centrifuged for two minutes at 8,000 xg with flow-through discarded. A “dry” centrifuge run at full spin followed to dry out the membrane. The column was transferred to a new sterile Eppendorf, and 40 µl RNase-free water was added directly to the spin column membrane. This was then incubated at room temperature for one minute and centrifuged for one minute at 8,000 xg. This was repeated once, with 5 µl aliquoted for quality assessment, while the remainder of the final samples were stored at −80°C. The 5µl aliquots were assessed for RNA purity, with 1.5 µl used for concentration quantification by measuring the absorbance at a wavelength of 206 nm with a spectrophotometer (NanoDrop 2000C, Thermo Fisher Scientific, USA), with the remainder run on a 1% agarose gel electrophoresis with a 1 Kb+ ladder (Invitrogen) and 2 µl of gel pilot loading dye (Qiagen) to check for clean bands and therefore a lack of degradation.

Samples with a 260/230 ratio of <1.5 underwent a clean-up to improve purity. A 1/10 sample volume of NaAc and 2.5x volume of 100% ethanol was added to each sample and stored at −80°C overnight. Samples were spun at 12,000 g for 30 minutes at 4°C, and the supernatant was removed. 200µl of 70% ethanol was added and samples were spun at 12,000 g for 10 minutes at 4°C. The supernatant was removed, and samples were air-dried for 5–10 minutes. Samples were resuspended in 20 µl RNase-free water and rerun in the Nanodrop to check for altered purity and concentration.

Candidate reference genes were first selected from Polistes genomes available on NCBI. Genes were identified in the V. vulgaris genome (Harrop et al., 2020) with diamond v0.9.24.125 using the command “diamond blastp –query proteins_of_interest.faa –db Vespula_vulgaris”. Proteins were matched to Polistes genome protein models and hits were further classified using visual alignment. Samtools faidx version 1.13-5-gb188dd8 was used to extract the predicted genes from the V. vulgaris genome where the proteins matched. Gene-specific primers were then designed using Primer Wizard (Benchling) and used to clone the open reading frame (ORF) of each reference gene. Primer sequences used for this study are shown in Table 1 .

Reverse transcription of samples into synthesized cDNA was conducted using the following reagents: 4 µl 5x VILO™ Reaction Mix, 2 µl 10x SuperScript^® Enzyme Mix, 1 µg RNA, and nuclease-free water to a total volume of 20 µl. Samples were then incubated at 25°C for 10 minutes, 42°C for 60 minutes, and terminated at 85°C for 5 minutes. Samples were stored at −20°C until ready for further testing. Primers were checked for specificity via a PCR. Two controls were used: a water-only control and 100 ng of genomic DNA from a gyne sample to verify that amplification from genomic DNA results in a bigger fragment. Reagents consisted of 5 µl NEBNext^® Ultra™ II Q5^® Master Mix (Ipswich, MA), 1 µl of each of the forward and reverse primers (10 µM), 1 µl of cDNA or 3 µl of gDNA, and nuclease-free water to a total of 10 µl. Conditions of the PCR amplification comprised of 3 min at 94°C, 35 cycles of 94°C for 10 s, 58°C for 30 s, 72°C for 20 s, and a final 7 min at 72°C. Samples were visualised on a 1% agarose gel electrophoresis with a 1 Kb+ ladder (Invitrogen) and 2 µl of gel pilot loading dye (Qiagen).

To measure the efficiency of the reference genes I performed a PCR on a Light Cycler 480 thermocycler (LC480-II; Roche Diagnostics, Rotkreuz, Switzerland). An RNA sample underwent three reverse transcription reactions at 20 µl each to a pooled total volume of 60 µl. This was diluted with 240 µl of nuclease-free water to act as the undiluted initial sample. A serial dilution of cDNA (1:1, 1:10, 1:100, 1:1,000 and 1:10,000) was made with 5 µl of dilutions for each well and one water control. The qPCR reactions comprised of a master mix with 2x SYBR green mix (10 µl), 10 pmol/µl each forward and reverse primer (0.6 µl), and 3.8 µl RNase-free water to a final 15 µl volume. Each dilution was done in triplicate for each gene. he reaction cycle was as follows: 3 min at 94 °C, followed by 35 cycles of 30 s at 94°C, 30 s at 58°C, and 20 s at 72°C, and a final 7 min at 72°C. Each plate included a water-negative control for each primer pair. Analysis of the amplification plot and melting curve followed each reaction.

All 10 primers underwent an initial RT-qPCR to test their amplification potential. This list was reduced to the six most effective primers, as four failed to produce stable curves across replicates.

Six reference genes with the most stable standard curves were selected for further analyses. The relative expression profiles of these genes were determined at three different developmental stages (larvae, pupae, and adults) and sexes (workers and drones). The stability of these six genes was assessed using Microsoft Excel-based software tools: BestKeeper, geNorm, NormFinder, and the ΔCt method. BestKeeper uses geometric means of Ct values to express the stability of reference genes. A standard deviation (SD) and stability value are used (SV); the lower SV represents more stable genes (21). geNorm automatically calculates the stability measure M by calculating the average pairwise variation of the reference gene against all other genes used in the analysis (16). NormFinder implements an ANOVA-derived model that calculates intra- and inter-group variation and then ranks the reference genes by an SV (22). The ΔCt method, where “pairs of genes” are compared and a rank order is determined based on mean delta Ct scores. In these analyses, genes with the lowest values are ranked as the most stably expressed for that given experimental condition. And finally, a comprehensive ranking of the stability of all candidate reference genes established from these four different statistical algorithms is employed using a Web-based analysis tool, RefFinder (https://www.heartcure.com.au/reffinder/; 22).

The primer specificity was validated by using 1% agarose gels ( Figure 1A ). Primer specificity of all candidate reference genes was further tested via the amplification and dissociation curve analysis of every gene, visualized with a single peak and no detectable signal in the negative control ( Figure 1B ).

Standard curves were generated for every gene based on a five-fold serial dilution of the pooled cDNA ( Figure 1B ). The correlation coefficient (R2) values ranged from 0.8–1, and the PCR efficiency values determined by the standard curve ranged from 81.41%–268.42% ( Table 1 ). From this data, the six following genes were selected for further analyses: CAPZB, KTUB, RPL18X3, AK1, EF1A, and TBP.

Expression levels were defined as the number of amplification cycles required to reach a fixed threshold in the exponential phase of the PCR reaction (26). Thresholds were automatically set by the instrument software in all genes to determine the Ct values. These expression levels were tested across two treatments in V. vulgaris, with variable Ct values highlighting a range of expression levels and varied expression patterns for these six reference genes ( Figure 2 ).

The expression levels of CAPZB, RPL18X3, and TBP showed greater variability within the developmental stage, while the Ct values for AK, KTUB, and TBP were the most varied in the sexes. EF1A and KTUB showed the narrowest range of variability between factors, suggesting a relatively stable expression at developmental stages and sexes ( Figure 2 ).

All candidate reference genes were ranked based on the stability of gene expression for both treatments via four statistical algorithms.

BestKeeper ranked the candidate reference genes AK1 as the most stably expressed candidate reference gene in both developmental stage and sex, and CAPZB as the least stably expressed candidate reference gene in both conditions. Based on the ΔCt method ranking system, the candidate reference genes TBP and CAPZB were the most stably expressed genes for developmental stage and sex, respectively. Analysis of candidate reference genes ranked EF1A as the most stable gene based on both developmental stage and sex expression profiling. Using the analysis tool NormFinder, genes TBP and KTUB were ranked as the most stable genes across developmental stages and sexes, respectively. The overall values of the stability of the six selected candidate reference genes can be viewed in Table 2 .

After the expression stabilities for developmental stage and sex were analysed by the above four methods, RefFinder was employed to calculate an overall stability ranking and therefore the most appropriate reference genes for these experimental conditions.

For developmental stages, analyses using the ΔCt method and NormFinder both ranked the candidate reference gene TBP as the most stable, while BestKeeper and geNorm ranked AK1 and EF1A as the most stable genes, respectively. In sex, the most stable gene varied for each form of analysis, with AK1 for BestKeeper, CAPZB for the ΔCt method, KTUB when ranked by NormFinder, and EF1A when using geNorm ( Table 1 ). These algorithms varied in the ranking of the most unstable genes, putting AK1 as the least stable in sex for all measurements bar BestKeeper, while NormFinder and geNorm consistently ranked AK1 as the least stable gene under the sex condition. RefFinder ranked the stability of these six genes from most to least stable in the developmental stage TBP>EF1A>RPL18X3>CAPZB> AK1>KTUB ( Figure 3A ) and as KTUB>EF1A>CAPZB>TBP>RPL18X3>AK1 under the treatment for sex ( Figure 3B ). geNorm attempted to determine an optimal number of reference genes required to normalise target gene data, but could not recommend an amount as the variability between sequential normalisation factors (based on the n and n+1 least variable reference targets) is relatively high (geNorm V >0.15). Additional candidate reference genes were recommended, or if not possible, to use the four reference targets with the lowest M value when running an experiment under the developmental stage condition, or three reference targets with the lowest M value for experiments concerning sex. This was recommended as using multiple (in this case, non-optimal) reference targets would result in more accurate normalization compared to using a single nonvalidated reference gene.

The findings from this study show a striking difference in the stability of candidate reference genes from one condition to the other in V. vulgaris. While a comprehensive ranking by RefFinder showed that KTUB expressed the highest stability between sexes, it also indicated that the same gene had the most varied stability under the developmental stages condition. TBP and EF1A were considered the most stably expressed gene for the developmental stages condition. The expression patterns of RPB7 across various life stages were inconsistent when normalized with these two most stable reference genes ( Figure 4 ). Transcripts of RPB7 Ct values showed some variation across all developmental stages, while normalization of transcripts using TBP and CAPZB highlighted peak RPB7 expression in pupae ( Figure 4 ). Overall, RPB7 had varied expression patterns, with over 400-fold more expression when normalized with the least stable KTUB gene, compared to a less than four-fold change in expression when normalized with EF1A ( Figure 4 ).