Transcriptome-Based Selection and Validation of Reference Genes for Gene Expression in Goji Fruit Fly (Neoceratitis asiatica Becker) under Developmental Stages and Five Abiotic Stresses

Goji fruit fly, Neoceratitis asiatica, is a major pest on the well-known medicinal plant Lycium barbarum. Dissecting molecular mechanisms of infestation and host selection of N. asiatica will contribute to the determination of best management practices for pest fly control. Gene expression normalization by Real-time quantitative PCR (qPCR) requires the selection and validation of appropriate reference genes (RGs). Hence, 15 candidate RGs were selected from transcriptome data of N. asiatica. Their expression stability was evaluated with five algorithms (∆Ct, Normfinder, GeNorm, BestKeeper, and RefFinder) for sample types differing in the developmental stage, sex, tissue type, and in response to five different abiotic stresses. Our results indicated that the RGs β-Actin + GST for sex, RPL32 + EF1α for tissue type, RPS13+ EF1α for developmental stages along with odor stimulation, color induction, and starvation-refeeding stresses, EF1α + GAPDH under insecticide stress, RPS13 + RPS18 under temperature stress, respectively, were selected as the most suitable RGs for qPCR normalization. Overall, RPS18 and EF1α were the two most stable RGs in all conditions, while RPS15 and EF1β were the least stable RGs. The corresponding suitable RGs and one unstable RG were used to normalize a target odorant-binding protein OBP56a gene in male and female antennae, different tissues, and under odor stimulation. The results of OBP56a expression were consistent with transcriptome data. Our study is the first research on the most stable RGs selection in N. asiatica, which will facilitate further studies on the mechanisms of host selection and insecticide resistance in N. asiatica.


Introduction
Lycium barbarum L. (Solanaceae, Lycium) has been used as a traditional Chinese medicinal herb for a period of more than 2000 years [1]. The populations of this herb are mainly distributed in Ningxia Hui Autonomous Region, China ( Figure 1A) [2]. The dried ripe fruits of L. barbarum (also named wolfberry, goji berry, gouqizi) ( Figure 1B) are rich in polysaccharides, carotene, amino acids, vitamins, and trace elements [3]. Studies have proven that consuming goji berries can enhance anti-aging and reduce age-related vision loss, which has been associated with prolonged life [3,4]. In recent years, people's demand for goji berry is growing every year domestically and internationally, which leads to increasing requirements for high yield and quality [5]. However, the yield and quality of goji berries at harvest are affected by the degree of pest insect infestation during multiple developmental growth stages of L. barbarum L. [6][7][8]. A major insect pest of L. barbarum is the goji fruit fly, Neoceratitis asiatica [8]. Infestation of L. barbarum with goji fruit flies can result in fruit-damage rates of up to 80% every year [8]. It is challenging to effectively the goji fruit fly, Neoceratitis asiatica [8]. Infestation of L. barbarum with goji fruit flies can result in fruit-damage rates of up to 80% every year [8]. It is challenging to effectively control this pest in the field because of its hidden behavior and inadequate knowledge of it ( Figure 1C,D) [8]. At present, most reports on N. asiatica are limited to the studies of biological characteristics and pest control strategies [8,9]. However, the mechanistic understanding of plant infestation and host plant selection by N. asiatica is still unclear. Further research on the understanding of the molecular mechanisms behind the behavior of N. asiatica will facilitate more effective management of goji fruit flies in the field. The study of gene expression is a key step in functional analysis and molecular biology research [10]. Real-time quantitative PCR (qPCR) is a commonly used tool to measure the relative mRNA expression levels of target genes [10][11][12]. The qPCR method has been widely used in scientific research, including insect classification and evolution, insecticide resistance, and insect-host-plant interaction [13][14][15]. When qPCR data are analyzed by the 2 −∆∆Ct method [16], the most suitable reference genes (RGs) must be selected for normalizing target gene expression to improve both the accuracy and reliability of quantitative results [10][11][12]17]. It is generally believed that ideal RGs should be stably expressed in uniform levels under various experimental conditions. In fact, however, studies indicated that there are no RGs that can be stably expressed under all conditions [12,17]. For example, actin commonly used as an RG in insects does not exhibit equally stable expression for the two closely related species Helicoverpa armigera and Helicoverpa assulta under the same experimental conditions [18,19]. In addition, the most stable RGs for the same insect The study of gene expression is a key step in functional analysis and molecular biology research [10]. Real-time quantitative PCR (qPCR) is a commonly used tool to measure the relative mRNA expression levels of target genes [10][11][12]. The qPCR method has been widely used in scientific research, including insect classification and evolution, insecticide resistance, and insect-host-plant interaction [13][14][15]. When qPCR data are analyzed by the 2 −∆∆Ct method [16], the most suitable reference genes (RGs) must be selected for normalizing target gene expression to improve both the accuracy and reliability of quantitative results [10][11][12]17]. It is generally believed that ideal RGs should be stably expressed in uniform levels under various experimental conditions. In fact, however, studies indicated that there are no RGs that can be stably expressed under all conditions [12,17]. For example, actin commonly used as an RG in insects does not exhibit equally stable expression for the two closely related species Helicoverpa armigera and Helicoverpa assulta under the same experimental conditions [18,19]. In addition, the most stable RGs for the same insect species (such as Bradysia odoriphaga) are different under distinct insecticide stresses. For example, Elongation factor 1α (EF1α), Actin (actin), and Ribosomal protein L10 (RPL10) were stably expressed after imidacloprid, chlorfluazuron, and phoxim treatments, respectively [20]. Therefore, the most suitable RGs for each experimental condition should be screened to make sure the accuracy of target gene expression in qPCR analysis. Substantial research has been conducted to evaluate and verify reference genes from different insect species [17,20,21]. The most appropriate RGs have been studied and selected in crop insect pests such as locusts, cotton bollworms, and rice stem borers for mechanistic studies of feeding, oviposition, reproduction, and avoidance of natural enemies [22][23][24]. However, relatively few studies have been conducted on suitable RGs screening for medicinal plant insects [25]. Up to now, there is no report on internal RGs screening related to N. asiatica.
In this study, 15 candidate RGs widely used for qPCR normalization in insects [17,20] were identified from the transcriptomic data of N. asiatica adults. Then the expression patterns of these genes were measured by qPCR at multiple developmental stages, for both sexes, for different tissues, and under five different abiotic stresses (odor stimulation, color induction, insecticide treatment, starvation-refeeding, different temperature). The expression stability of those candidate RGs was evaluated using five different algorithms (∆Ct, NormFinder, BestKeeper, GeNorm, and RefFinder). The most stable and optimal combination of RGs was determined by GeNorm and RefFinder analyses. A target gene, odorant-binding protein 56a (OBP56a) involved in odor recognition [26], was used to further verify the most stable or unstable RGs in both sex's antennae and different tissues and under odor stimulation by qPCR analysis. The proposed research provides the most stable reference genes of N. asiatica in different experimental conditions for further studies on the mechanisms of host selection and insecticide resistance in N. asiatica.

Verification of Primer Specificity and PCR Amplification Efficiencies
Fifteen candidate RGs (α-TUB, β-Actin, EF1α, EF1β, GAPDH, G6PDH, UBC, AK, GST, SDHA, TBP, RPL32, RPS13, RPS15, and RPS18) and one target gene (OBP56a) of N. asiatica were selected for qPCR normalization (Table 1). Sequence details and best BLAST hits for these candidate genes are presented in Table S1. The specific amplification of each primer pair of candidate RGs (Table 1) was confirmed with qPCR. The specificity of each RG primer was validated by melting curve analysis. As shown in Figure 2, the melting curve of all primer sets exhibited a single amplification peak ( Figure 2). The size of amplicons ranged from 92 to 198 bp. The amplification efficiencies (E) for these genes varied from 95.1% for EF1β to 104.2% for RPL32, and the correlation coefficients (R 2 ) varied from 0.994 (RPL32) to 1.000 (TBP) ( Figure S1 and Table 1).  Figure 3 and Table S2.
The Ct values displayed a wide range from 14.51 (AK) to 31.38 (RPS15) in all samples, and

Expression Profiles of Candidate Reference Genes in N. asiatica
The raw cycle threshold (Ct) values of the 15 candidate RGs for qPCR under eight different experimental conditions were collected and are shown in Figure 3 and Table S2.
The Ct values displayed a wide range from 14.51 (AK) to 31.38 (RPS15) in all samples, and the average Ct ranged from 16.51 ± 0.06 (GAPDH) to 25.57 ± 0.25 (TBP), indicating that GAPDH had the highest expression level, whereas TBP had the lowest expression level in the samples tested. The mean Ct values of EF1α (17.12 ± 0.01) in each or total samples with the minimum standard error (SE), indicating that this gene was the most stable of the genes tested under different experimental conditions. Additionally, the most unstable gene was TBP in all samples (25.57 ± 0.25) (Figure 3 and Table S2). The results also showed that EF1α exhibited the smallest variation in expression levels, whereas TBP displayed the largest variation among all samples.

Stability of Candidate Reference Genes under Eight Experimental Conditions
Firstly, the expression stability of 15 candidate RGs was evaluated by four algorithms (∆Ct method [27], GeNorm [28], NormFinder [29], and BestKeeper [30]), respectively. Then the RefFinder program [31] combined the results of the above four algorithms to calculate their geometric mean to obtain a comprehensive reference gene ranking. The optimal number for normalization of qPCR analysis was determined by GeNorm results (Pairwise value Vn/(n + 1) < 0.15) [20]. The optimal combination RGs were selected by GeNorm and RefFinder results [20,31].

Odor Stimulation
Based on the ∆Ct and NormFinder analyses, RPS13 and RPL32 were the most stable genes, while EF1α and RPS18 were by GeNorm and BestKeeper analyses ( Tables 2 and 3). All four algorithms showed that RPS15 and TBP were the least stable genes. GeNorm analysis presented the pairwise value V2/3 as less than 0.15, indicating that the two RGs were the optimal number for the normalization of qPCR data ( Figure 4). Comprehensive ranking order for gene stability using RefFinder was as follows (from the most to least stable): RPS13, EF1α, RPS18, RPL32, GST, GAPDH, β-Actin, SDHA, α-TUB, G6PDH, UBC, AK, EF1β, TBP, RPS15 ( Figure 5 and Table S3). Therefore, a combination of RPS13 and EF1α was the optimal RGs for the normalization of target gene expression under odor stimulation in N. asiatica (Table 4).

Stability of Candidate Reference Genes under Eight Experimental Conditions
Firstly, the expression stability of 15 candidate RGs was evaluated by four algorithms (∆Ct method [27], GeNorm [28], NormFinder [29], and BestKeeper [30]), respectively. Then the RefFinder program [31] combined the results of the above four algorithms to calculate their geometric mean to obtain a comprehensive reference gene ranking. The optimal number for normalization of qPCR analysis was determined by GeNorm results (Pairwise value Vn/(n + 1) < 0.15) [20]. The optimal combination RGs were selected by GeNorm and RefFinder results [20,31].

Odor Stimulation
Based on the ∆Ct and NormFinder analyses, RPS13 and RPL32 were the most stable genes, while EF1α and RPS18 were by GeNorm and BestKeeper analyses (Tables 2 and 3). All four algorithms showed that RPS15 and TBP were the least stable genes. GeNorm analysis presented the pairwise value V2/3 as less than 0.15, indicating that the two RGs were the optimal number for the normalization of qPCR data ( Figure 4). Comprehensive ranking order for gene stability using RefFinder was as follows (from the most to least stable): RPS13, EF1α, RPS18, RPL32, GST, GAPDH, β-Actin, SDHA, α-TUB, G6PDH, UBC, AK, EF1β, TBP, RPS15 ( Figure 5 and Table S3). Therefore, a combination of RPS13 and EF1α was the optimal RGs for the normalization of target gene expression under odor stimulation in N. asiatica (Table 4). EF1α

Color Induction
The results for color induction indicated that EF1α and RPS13 were the most stable genes, on the contrary, TBP and RPS15 were the least stable genes (Tables 2 and 3 Table S3). As shown in Figure 3, all pairwise values were less than 0.15 by GeNorm analysis (Figure 4). Combined with the analysis results of RefFinder and GeNorm, the optimal combination RGs were EF1α and RPS13 for normalization of qPCR target gene expression under color induction in N. asiatica (Figures 4 and 5 and Table 4).

Insecticide Treatment
The results of ∆Ct showed that EF1α and GST were the most stable genes, however, EF1α and GAPDH were with the other three algorithms (GeNorm, NormFinder, and BestKeeper) (Tables 2 and 3). The TBP gene was the most unstable gene with all analysis methods used (Tables 2 and 3). The RefFinder results indicated the rank order for gene stability was as follows (from the most to least stable): EF1α, GAPDH, GST, RPS18, SDHA, RPS13, RPL32, RPS15, α-TUB, AK, β-Actin, G6PDH, EF1β, UBC, TBP ( Figure 5). Among them, EF1α and GAPDH were considered as the optimal combination RGs for qPCR normalization under insecticide (imidacloprid) stress in N. asiatica, as determined from GeNorm results (V2/3 < 0.15) (Figures 4 and 5 and Table 4).     Table 1. The most stable genes are listed on the left, while the least stable genes are listed on the right.   Table 1. The most stable genes are listed on the left, while the least stable genes are listed on the right.

Starvation-Refeeding
The EF1α was the most stable gene in all analysis methods used except for the NormFinder method which showed the most stable gene to be RPS18 (Tables 2 and 3). All analyses showed that RPS15 and EF1β were the least stable genes (Tables 2 and 3). Based on the RefFinder, the comprehensive bank order of gene stability was as follows (from the most to least stable): EF1α, RPS13, RPS18, GAPDH, RPL32, TBP, α-TUB, G6PDH, GST, AK, SDHA, β-Actin, UBC, EF1β, RPS15 ( Figure 5). GeNorm results (V2/3 < 0.15) determined that a combination of two RGs (EF1α and RPS13) was a better choice for qPCR normalization under starvation-refeeding treatment in N. asiatica (Figures 4 and 5 and Table 4).

Different Temperature Treatments
EF1α and RPS18 were the most stable gene in the GeNorm and BestKeeper analyses, whereas RPS13 in ∆Ct analysis and AK in NormFinder analysis (Tables 2 and 3). Four analysis results indicated the least stable genes were TBP and G6PDH (Tables 2 and 3). The RefFinder results revealed that RPS13 (2.63) and RPS18 (2.63) had a lower geometric mean and were the most stable genes ( Figure 5 and Table S3). Based on GeNorm analysis, all pairwise values were less than 0.15, except for V14/15 ( Figure 4). Thus, RPS13 and RPS18 were the optimal combinations for qPCR normalization under different temperature treatments in N. asiatica (Figures 4 and 5 and Table 4).

Different Developmental Stages
GAPDH and RPS13 were the most stable genes in the ∆Ct and NormFinder analysis, while EF1α was the most stable gene in the GeNorm and BestKeeper analyses (Tables 2 and 3). All analyses showed that TBP and G6PDH were the least stable genes (Tables 2 and 3). The bank order of gene stability by RefFinder was from the most to least stable genes as follows: RPS13, EF1α, GAPDH, RPS18, GST, UBC, EF1β, RPL32, β-Actin, SDHA, RPS15, α-TUB, AK, TBP, G6PDH ( Figure 5). In view of the GeNorm results (Vn/(n + 1) < 0.15), RPS13 and EF1α were proposed as the optimal combination for qPCR normalization in developmental stages of N. asiatica (Figures 4 and 5 and Table 4).

Both Sexes
The GeNorm and NormFinder results showed that β-Actin and EF1α were the most stable genes, whereas the other two analyses had the same gene found to be the most stable RPS13 (Tables 2 and 3). TBP was identified as the most unstable gene in all analyses (Tables 2 and 3 and Figure 5). From the results of RefFinder (from the most to least stable genes): β-Actin, GST, RPS13, EF1α, RPL32, SDHA, G6PDH, UBC, GAPDH, EF1β, AK, RPS18, RPS15, α-TUB, TBP) and from GeNorm analysis (Vn/(n + 1) < 0.15), the optimal combination RGs for qPCR normalization in the male and female of N. asiatica adults were β-Actin and GST (Figures 4 and 5 and Table 4).

Different Tissues
The ∆Ct and NormFinder data indicated RPL32 and GST were the most stable genes. However, GeNorm and BestKeeper results showed that EF1α and RPS18 were the most stable genes (Tables 2 and 3). All analyses showed that EF1β was the least stable gene (Tables 2 and 3). The rank order for gene stability determined by RefFinder was as follows (from the most to least stable genes): RPL32, EF1α, GST, RPS18, RPS15, GAPDH, UBC, AK, β-Actin, SDHA, RPS13, TBP, α-TUB, G6PDH, EF1β ( Figure 5). Hence, RPL32 and EF1α were considered as the optimal combination RGs for qPCR normalization in different tissues of N. asiatica, as determined from GeNorm results (V2/3 < 0.15) ( Figure 4 and Table 4).

Comprehensive Ranking Analysis of Candidate Reference Genes
The EF1α was the most stable gene in total samples by GeNorm and BestKeeper analyses, while RPS18 and GST were the most stable genes in ∆Ct and NormFinder analyses, respectively (Tables 2 and 3). The RPS15 and EF1β were the least stable genes in all analyses (Tables 2 and 3). The overall rank order for gene stability determined from RefFinder results was as follows (most to least stable): RPS18, EF1α, GAPDH, RPS13, GST, β-Actin, UBC, RPL32, SDHA, AK, α-TUB, G6PDH, TBP, EF1β, RPS15 ( Figure 5). The GeNorm analysis showed that all pairwise variation values were less than 0.15, except for V13/14 and V14/15 (Figure 4). Based on the RefFinder data, RPS18 and EF1α are the most suitable internal RGs for normalizing target gene expression in N. asiatica (Figures 4 and 5 and Table 4).

Validation of the Selected Candidate Reference Genes
OBPs can recognize and transport odor molecules to facilitate insects in host plants' locations [32]. In this study, OBP56a of N. asiatica was selected as the target gene to verify the expression stability of potential RGs by qPCR at both sexes' antennae, different tissues, and under odor stimulation. The sequence analysis and amplification efficiencies of OBP56a by qPCR were listed in Supplementary Material (Figures S2-S4). According to our results (Table 4), a single most stable RG or optimal combination RGs and the corresponding most unstable RG were provided for normalizing OBP56a expression in N. asiatica (Figures 6 and 7).     As shown in Figure 6, the transcript levels of OBP56a in the antennae (male and female insects) and different tissues (antennae, legs, and ovipositor of female) of N. asiatica adults were firstly shown by FPKM values of transcriptome data ( Figure 6A,B). Further verification by qPCR indicated that target OBP56a had a similar expression level by using either single RG or in combination RGs (β-Actin + GST in sex, RPL32 + EF1α in tissues) to normalize its expression in both sexes and different tissues, respectively. However, the expression profiles of OBP56a differed significantly when using the least stable RG (TBP or EF1β) for normalization (p < 0.05) ( Figure 6C,D). By comparison, it was found that the expression patterns (mainly in female antenna) of OBP56a by using the most stable RGs (β-Actin + GST, RPL32 + EF1α) for normalization were consistent with the expression trend of transcriptome data in male and female antennae and among different tissues (FPKM) (Figure 6C,D).
Compared to the unstimulated control, the expression of OBP56a based on FPKM values were upregulated in the antennae, legs, and ovipositor of female after odor stimulation (host volatiles), especially in the ovipositor ( Figure 7A). Meanwhile, the qPCR results for single stable RG (RPS13 or EF1α) or a combination stable RGs (RPS13 + EF1α) normalization showed a similar upregulation trend after odor stimulation ( Figure 7B-D). Nevertheless, the OBP56a expression was significantly downregulated in the ovipositor under the least stable RG (RPS15) normalization (p < 0.01) ( Figure 7E). After odor stimulation, relative expression levels of OBP56a in antennae, legs, and ovipositor of females had a significant difference between the most stable RGs and the least stable RG normalization (p < 0.05) ( Figure 7F).

Discussion
L. barbarum is a traditional medicinal herb that has been used as a powerful anti-aging agent due to its health benefits [1,4]. Infestation of goji fruit fly N. asiatica not only causes the decline of L. barbarum's production and quality but also restricts the healthy development of the L. barbarum industry [8]. However, the mechanisms underlying the host selection of N. asiatica are still unclear. The method of qPCR is a rapid and accurate method for detecting and analyzing expression levels for genes of interest [11]. RG normalization is a key step needed for enhancing the accuracy of expression level analysis of target genes by qPCR [12]. So far, no study on reference gene screening has been reported in N. asiatica. Therefore, we identified 15 commonly used candidate RGs from the transcriptome data of N. asiatica. Their expression stability was evaluated by different algorithms to clarify the optimal RGs in eight different experimental conditions for further studies of gene function.
To sum up, our study demonstrated that the expression stability of RGs differed under different experimental treatments and no RG has stable expression in all conditions in N. asiatica. This result is in line with previously reported results [10][11][12]. Actin was commonly selected as a stable RG for gene expression analysis in developmental stages, sex, tissues, and other different conditions in insects, especially in lepidopteran species [33,34]. However, the current study indicated that Actin was only expressed stably in both sexes rather than across developmental stages (RPS13) and tissue types (RPL32) of N. asiatica. In addition, studies on two closely related species of N. asiatica showed that three RGs (TUB, GAPDH, GST) for Bactrocera minax and RPL60 for B. cucurbitae were the most stable RGs for qPCR analysis under tissue-development condition, however, actin was shown to be the most unstable RG [35,36]. It is a normal phenomenon that there are differences in the selection of appropriate RGs under the same conditions among diverse insect species [20,35,36]. Sometimes RGs differ greatly in different developmental stages even in the same tissue [20,36]. Two or more RGs, in fact, can improve the accuracy of qPCR results compared to a single RG [20]. In our study, GeNorm data showed that the pairwise variation value at V2/3 below the proposed 0.15 cut-off value, which indicated that the combination of two stable RGs would be a better choice for qPCR normalization, rather than the use of a single RG under eight experimental conditions.
In the present study, the stability of RGs in N. asiatica could differ under various experimental conditions. Our results were similar to previous studies [17,20,21]. RPS18 and EF1α were the most stable RGs under sex, tissues, developmental stages, and five abiotic stresses in N. asiatica. Among them, RPS13 and EF1α were the suitable combination of RGs under developmental stages, odor stimulation, color induction, and starvation-refeeding stresses. β-Actin + GST and RPL32 + EF1α were stably expressed in both sexes and in the different tissues of N. asiatica, respectively. Moreover, EF1α + GAPDH and RPS13 + RPS18 were selected as the optimal combination RGs to normalizing gene expression under insecticide treatment and different temperature stresses, respectively. Previous research has shown that the expression stability of RGs is different under species, growth stage, and biotic and abiotic stresses [17,20,35]. In the larva and adult stages of B. cucurbitae, RPL32 and RPS13 were the most stable RGs under different temperature stress [37]. In addition, α-TUB and RPL32 were found to be suitable RGs for normalizing the qPCR data in Phenacoccus solenopsis under temperature stress [38]. Actin2 and α-TUB were the most appropriate RGs in both males and females of B. dorsalis, while Actin5 + α-TUB and Actin3 + α-TUB were considered the best RG combinations for females and males, respectively [39]. Meanwhile, in Anastatus japonicus, Actin and EF1α were optimal RGs in male and female adults [40].
There are a few differences between our results and previous reports for both sexes and under temperature stress [37,39]. When insects were treated with different insecticides, their reference genes could not always express stably [20,41]. For instance, EF1α was the most stable RG in Bradysia odoriphaga under imidacloprid treatment, while Actin was expressed most stably for chlorpyrifos and chlorfluazuron treatments [20]. Furthermore, a previous study also indicated RPS15 was found to be the most stable reference gene in Br. odoriphaga under different developmental stages and different temperature conditions [41]. This result, however, was the opposite of the finding in our study, that RPS15 was the least stable RG in N. asiatica. These findings further confirmed that screening and verification one of the expression stability for RGs, are very important for gene research under each specific experimental condition.
Randomly selected or unvalidated reference genes are not recommended for direct qPCR normalization of target gene expression, because they cannot ensure accuracy relative to quantification [20,36]. For example, in this study, the expression level of target OBP56a in the ovipositor of female N. asiatica after stimulation of host plant volatiles was significantly up-regulated using EF1α and RPS13 (the most stable genes) as internal control, while was significantly down-regulated using RPS15 (the least stable gene) as a reference gene (Figure 7). Validated stable reference genes for normalizing gene expression could obtain better quantitative results [35][36][37][38][39][40][41]. It has also been documented that at least two or more reference genes should be selected for gene quantitative normalization [10]. In this study, OBP56a mainly involved in odor recognition with a conserved six cysteine residues pattern [26,32] was used for reference gene verification ( Figure S3). In N. asiatica, β-Actin and GST were the best choices for normalizing OBP56a expression in males and females via GeNorm (V2/3 < 0.15) and RefFinder algorithms. qPCR results showed that OBP56a had a higher expression in female antennae than in male antennae, which was consistent not only with the transcriptome data of N. asiatica, but also with a previous study that used RPL18, RPS17, and EF1α as RGs in Anastrepha obliqua [42]. However, using the unstable TBP and RPS15 as RGs failed to detect a significant effect on OBP56a expression under both sexual antennae and odor stimulation. Therefore, optimization of reference genes is crucial for the accurate normalization of gene expression, especially for the identification of the relatively subtle difference. To improve the accuracy of target gene quantification, the appropriate reference genes must be screened and verified under the experimental conditions before qPCR detection.
In conclusion, this research is the first to validate reference genes and detect the expression profiles of OBP56a in both male and female antennae, different tissues, and odor stimulation in N. asiatica. We concluded the corresponding optimal reference genes under eight different experimental conditions based on five algorithms (Table 4). Our findings provide a scale for quantitative gene expression analysis by qPCR under various developmental stages and abiotic stresses in the pest species N. asiatica. The purpose of the research is to lay a foundation for future elucidation of the mechanisms for host selection and insecticide resistance in N. asiatica. Additionally, our results further emphasize that the most suitable reference genes should be screened and validated under different experimental conditions to ensure the accuracy of gene quantification.

Insect Rearing
N. asiatica's larvae and pupae collected from the main production area of wolfberry in Ningxia, China, were maintained at the Institute of Medicinal Plant Development, Chinese Academy of Medical Sciences, and Peking Union Medical College. After emergence, adults were fed with 10% honey water and reared in a cage (50 cm × 50 cm × 50 cm) made of gauze. The female adults after mating were transferred to a new cage with wolfberry young fruits for further oviposition. The conditions were controlled by an artificial climate cabinet (PXZ-430B, Ningbo Jiangnan Instrument Factory, Ningbo, China) with a 14 Light (L):10 Dark (D) photoperiod at 25 ± 2 • C and 40% relative humidity.

Collection of Samples under Different Experimental Conditions
The N. asiatica adults (males and females) from three days after emergence were used in the experiments. Each experiment was set as three biological replicates. The collected methods used in this study were similar to previous studies [20,33]. The following is detailed information on the sample collection under different experimental conditions. (1) Odor stimulation. The 20 N. asiatica adults were exposed to three odor compounds (1 mM methyl salicylate, 1 mM methyleugenol, and 1 mM 4-(p-Acetoxyphenyl)-2-butanone) for 2 h (h) respectively. As a control, the adults were unstimulated with odor compounds. The stimulated and unstimulated adults were collected separately. (2) Color induction. The 20 N. asiatica adults were kept in three colored plastic boxes (20 cm × 10 cm × 10 cm, green, yellow, and red) for 24 h, respectively. The other 20 adults were induced without color as a control group. After 24 h, N. asiatica adults for the corresponding treatments were collected separately. (3) Insecticide treatment. The N. asiatica adults were treated with different concentrations of imidacloprid by spraying method for 24 h. These concentrations included 0.1 mg/L, 0.01 mg/L, and 0.001 mg/L. Twenty individuals were used for each concentration. The control group was 20 N. asiatica adults treated with water for 24 h. Then these adults were collected in different 1.5 mL microcentrifuge tubes. (4) Starvation and refeeding treatment. The 20 N. asiatica adults were starved for 24 h. The other 20 adults starved for 24 h were refed with 10% honey for 24 h. The normally reared 20 adults were used as the control group. After 24 h, these samples were collected separately. (5) Different temperature treatments. The N. asiatica adults were exposed to five different temperatures (including 15 • C, 20 • C, 25 • C, 30 • C, and 35 • C) for 2 h. Twenty individuals were collected for each temperature. (6) Different developmental stages. One hundred eggs (3 days old), 20 3rd instar larvae, 10 pupae (5 days old), and 20 adults (a mixture of males and females, male:female = 1:1) were collected separately. (7) Both sexes. Twenty individuals from each adult female and male N. asiatica were collected. (8) Different tissues. Tissues were dissected from N. asiatica adults such as 300 antennae, 30 heads without antennae, 30 thorax, 30 abdomen, and 180 legs, kept in different microcentrifuge tubes (1.5 mL, no RNase). All samples collected above were quickly frozen in liquid nitrogen and stored at −80 • C.

RNA Extraction and cDNA Synthesis
The technique and methods used in this study were similar to our previous studies [33,34]. The RNA and first-strand cDNA of these above-collected samples were extracted and synthesized using the Invitrogen TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) and the primeScriptTM RT reagent Kit with gDNA Eraser (Takara, Dalian, China), respectively. The total RNA quantity and quality were determined by a NanoDrop 2000 spectrophotometer (NanoDrop, Wilmington, DE, USA). Then 1 µg of total RNA was used to synthesize the cDNA.
The technique and methods of qPCR used in this study were similar to our previous studies [33,34]. Each qPCR reaction was conducted in a 20 µL reaction: 10 µL of SYBR ® Premix Ex TaqTM (TliRNase H Plus) manual (Takara, China), 0.5 µL of each primer (10 µM), 1 µL of sample cDNA, and 8 µL of sterilized ddH2O. The reactions were performed on a StepOne thermocycler (Bio-Rad CFX, USA) using the two-step method (run as follows: 94 • C for 2 min, followed by 40 cycles of 94 • C for 15 s, 60 • C for 30 s, 60 • C for 1 min) and were analyzed with a melting curve analysis program (heated to 95 • C for 30 s and cooled to 60 • C for 15 s). The Ct values were obtained from qPCR results analyzed by the Bio-Rad CFX 3.0 software.

Determining the Expression Stability of Candidate Reference Genes
The qPCR data for each of the eight experimental samples were analyzed independently. Each sample for one RG was set as three biological and three technical replicates by qPCR analysis. Four algorithms (∆Ct method [27], GeNorm [28], NormFinder [29], BestKeeper [30]) were used to evaluate the expression stability of 15 candidate RGs. The comprehensive rank was calculated by RefFinder [31], which could unify and merge the above four algorithms. A suitable number of RGs for qPCR normalization was determined by GeNorm analysis [20,28]. Analyses of these algorithms used in this study were similar to previous studies [10][11][12]20,21].

Validation of the Selected Reference Genes
To confirm the reliability of the potential reference genes, the relative expression of target OBP56a was measured by qPCR at the antennae of both sexes (male and female) and different tissues (antenna, leg, ovipositor) of N. asiatica female adults and under odor stimulation (host plant volatiles). The OBP56a was identified from transcriptome sequencing of N. asiatica adults, whose details (including sequence analysis, and primer evaluation) were listed in Figures S2-S4 and Table S1. The OBP56a expression was normalized with the most stable and least stable RGs obtained from the comprehensive analysis of GeNorm and RefFinder. The study method of the OBP56a by qPCR was described in detail in Section 4.4. Each experimental sample was set as three biological and three technical replicates in qPCR analysis. The relative expression level of target OBP56a was calculated according to the 2 −∆∆Ct method [16]. The FPKM value of OBP56a in transcriptome data was used as a reference for its expression pattern.

Statistical Analysis
All statistic comparison was determined using SPSS 22.0 software (SPSS Inc., Chicago, IL, USA). Data multiple comparisons were assessed by ANOVA following Tukey's HSD multiple comparison test (p < 0.05). The statistical significance of the difference between the two treatments was analyzed using a pairwise Student's t-test (p < 0.05) [33,34].  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data in this study will be available from the corresponding author upon reasonable request.