The symbiotic S. oryzae population was sampled in Azergues valley, France in 1984, and has been reared in laboratory conditions ever since. The aposymbiotic strain was obtained by heat treatment in 2010, following a protocol described previously (Nardon, 1973). S. oryzae symbiotic and artificially obtained aposymbiotic weevils were reared in a climate chamber (27.5°C, 70% relative humidity, no light source) on wheat grains as a food source and substrate for egg laying. In these conditions, the timespan between egg laying and the emergence of adults from the grain is 1 month for symbiotic weevils, and 5 weeks for aposymbiotic ones. Adults emerging from the grain are fully formed, although their cuticle is reinforced during the 5 days following emergence from the grain, together with sexual maturation (Maire et al., 2019). The time between the end of metamorphosis and emergence has been calculated to be around 3 days (Vigneron et al., 2014). Since our studies are mainly focused on S. oryzae transcript levels in the gut, from larval stages to the symbiotic clearance (around day 9 of adulthood), we tested candidate reference genes on gut samples (which contains the bacteriomes) harvested from the last larval stage (L4) to day 9 of adulthood (D9), including pupae (P), and intermediate adult stages (D1, D3, D5, and D7).

For sampling, weevils were immersed in buffer TA (35 mM Tris/HC_l, 25 mM KC_l, 10 mM MgCl₂, 250 mM sucrose, pH 7.5) and dissected under a stereomicroscope. Three biological samples per time point were collected, each composed of five guts and their bacteriomes. All samples were collected on ice in 1.5 mL RNAse-free tubes and stored at −80°C until use.

RNAseq datasets were processed in our previous publication (Ferrarini et al., 2023), and raw reads are available at BioProject PRJNA918957. Briefly, quality check was performed with FastQC v0.11.8 (Andrews, 2010), and raw reads were quality trimmed using trim_galore from Cutadapt v0.6.7 (Martin, 2011), then mapped to the S. oryzae’s genome (GCA_002938485.2 Soryzae_2.0) using STAR v2.7.3a (Dobin et al., 2013). Then, uniquely mapping reads were counted with featureCounts v2.0.1 (Liao et al., 2014), and mapping quality was assessed by multiqc v1.13 (Ewels et al., 2016). Differential expression analysis was conducted with SarTOOLS, taking advantage of a global linear model from Deseq2 (Love et al., 2014, 2), and log2 fold changes and adjusted p-values were obtained for sequential pairwise comparisons, i.e., L4 vs. D1, D1 vs. D3 etc., and can be found in the Supplementary material of our previous report (Ferrarini et al., 2023). Most graphic outputs are either performed in R, using ggplot2 (Wickham, 2009) or with GraphPad Prism software.

Primers for each candidate reference gene and GAPDH were designed with Primer3 software (Kõressaar et al., 2018, 3), with “best annealing” at 55°C, amplicon length between 70 and 120 nucleotides (Table 2). Forward and reverse primers were always located in two different exons in order to span at least one intron. In the case of multiple predicted isoforms for each transcript, the primers were designed so that they could amplify all of them. Primer sequences were checked by BLAST (Camacho et al., 2009) to ensure their specificity in the target genome.

Total RNA was purified from samples using the RNAqueous Micro kit (Ambion) and its DNAse treatment, following the manufacturer’s instructions. Final RNA concentration was measured with a Nanodrop^® spectrophotometer (Thermo Scientific), and RNA quality was checked using agarose gel electrophoresis. Reverse transcription to cDNA was carried out using the iScript cDNA Synthesis Kit (Bio-Rad), which relies on a blend of oligo(dT) and random primers. Reverse transcription was performed starting from 1 μg of total RNA. Transcript amplification by qRT-PCR was performed with a CFX Connect Real-Time detection system (Bio-Rad) using the LightCycler Fast Start DNA Master SYBR Green I kit (Roche Diagnostics) as previously described (Maire et al., 2019). Each reaction contained 5 μl of Master Sybr Green mix, 0.5 μl of each primer (10 μM), 1.5 μl RNAse Free water, and 2.5 μl of 1:5 diluted cDNA. Primers for amplification are listed in Table 2. After 5 min at 95°C, the cycling conditions were: 45 cycles at 95°C for 10 s, 56°C for 20 s, and 72°C for 30 s. A melting curve was obtained at the end of each PCR by heating for 30 s at 66°C, then increasing the temperature up to 95°C (increment rates of 0.11°C/s). For every candidate gene, three biological and two technical replicates were tested for each developmental stage, together with a standard curve made with PCR products as templates, at concentrations spanning from 0.2 fg/μl to 2 pg/μl of amplicon, to calculate the reaction efficiency and ensure linearity of the amplification. The data from each transcript were collected within the same plate.

To evaluate the reliability of each candidate reference gene, we compared the average cycle threshold (Ct) of each developmental stage and the percentage covariance (CV%) value for each gene (Boda et al., 2009; Sundaram et al., 2019). We then compared the rankings of various web-based algorithms: GeNorm (Vandesompele et al., 2002), the ΔCt method (Silver et al., 2006), BestKeeper (Pfaffl et al., 2004), and Reffinder (Xie et al., 2012). Each candidate reference gene was compared, ranked and evaluated through the tools described.

To find S. oryzae reference genes to be used across gut and bacteriomes development (from larval stage 4 to 9-day-old adults), we scanned a comprehensive transcriptomic dataset spanning S. oryzae life cycle (Ferrarini et al., 2023). Candidate reference genes should be stable between conditions, highly expressed, and present transcript isoforms suitable for primer design. There are 17 957 genes in S. oryzae (Parisot et al., 2021) and our previous analysis showed that 11 908 were differentially expressed between sequential conditions (Ferrarini et al., 2023). In order to find stable genes across development, we removed genes presenting an adjusted p-value < 0.05 between time points (∼1 500 remaining genes), along with a significant log2 fold change lower than −1 or higher than 1 between two consecutive stages, yielding 100 candidate genes. A selection was made by prioritizing highly expressed genes (mean expression level higher than 1000), in order to improve qRT-PCR detection, decreasing the candidate list to only 19 genes. Among the resulting pre-selected 19 genes, we visually selected those with the most steady expression along the development. We eliminated genes encoding regulators, which are more likely to change in expression according to the experimental conditions. Among the remaining genes, we avoided uncharacterized genes, abundant in the S. oryzae genome. These filtering steps resulted in selecting nine potential reference genes for experimental validation (Table 1 and Figure 1). The selected genes spanned various predicted biological functions, to reduce the risk of co-regulation among them and to facilitate the future use of the same genes in other organisms.

Finally, we also included GAPDH, a very common reference gene in our model and similar model organisms, and the homolog of syntaxin-1 (SYN1A), a protein involved in intracellular vesicle traffic, especially in neurons (Zhou et al., 2000), previously selected as a suitable reference gene by Lord et al. (2010) for Tribolium castaneum, the most widely studied member of Coleoptera. These genes were included in the study as a point of comparison for the newly selected genes. Both genes showed little variation between time points (Figure 1). Altogether, a set of 11 candidate genes were then submitted for further analysis.

While the transcriptomic analysis suggests the candidate genes are stable across S. oryzae’s gut development, it is important to verify that it is indeed possible to amplify such transcripts in an accurate manner by qRT-PCR. A preliminary PCR test allowed amplification of all predicted transcripts from gut cDNA pools, obtaining fragments compatible with the predicted amplicon sizes. Subsequent qRT-PCR quantification from gut samples (larval stage L4, pupal stage P, and adult stages D1, D3, D5, D7, and D9) was successful for nine out of eleven genes, leading to acceptable efficiency and linearity (Table 2). Two exceptions were the primers for KIAA and RAP1, which did not lead to acceptable fragment amplification under standard conditions (efficiency values of 61.6 and 78.4, respectively; average Ct values of 38.01 and 32.58, respectively); and were therefore excluded from the list of potential candidates (Table 2). Moreover, although successful, YTH3 expression levels appear lower than expected based on the transcriptomic analysis. Discrepancies between amplification in qRT-PCR and normalized read counts could be due to low primer efficiency, a tendency to primer dimerization, or suboptimal amplification of transcripts in our standardized qRT-PCR protocol, as well as differences in normalization strategies between the two techniques. One should also note that for qRT-PCR, larval samples consisted of bacteriomes and whole guts, while in the previous transcriptomic dataset on which the candidate reference gene selection was based, only larval bacteriomes were dissected (Ferrarini et al., 2023).

The mean Ct values of amplicons for each candidate reference gene at each developmental stage are reported in Figure 2. The gene with the lowest average Ct was GAPDH (mean for all stages: 20.0), followed by PP12A (mean for all stages: 23.99) and PP43 (mean for all stages: 24.66). All average Ct values are considered low (GAPDH) or medium (all other genes) and are therefore suitable for estimating gene expression of various target genes with different expression profiles. Statistical analysis of Ct values revealed that the mean expression level for the majority of the genes is subjected to mild variations between developmental stages (see asterisks in Figure 2), except for GAPDH and SYN1A.

To identify the most suitable reference gene to compare transcript expression across S. oryzae development, we relied on the comparison of statistical methods recently conducted by Sundaram et al. (2019). After describing the merits and limitations of various comparison algorithms, the authors suggested the following as the best pipeline for the identification of reliable reference genes: (i) visual estimation of expression variation across stages, (ii) analysis of CV% values [(Boda et al., 2009); a cut-off of 50 is recommended]; and (iii) ranking of the remaining candidate reference genes by the NormFinder algorithm (Andersen et al., 2004).

Since the visual estimation of expression variation across stages was already performed and showed mild variation for most genes (Figure 2), we calculated the CV% value for each candidate reference gene. The CV% values of four out of nine genes (YTH3, CYCT1, PAN3, and TFIIF) across all developmental stages are indeed higher than 50 (Figure 3A), reflecting variations observed in particular at the pupal stage (Figure 2). When only adult stages were considered (Figure 3B), only two genes still displayed CV% values higher than 50 (CYCT1 and TFIIF). After excluding genes with CV% higher than 50, we ran the NormFinder algorithm together with other common reference gene ranking algorithms: GeNorm (Vandesompele et al., 2002), the ΔCt method (Silver et al., 2006), BestKeeper (Pfaffl et al., 2004), and Reffinder (Xie et al., 2012). To note, all ranking systems used here are relative to the set of genes that have been chosen and input in the system. Hence, the important information here is the rank that the genes obtain (i.e., their classification), relative to the others. Including the GAPDH and SYN1A genes that are routinely used as “standard” reference genes ensure a point of comparison with the candidate genes we found in our process. Rankings according to all prediction tools are reported in Table 3 for all stages, and in Table 4 for adult stages only. Overall, each of these methods provided good stability values for all the preselected candidate genes, with a consensus for TAO and SYN1A at all stages and for TAO, PP12A, and YTH3 in adult stages only. GAPDH was often ranked last among these set of candidate reference genes.

The pairwise variation value (V) between the ranked genes was calculated with the geNorm algorithm and used to estimate the stability of the normalization factor with the addition of one normalization gene at a time. The cut-off value for pairwise variation used was 0.15 (Vandesompele et al., 2002), meaning that as soon as the value of V drops below 0.15 the number of reference genes was considered suitable for correct standardization of gene expression. The V value for the set of candidate reference genes was always lower than 0.15 (Figure 4). This indicates that two reference genes, which is the minimal number of reference genes accepted by the MIQE guidelines (Bustin et al., 2009), are sufficient for standardization of gene expression across S. oryzae development in the gut, from the L4 larval stages to D9 adults. It also shows that any combination of the proposed nine reference genes is suitable for normalization.

Weevils artificially deprived of symbionts by heat treatment (aposymbiotic) do survive under laboratory standard conditions, although their fitness is partially compromised, e.g., the larval development takes longer, fertility is lower (Nardon, 1973), and aposymbiotic adults are unable to fly (Nardon et al., 1998). In aposymbiotic weevils, the insect cuticle, which is derived from aromatic amino acids mainly supplied by the bacteria, does not fully develop but remains thinner and more delicate (Vigneron et al., 2014). In order to pinpoint symbiosis-specific transcriptional changes, the transcriptomic profiles obtained for symbiotic weevils (Ferrarini et al., 2023) should be compared with the expression profiles obtained at the same stage in aposymbiotic weevils. We took advantage of existing transcriptomic datasets produced from adult symbiotic and aposymbiotic weevils at different adult stages, to compare the expression of all candidate genes between symbiotic stages (Figure 5; Ferrarini et al., 2023). Differential expression analysis of all nine candidate genes between symbiotic and aposymbiotic weevils shows only three genes are equally expressed between symbiotic conditions: PAN3, TAO, and YTH3. Two of these genes were chosen as suitable candidates for symbiotic analysis, TAO and YTH3 and should constitute a good set for symbiotic/aposymbiotic comparisons.