Identifying optimal reference genes for gene expression studies in Eurasian spruce bark beetle, Ips typographus (Coleoptera: Curculionidae: Scolytinae)

Eurasian spruce bark beetle (Ips typographus [L.]) causes substantial damage to spruce forests worldwide. Undoubtedly, more aggressive measures are necessary to restrict the enduring loss. Finishing genome sequencing is a landmark achievement for deploying molecular techniques (i.e., RNA interference) to manage this pest. Gene expression studies assist in understanding insect physiology and deployment of molecular approaches for pest management. RT-qPCR is a valuable technique for such studies. However, accuracy and reliability depend on suitable reference genes. With the genome sequence available and the growing requirement of molecular tools for aggressive forest pest management, it is crucial to find suitable reference genes in Ips typographus under different experimental conditions. Hence, we evaluated the stability of twelve candidate reference genes under diverse experimental conditions such as biotic (developmental, sex and tissues) and abiotic factors (i.e., temperature and juvenile hormone treatment) to identify the reference genes. Our results revealed that ribosomal protein 3a (RPS3-a) was the best reference gene across all the experimental conditions, with minor exceptions. However, the stability of the reference gene can differ based on experiments. Nevertheless, present study provides a comprehensive list of reference genes under different experimental conditions for Ips typographus and contributes to “future genomic and functional genomic research”.

www.nature.com/scientificreports/ and reproducibility compared with other traditional methods molecular techniques [18][19][20][21][22][23] . It is especially useful to detect low-abundance mRNAs in limited samples 24,25 . Nevertheless, RT-qPCR data are influenced by many factors such as initial RNA sample quantity and quality, the efficiency of cDNA synthesis, mRNA recovery, primer and PCR efficiency [26][27][28][29][30][31] . Furthermore, the reliability of RT-qPCR data is highly dependent on the appropriate reference genes as internal controls from the same samples across various biotic and abiotic stresses and treatments 25 . A literature search indicates that reference genes that are constitutively expressed under different environmental factors and maintain the essential cellular functions have been used extensively as internal controls for expression normalization of target genes 32,33 . It is pretty clear now that a single reference gene is not appropriate for the more comprehensive experimental conditions, and it can generate an error in the gene expression estimations causing nonoptimal interpretation of the data 25,34,35 . For mitigation, it is often recommended to use multiple reference genes to minimize variations by RT-qPCR normalization 23,36,37 . Several research studies revealed that most reference gene expression depends on samples/experimental conditions, suggesting no universal reference gene is available for all experimental conditions 38 . Hence, for accurate gene expression normalization, it is essential to evaluate the stability of reference genes for different environmental conditions, life stages, sex-specific and tissue-specific stages for each insect 25,31,32,[39][40][41][42][43] .
In the present study, we analyzed the expression level of 12 commonly occurring reference genes in different coleopteran insects based on published articles 38,43,44 . These genes mostly perform conserved cellular functions in the coleopteran insects, hence expected to be constitutively and stably expressed in all tissues and cells under different experimental conditions 38 . Our objective is to find suitable reference genes for future gene expressions studies in Ips typographus. Candidate reference genes were evaluated across different experimental conditions are elongation factor 1α (EF-1α), ribosomal protein L13a (RPL13a3), arginine kinase isoform X1(ArgK), ribosomal protein L7 (RPL7), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), β-actin (Actin), ribosomal protein 3a (RPS3-a), Tubulin beta-1 chain (β-Tubulin), ubiquitin C variant (UbiQ), V-type proton ATPase catalytic subunit A (V-ATPase-A), ribosomal protein S7 (RPS7), ribosomal protein L6 (RPL6). Reference genes were tested across developmental stages (first, second and third instar larvae, pupa, adult male and female) and tissues (head, midgut, and fat body). Additionally, adult beetles were exposed to juvenile hormone treatment (JHIII) and varied temperatures to identify reference genes for those conditions. Wild and lab-reared beetles were also evaluated for the same intent. Nevertheless, the present study tried to identify reference genes suitable for various tested experimental conditions for I. typographus.

Material and methods
Beetles. Ips typographus were obtained from Kostelec nad Cernými lesy (50° 00′ 07.2″ N 14° 50′ 56.3″ E, under School Forest Enterprise) in summer 2021. School Forest Enterprise in Kostelec nad Černými lesy (Czech University of Life Sciences Prague) is located in the Central Bohemia region, 40 km southwest from Prague. It is a relatively warmer and dryer area, mean annual temperature of 7-7.5 °C, the mean annual sum of precipitation 600 mm, length of vegetation season 150 days. Range of elevation 350-520 m a.s.l. Natural vegetation was composed mainly of beech and oak mixed forest with fir. However, the current species composition is dominated by spruce (50%) and pine (18%) 45 . Due to the extreme drought in 2018, the whole area was under tremendous pressure from the bark beetle outbreak, mainly I. typographus 46 . Collected beetles from the forest were maintained in an insect rearing chamber with fresh Norway spruce logs (Picea abies) at 27 ± 1 °C under 70 ± 5% humidity and a 16:8-h light/dark (L:D) photoperiod at Faculty of Forestry and Wood Sciences, Czech University of Life Sciences, Prague. The wild beetle population was supplied with fresh spruce logs for the next generation and maintained for 35-40 days.
Experimental conditions. Biotic factors (life stages and tissue types). Developmental stages such as three larval stages, pupae, callow male and female, and fed adult male and female were collected from I. typographus. Tissues, including head, fat body, and gut, were dissected from the callow and fed adult males and females (Fig. 1). The samples details are summarized in Table 1. Each replicate was derived by pooling tissues from ten beetles for tissue-specific comparisons. Four biological replicates were used for each bark beetle developmental stages and tissues.

Abiotic factors (temperature, JHIII treatment and wild beetles vs lab-reared beetles).
To examine the effect of temperature treatment, the freshly emerged adults were exposed to 4, 27, and 37 °C for 72 h. For the juvenile hormone treatment (JHIII), fed adult male and female beetles were treated topically on the ventral surface of the abdomen with 10 µg JHIII (Sigma-Aldrich, St. Louis, MO, USA) and acetone as control 47,48 . Beetles were maintained in similar conditions as mentioned above for 72 h. Wild beetles (adults) were collected from the infested Norway spruce stands from Kostelec nad Cernými lesy, Czech Republic, to examine the climatic effect. Similarly, laboratory-reared adults emerged from eggs produced by the wild-beetles. The lab-reared beetles were collected after three generations. Four independent biological replications for each experimental condition were used. The experimental design is summarized in Fig. 1. Selection of candidate reference genes for evaluation. The reference genes selected for this study were reference genes in other Coleoptera species and showed stable expression in our in-house I. typographus genome data 44 (Data not shown). The primers for those selected genes were designed via IDT (Integrated DNA Technologies), and the primer efficiencies (E) and correlation coefficients (R2) were also calculated ( Table 2).
Total RNA extraction and RT-qPCR analysis. Total RNA was extracted from each tissue using TRI-zol® (Invitrogen, Carlsbad, CA) following the manufacturer's protocol. Isolated RNA was further treated with Data analyses. Four different algorithms, namely geNorm 25 , NormFinder 49 , BestKeeper 39 and ΔCT method 50 , were used for measuring reference gene expression stability. geNorm computes the expression stability value (M) and pairwise variation comparison. NormFinder ranks the set of reference genes based on their expression stability within a given set of samples. BestKeeper is also a freely available algorithm that considers the Cq (quantification cycle) value of all reference genes to evaluate standard deviation and correlation coef-  www.nature.com/scientificreports/ ficient. ΔCT method directly evaluates the relative expression of 'gene pairs' within each sample. The mean Cq values of each reference gene from each experiment are offered as input data and subsequently processed using the web-based tool RefFinder (https:// www. heart cure. com. au/ reffi nder/), which delivers a comprehensive stability index that ultimately ranks each reference gene. Pairwise variation (V), estimated by geNorm, was used to determine the optimal number of reference genes for precise RT-qPCR normalization. The Vn/Vn + 1 value exhibited the pairwise variation between two sequential normalization factors 25 .
Validation of selected reference genes. In insects, cytochrome P450 monooxygenases (P450s) are key enzymes that detoxify a broad spectrum of xenobiotics such as plant allelochemicals and synthetic insecticides.
To validate the selected reference gene, we analysed relative expression levels of cytochrome P450 (CYP03903; www.nature.com/scientificreports/ Table 2) in life stages and various tissue types of the beetle according to the 2 − ΔΔCT method for gene normalization using most and least stable reference genes 51 . Target genes mRNA expression level was analysed using one-way ANOVA using GraphPad Prism software. P-value < 0.05 was deployed as the cut-off to indicate significant differences between tested samples.

Results
Candidate reference gene selection and PCR efficiency. Total twelve candidate reference genes, EF-1α, RPL13a, ArgK, RPL7, GAPDH,Actin, RPS3-a, β-Tubulin, UbiQ, V-ATPase-A, RPS7, and RPL6 were selected for identifying suitable reference genes from I. typographus. Each reference gene was produced a single amplicon, as deducted by agarose gel electrophoresis ( Figure S1) and melting curve analyses ( Figure S2). The amplification efficiency of each primer pair fluctuated from 92.2 to 119.24%, and the correlation coefficients (R 2 ) were larger than 0.94 ( Table 2). The Ct values of the twelve candidate reference genes ranged from 19.87 to 34.78 and covered all experimental conditions (Fig. 2). While most Ct values ranged from 19 to 27, Actin, eEF2, β-Tubulin, and RPS3 were the most abundant transcripts under almost all experimental conditions. The least frequently expressed reference genes were NADH, RPL17, and HSP83. The five remaining reference genes were expressed at moderate levels.
Stability of candidate reference genes. Four different algorithms (geNorm, NormFinder, BestKeeper, and delta-CT) were used to find stable reference genes. Gene expression stability was assessed under different experimental conditions such as biotic factors (Table 1; "Section Biotic factors (life stages and tissue types)") and abiotic factors ("Section Abiotic factors (temperature, JHIII treatment and wild beetles vs lab-reared beetles)" using the web-based tool RefFinder.
Biotic conditions. For different developmental stages that include three larval stages (from the first instar to third instar; L1-L3), one pupal stage (P), and two adults (callow and fed), the top three most stable candidates were RPS3-a, EF-1a, and RPS7 based on Ct method, Bestkeeper, RefFinder and Normfinder ( Table 3). The stability ranking order of the first three most stable reference genes that was obtained from four programs, was inconsistent ( Table 3). The geNorm ranking of the top most stable reference genes was RPS3-a, RPL6, RPS7, and RPL13a (Fig. 3A). Integrating the results of all programs identified the consensus top three candidates, RPS3-a, EF-1α, and RPS7, as the most stable reference genes across the developmental stages. Alternatively, GADPH, UbiQ, and RPL7 were the least stable genes (Table 3 and Fig. 3A).  www.nature.com/scientificreports/ For sex-specific tissue comparisons (head, fat body, gut tissues of male and female beetles), the consensus top three candidates were calculated separately for male and female [callow male (CM), female (CF), fed adult male (AMF) and female (AFF)]. The top three most stable genes in male tissues were RPS3-a, EF-1a, and RPS7. The least stable genes defined by four programs were V-ATPase-A, Actin, and ArgK (Table 3). However, based on RefFinder, the most stable reference genes were EF-1α, GADPH, ArgK, and RPL13a, and the least were Actin, V-ATPase-A, and RPL7 (Fig. 3B). In contrast, β-Tubulin, RPS3-a, and RPL6 were the most stable genes in females. The least stable genes were V-ATPase-A, Actin, and RPL7 ( Table 3). The stability ranking for most stable reference genes was constant for females, i.e., RPS7 and RPL6 in the RefFinder (Fig. 3C).
Tissue-specific expression stability of candidate reference genes was calculated for various bark beetle tissues (head, gut, and fat body). In head tissues, the most stable genes were RPL7, RPL6, and RPS3-a (Table 4). RPS3-a and RPL6 were the best combinations of reference genes based on RefFinder (Fig. 4A). Similarly, RPS7, GADPH, ArgK and EF-1α, RPS3-a, RPL6 were the most stable genes within the fat body and gut tissues, respectively (Table 4). In contrast, the best combination of the genes was RPS3-a/RPL6 and V-ATPase-A/RPL6, respectively, as per RefFinder (Fig. 4B,C). The least stable genes defined by four programs were RPL13a, RPL7, and UbiQ in the head, fat body, and gut tissues (Table 4). One interesting unanimous observation was that RPS3-a was the most stable reference gene except for the fat body under biotic conditions. www.nature.com/scientificreports/ Abiotic conditions. Abiotic conditions such as temperature, juvenile hormone III (JHIII) treatment, and laboratory-reared vs wild beetles were considered to evaluate the suitable reference genes using the aforementioned  Table 5). The best combination genes with the lowest expression stability value (M) or highest expression stability for different temperature exposure were RPS3-a and RPS7 (Fig. 5A). The top three stable genes for juvenile hormone III treatment were EF-1a, RPS7, and UbiQ (Table 5). Based on RefFinder comprehensive ranking, Actin, RPL6 and V-ATPase-A (Fig. 5B) were the most stable reference genes. The most stable reference gene expressions between wild and laboratory-reared beetles were EF-1a, RPS3-a, and RPL13a (Table 5), whereas the top three most stable reference genes via Ref-Finder were RPL13a, RPL7 and RPS3-a (Fig. 5C).

Selection of optimal reference genes for normalization.
To determine the consistent and more accurate results of the optimal number of reference genes in each experimental condition, pairwise variation (V) was assessed using geNorm 25 . The optimal number of reference genes by calculating pairing variable value Vn/n + 1 with a cut-off value of 0.15, where the n + 1 reference gene is unnecessary with values under the threshold. Alternatively, the first references gene is sufficient to normalize the target gene expression in those cases (Fig. 6). For biotic conditions minimum of two genes is required for normalization as per V-value above 0.15 in V2/3. in the developmental stages and tissue types were normalized with single reference genes and gene combinations recommended by geNorm (Fig. 6). The relative expression pattern of CYP450 was normalized with the expression level of most stable (RPS3-a, RPS3-a/RPL6) and the least stable reference genes (UbiQ for gut, RPL13a for head tissues). The results showed that CYP450 expression in third instar larva, pupa, and female tissues increased after normalization with the most stable reference gene alone or in combination (Fig. 7A). In compari- www.nature.com/scientificreports/ son, when the least stable UbiQ was used as a reference gene, the CYP450 expression patterns were inconsistent across different stages (Fig. 7A). Similarly, the expression patterns of CYP450 after normalization based on the single (RPS3-a) and combined (RPS3-a/RPL6) genes expression in gut and head tissues were showed a similar pattern (Fig. 7B,C). On the contrary, UbiQ expression-based normalization reduced the expression of CYP450 from 2.9 to 0.2-fold in the gut, whereas RPL13a based normalization increased 4.5-fold in head tissue. Hence, our result demonstrates the importance of selecting and validating reference genes to avoid misinterpretation of gene expression results.

Discussion
RT-qPCR is one of the effective and reliable techniques for quantifying the expression of mRNA levels under different experimental conditions. However, multiple factors, such as RNA extraction, cDNA synthesis, primers, and materials handling, impact the RT-qPCR results. The reliable reference gene can overcome confounding variations in the data. Therefore, the stability of the reference gene must be evaluated for each experimental condition to accurate and reliable data interpretation 37,42,49,52,53 . Ideally, the reference gene should have a steady and unaltered expression level during the entire experimental conditions 25 . Hence, screening the stable reference genes in different experimental conditions before examining the target gene expression is necessary for working on any insect system. www.nature.com/scientificreports/ Although I. typographus is one of the most destructive pests of Norway spruce (Picea abies) 1,3,4,54 , environmental stress such as drought significantly impacts its host colonization. It can provoke transitions from endemic to an epidemic bark beetle population, suitable for the mass attack to healthy trees 55,56 . Unravelling the www.nature.com/scientificreports/ mechanisms underlying I. typographus and host interaction was a daunting task until a decade ago. However, recent advances in genomic and transcriptomic studies on I. typographus 44 have managed to demonstrate the power of gathering valuable genetic information in delineating the beetle-host interaction dynamics. However, it is just the beginning, and without reference genes, there is no way forward for gene expression and RNAi based functional genomics studies. Our study delivers a catalogue of reference genes for the impending genomic and functional studies on I. typographus. In this study, twelve candidates reference genes were tested for the first time in the bark beetle I. typographus across various experimental conditions, including developmental stage, sex-specific, tissues specific, and exposure to abiotic conditions (Tables 3,4,5;Figs. 3,4,5). The results showed that all reference gene primers acquire good amplification efficiency (92.2% to 119.24%) and regression coefficient (0.944 to 0.998). Among the twelve reference genes in this study, we found that RPS3-a, EF-1α, RPL6, RPS7, and RPL7 were the most stable reference genes for I. typographus under various experimental conditions (Fig. 8).
Ribosomal proteins (RPs: RPL-Ribosomal protein large subunit and RPS-Ribosomal protein small subunit) are a large group of proteins that play a crucial role in the cellular process such as protein synthesis, cell growth, and development [57][58][59] Lee et al. 60 demonstrated that ribosomal mutation controls the cellular processes but not direct consequences of ribosome depletion. Several reports indicated that ribosomal protein was widely used as a reference gene in insect functional genomics 31,38,61 . Not surprisingly, ribosomal protein (RPL and RPS) was consistently expressed throughout most of the experimental conditions of the insect species. For example, RPL10 exhibited the most stable expression in different tissues, different diets conditions, and populations of Spodoptera litura, whereas RPS3 was expressed most stably in larvae after starving 62 . RPS13 and RPS7 exhibited the most stable expression under larval-crowding conditions in Mythimna separata 63 . Similarly, different developmental stages (RPL32, RPS18) and different tissues (RPS18) of Spodoptera frugiperda showed the most stable expression 64 , whereas RPS26 and RPL32 genes showed the same in Thermobia domestica 65 ; RPS11, RPL28, and RPL10 genes showed Tuta absoluta 58 ; RPS18 and RPL13 genes showed Rhopalosiphum padi tissues 66 . Sellamuthu et al. 43 reported that RPS3 is stably expressed in the developmental stage, sex-specific and tissue-specific conditions in pine-feeding bark beetle, I. sexdentatus. Furthermore, RPL36 displayed stable expression in the developmental stage, sex-specific and tissue-specific conditions of Glenea cantor 67 . Developmental stage and tissue-specific expression of RPS18 and RPS11 were highly stable in Cicindela campestris 36 . Similarly, RP49 expression was the most stable in tissue and sex-specific conditions in Harmonia axyridis 33 . Wang et al. 68 reported that RPL22e was stably expressed in sex-specific conditions of Mylabris cichorii. RPS23 and RPL12 were identified as reference genes in Anthonomus eugenii 22 . Similarly, RPL32 gene was found as the most stable gene in adult (male and female) tissues of Anoplophora glabripennis 69 . It was reported that RPL13A ribosomal protein was more stable expression under low-temperature treatments in Thitarodes armoricanus 70 and Paederus fuscipes 71 .
Other ribosomal proteins such as RPL27 and RPS7, RPL7 of S. frugiperda, were the most stable reference genes under low-temperature and high-temperature 64 . Tao et al. 72 showed that RPL-33 and RPS-26 are the most stable Figure 6. The optimal number of reference genes for the normalization of I. typographus under different experimental conditions [developmental stages, sex-specific (male and female), tissue-specific (head, midgut, fat body, temperature; Juvenile hormone III; wild-collected vs long-term laboratory-reared beetles]. Based on geNorm analysis, average pairwise variations were calculated between the normalization factors NFn and NFn + 1. Values < 0.15 imply that n + 1 genes were unnecessary to normalise gene expression. www.nature.com/scientificreports/ reference genes from microarray data and RPS-2 and RPS-4 from RNA-seq of Caenorhabditis elegans out of thirteen ribosomal proteins. Moreover, RPS20 was detected as the least stably expressed gene for analyzing Plutella xylostella under different conditions 73 . Our results also demonstrated stable ribosomal gene expression in biotic and abiotic conditions. According to developmental, sex-specific, and tissue-specific stages, RPS3-a, RPL7, and RPS6 showed higher expression stability. Furthermore, RPS3-a and RPS7 were stably expressed as second top rank in various abiotic conditions. These results also suggest that no single reference gene for different experimental conditions is observed in the present study. Elongation factor 1 (EF-1) plays a central role that promotes the delivery of aminoacyl-tRNA to the acceptor site of the ribosome during protein biosynthesis 32,74 . EF-1 has been a stable reference gene for many years and Figure 7. The relative expression levels of the target gene CYP450 were studied under developmental and various tissues after normalizing with the least and most suitable reference genes separately or in combination. The most stable reference gene RPS3-a and RPL6 for developmental stages, head, and gut tissue. The least stable reference gene UbiQ was used for both developmental and gut. For head tissue, RPL13a was used as the least stable reference gene. (A) Developmental stage, (B) Gut. (C) Head. Data represent mean values ± SD of four biological replicates. Asterisks indicate significant differences in the expression of the target gene normalized separately by different reference genes (****P < 0.0001***P < 0.001, **P < 0.01, *P < 0.05, ns indicate no significant difference). www.nature.com/scientificreports/ is widely used for reference genes in many insects [61][62][63]70 . Our results exhibited that elongation factor 1 was the second most stable expression in biotic and abiotic conditions. Similarly, Su et al. 67 reported that EF-1 was the most appropriate reference gene for all samples of Glenea cantor. Teng et al. 75 reported elongation factor2 (eEF2) as the most stably expressed gene in different developmental stages of Plutella xylostella. The expression stability of eEF2 in developmental stages and tissues was also documented in Agrilus planipennis 76 32 reported that the elongation factor was unsuitable for normalization of relative expression due to the different expression variability in different experimental conditions in Drosophila melanogaster. Nevertheless, our results confirmed that EF-1α had sable expression between JHIII treated vs control bark beetles, wild vs lab-reared beetles, and within gut tissues, male tissues, and developmental stages. To further validate our findings and confirm the stable reference gene, we analyzed the expression of cytochrome P450 (CYP03903) in the developmental stages, head, and gut tissues of the adult beetle. P450s is a large class of enzymes that often plays an essential role in detoxifying xenobiotics 79 . It evolves during the insecthost interaction to metabolize a wide range of plant allelochemicals 80 . Our results demonstrated the expression trends in developmental and tissue-specific conditions using single reference gene or gene combinations. The CYP03903 transcript was increased in pupa and female tissue compared to other stages. However, the same expression profile of CYP03903 significantly differs when normalized with the least stable reference gene showing the importance of having a suitable reference gene. For instance, the expression of CYP450 demonstrated just the opposite pattern based on the most or least stable gene was used for normalization. Nevertheless, the CYP03903 expression profile between the gut and head tissues is different, suggesting different functions in the respective tissues 79,81 .
To summarise the present study, twelve candidate reference genes of I. typographus were selected and systematically evaluated for their expression stability using four widely used software programs under different experimental conditions to obtain the best reference genes for each condition. The results showed that the most suitable candidate reference genes were RPS3-a, EF-1a, RPS7, RPL7, and RPL6 under different experimental www.nature.com/scientificreports/ conditions. Based on our comprehensive analysis, we recommended a list of reference genes and combinations of reference genes be used to normalize gene expression in I. typographus subjected to different experimental conditions. This is a much-awaited reference gene validation work on I. typographus, setting the foundation for future molecular (i.e., gene expression) and functional genomics studies (i.e., RNAi). The same genes can be further evaluated to identify suitable reference genes for other Ips species (Coleoptera: Curculionidae: Scolytinae).

Data availability
The raw data supporting the conclusions of this manuscript will be made available by the authors. www.nature.com/scientificreports/