Reference genes for gene expression studies by RT-qPCR in Brevipalpus yothersi (Acari: Tenuipalpidae), the mite vector of citrus leprosis virus

Quantitative reverse transcription PCR (RT-qPCR) is a high-throughput method to analyze the transcriptional expression of genes. Currently, no reference genes have been described for evaluating gene expression in Brevipalpus yothersi, the false spider mite, a polyphagous that act as vector of the citrus leprosis virus C (CiLV-C), an important citrus disease. This study aimed to identify the most stable reference genes in B. yothersi. The RT-qPCR expression data for selected genes were evaluated from three conditions: different developmental stages, plant hosts and acquisition of CiLV-C. To analyze the stability of the candidate reference genes we used ΔCq method, GeNorm, NormFinder, BestKeeper and RefFinder. Ubiq and GAPDH are best suited for normalizing gene expression data in viruliferous and non-viruliferous mites. Ubiq, EF1α and GAPDH are the most stable for different developmental stages. RPL13 and RPL32 are the best reference genes for approaches to B. yothersi in different host plants. Considering all the experimental conditions, Ubiq, EF1α, and GAPDH were the most stable genes. Here we developed an accurate and comprehensive RT-qPCR strategy for use in B. yothersi gene expression analysis. These results will improve the understanding of the biology of the false spider mites and their role as virus vectors.


B. yothersi from different host plants.
In samples of mites that were reared in different plant hosts, the gene pairs RPL13/ RPL32 and GAPDH/ RPL13 were ranked as the most stable by the ΔCq method and GeNorm, respectively. On the other hand, NormFinder and BestKeeper attributed GAPDH and RPL32 as the genes with highest stability in samples from it same treatment. The least stable genes in false spider mites from different plant hosts were Tub/EF1α according to GeNorm and NormFinder, and Act/Tub using ΔCq and BestKeeper.
Finally, considering all experimental conditions (virus acquisition, development stages, and plant hosts) in global analysis, distinct algorithms resulting on dissimilar rankings to reference genes. Both ΔCq method and BestKeeper identified EF1α as the most stable gene, NormFinder ranked Ubiq as the one with lowest variation between all samples, and GeNorm ranked GAPDH/Act genes as the best between all candidates. In the same way, the least stable genes were RPL13 to ΔCq and GeNorm, EF1α using NormFinder, and Tub according BestKeeper results.
Ranking of best reference genes using RefFinder. The RefFinder tool was used to define the overall final ranking of reference genes (Table 1), by calculating the geometric mean (GM) of the weights attributed to each candidate gene for each software described above 29 . According to RefFinder parameters, the genes Ubiq (GM = 1.68) and RPL13 (GM = 6) were the most and least stable, respectively, when samples from all experimental conditions were analyzed together. In viruliferous and non-viruliferous mites, the stability of expression (from highest to lowest) ranked the genes as follows: www.nature.com/scientificreports www.nature.com/scientificreports/ different development stages of the false spider mites, the gene reference order recommended by RefFinder (starting from the best candidate) was: Ubiq > EF1α > GAPDH > RPL32 > Act > Tub > RPL13. Finally, when mites were reared in different host plants, the genes were ranked from the most to the least stable in the following order: Number of reference genes for normalization by Rt-qpCR. For robust and consistent RT-qPCR results, at least two reference genes are required. Thus, the pairwise variation (V n /V n+1 ) was performed using GeNorm to determine the ideal number of reference genes for each experimental condition (Fig. 2). The cut-off value of 0.15 was used to indicate the need for inclusion of additional reference genes in a normalization factor 19 . In the experimental assays with non-viruliferous and viruliferous mites, or mites from different plant hosts, the values of V2/3 were lower than 0.15, indicating that two genes are reliable to accurate normalization of RT-qPCR. On the other hand, in false spider mites from different developmental stages, or samples from all conditions merged, the V2/3 parameter were 0.179 and 0.186, respectively, indicating that the addition of a third reference gene would improve the reliability of normalization.
Validation of reference genes. The expression of the V-type proton ATPase catalytic subunit A (V-ATPase), a false spider mite target gene, was evaluated to validate reference genes from different developmental stages of B. yothersi. The V-ATPase transcript levels were normalized with the three most stable (Ubiq, EF1α e GAPDH) and the two most unstable (RPL13 e Tub) reference genes (Fig. 3), as determined by RefFinder ( Table 1). The normalization of V-ATPase expression with the best or the three best reference genes resulted in nearly invariable expression between mite stages. However, when the two worst unstable genes were used, individually or combined as a normalization factor, the V-ATPase fold-change increased in larvae, protonymph and deutonymph relative to the egg control (Fig. 3).

Discussion
Gene expression analysis is a powerful way to understand the uncommon reproduction system of Brevipalpus and the complex interaction between these mites with plants and viruses. Nevertheless, information on the expression of false spider mites genes is scarce, particularly due to the absence of standard protocols. RT-qPCR has been a sensitive, specific and reproducible tool for gene expression analysis 30 . To obtain feasible RT-qPCR Figure 2. Determination of the optimal number of reference genes for RT-qPCR normalization. The pairwise variation (V n /V n+1 ) was analyzed for normalization factors NF n and NF n+1 by GeNorm software to determine the optimal number of references genes required for RT-qPCR data normalization. Values < 0.15 indicate that additional genes are not required for the normalization of gene expression. www.nature.com/scientificreports www.nature.com/scientificreports/ measurements, the expression of target genes must be normalized by reference genes with stable expression in the different biological conditions that the target genes will be tested 31,32 .
In the last years, an increasing number of studies has established suitable reference genes for arthropods [33][34][35] , in case of phytophagous mites, the studies of reference genes includes xenobiotic toxicity, development stage, abiotic stress and shift plant host 22,23,36,37 . GAPDH was found as best housekeeping genes for T. urticae in case of plant shift studies and RPL13 was cited one as the best reference genes for development studies for this mite 36 . While for Panonychus citri, in development studies and abiotic stress, besides GADPH, EF1α was the most stable genes 23 . RPL13 was chosen as best in development stages studies a for T. cinnabarinus 22 . Whereas for Brevipalpus species there are no reference gene information.
Here, we evaluated the expression stability of seven B. yothersi candidate reference genes under different conditions: CiLV-C acquisition (viruliferous and non-viruliferous populations), different developmental stages (eggs, larvae, protonymphs, deutonymphs and adult females), and rearing on distinct host plants (sweet orange and common bean). Our study show that all seven candidates are stable under the tested conditions and can be immediately used for normalization of RT-qPCR studies, because all tested genes exhibited high expression stability with M < 0.86. We also present the best ones for each biological condition and the recommended number of reference genes to obtain reliable results. We used four different algorithms (Comparative ΔCq, NormFinder, GeNorm and BestKeeper) to identify the most stable genes. Despite small differences in their results, as expected 38,39 , similar sets of reference genes were identified for each experimental condition. The RefFinder tool was used to combine the results from the distinct algorithms in a final classification.
The RefFinder classified the Ubiq and GAPDH genes as the most stable across viruliferous and non-viruliferous false spider mite populations ( Table 1). The main function of ubiquitin proteins is the selective degradation of short-lived protein in eukaryotic cells, performing a housekeeping role in the control in numerous cellular process 40 . GAPDH encodes an enzyme essential for glycolysis and glycogenesis pathways, among others functions in several cellular processes 41,42 . Ubiq gene was previously identified as the most stable candidate gene in Delphacodes kuscheli infected with Mal de Río Cuarto virus (MRCV) 43 and Bombus terrestres infected with Israeli acute paralysis virus 44 . At the same viruliferous and non-viruliferous comparison, we identified Act and RPL32 as the least stable genes. Both genes were revealed as the most reliable ones in Rhopalosphum padi non-viruliferous and viruliferous for Barley yellow dwarf virus 45 , supporting the notion that gene stability must be checked for each specie and condition. The pairwise variation analysis revealed that two reference genes perform trustfully (V2/3 < 0.15) across mite samples with and without CiLV-C, and hence adding more genes would not increase the confidence of the analysis (Fig. 2). We suggest using Ubiq and GAPDH in transcriptome studies with viruliferous and non-viruliferous Brevipalpus mites.
At B. yothersi mites from different development stages the Ubiq, EF1α, and GAPDH genes were identified as the most stable, whilst Tub and RPL13 were the least stable ones. Ubiq has been used to normalize RT-qPCR data from different developmental stages of the bettle Agrilus planipennis 46 . GAPDH was identified as a reliable reference gene for Spodoptera litura 47 and Sesamia inferens 48 caterpillars, and Lipaphis erysimi 35 thrips. EF1α, coding for a protein involved in the delivery of aminoacyl tRNAs to the acceptor site on ribosomes during protein translation 49 , was suggested as a reference gene for Aphis glycines 50 , Diabrotica virgifera virgifera 51 , Panonychus citri 23 , Locusta migratoria 52 , Frankliniella occidentalis 53 and Coleomegilla maculata 24 . Even though Tub and RPL13 were the least stable in B. yothersi mites from different development stages, they were pointed as the most stable ones in Lipaphis erysimi 35 , Tetranychus urticae 37 and T. cinnabarinus 22 . The pairwise variation analysis indicated the use of at least three reference genes to ensure reliable RT-qPCR analysis (V2/3 > 0.15) (Fig. 2). Thus, we suggest using Ubiq, EF1α, and GAPDH to normalize RT-qPCR data comparing Brevipalpus mites from distinct developmental stages.
In RefFinder analyses, we find the RPL13 and RPL32 ribosomal protein-coding genes as the most stable genes across B. yothersi populations fed on beans or sweet orange fruits. A limited number of studies have tested reference genes from arthropods feeding on distinct host plants or artificial diets. These works identified EF1a as the best reference for Drosophila melanogaster 31 and Lipaphis erysimi 35 , and GAPDH as the most stable in T. urticae populations 36 . Nevertheless, ribossomal protein genes have been extensively used for RT-qPCR normalization of data from a wide range of species and conditions 39,48,50,54 . The pairwise variation analysis indicated that the use of two reference genes is enough to ensure the reliability of RT-qPCR normalization (V2/3 < 0.15) in mites reared on beans or sweet orange fruits. Hence, we suggest using RPL13 and RPL32 in comparisons of Brevipalpus mite populations reared on different host plants.
Combining all the experimental conditions applied to B. yothersi mites, Ubiq, EF1α, and GAPDH were classified as the most stably expressed reference genes. The pairwise variation analysis indicated that using three genes (V3/4 < 0.15) improve the data confidence (Fig. 2) when mites samples comprising all sources of variability analyzed here (virus presence, different developmental stages and feeding from distinct host plants) are assayed. However, when the distinct experimental conditions are individually analyzed, we identified different sets of genes as the most stably expressed, reinforcing the idea that there is no universal reference gene for every situation and suggesting the use of the best gene set at each condition.
Finally, to validate the reliability of the reference genes identified here, we evaluated the normalized expression profile of V-ATPase gene in different development stages of B. yothersi. We found variable gene expression between treatments when we normalize the transcript levels using the most unstable gene (RPL13) or a combination of the two most unstable ones (RPL13/TUB) (Fig. 3). Results using inferior reference genes led to the conclusion that proto or deutonymphs express V-ATPAse at a higher level that the other mite phases. On the other hand, using the most stable gene or the three more stables ones combined (Ubiq/EF1α/GAPDH) to normalize the RT-qPCR data, no difference in V-ATPAse expression was observed between the development stages (Fig. 3).
These results indicate that the use of unstable reference genes might mask or introduce artificial changes in gene expression, leading to misinterpretations of biological phenomena. Similar results have been found in RT-qPCR www.nature.com/scientificreports www.nature.com/scientificreports/ studies from other systems, such as Drosophila melanogaster 55 or citrus plants 56 , emphasizing the need for testing a set of reference genes for each experimental condition assayed.
This study is the first evaluation of the expression stability of B. yothersi genes and provides a standardized procedure for gene expression analyses of false spider mites. The results presented here will allow the precise normalization in RT-qPCR studies involving Brevipalpus mites gene function and mite/virus interactions and will contribute to improve our understanding of the molecular mechanisms involved in the citrus leprosis and others viral diseases vectored by false spider mites.

Materials and Methods
Mite rearing and synchronization. B. yothersi isoline has been reared on healthy sweet orange fruits (Citrus sinensis L. Osbeck) at 25 ± 5 °C, 60 ± 10% RH and 14:10 h (light: dark) photoperiod conditions since 2008, on at the Acarology Laboratory from the Sylvio Moreira Citrus Center, Agronomic Institute (IAC), Cordeirópolis (São Paulo State, Brazil). To synchronize mite ages, adult females were placed on sweet orange fruits for egg-laying and removed after 24 h. Developmental stage. Five different stages (egg, larvae, protonymphs, deutonymphs and adult females) were collected from B. yothersi age-synchronized colonies onto healthy sweet orange fruits. Pools of 500 eggs, 400 larvae, 400 protonymphs, 300 deutonymphs, and 300 females were transferred directly to microcentrifuge tubes. To each developmental stage, three replicates were independently collected. The samples were flash frozen in liquid nitrogen and stored at −80 °C until use.
Mites host. B. yothersi mites were reared on plants of sweet orange and common bean (Phaseolus vulgaris L.) and maintained at 25 ± 5 °C, 60 ± 10% RH with a photoperiod of 14:10 h (light: dark). After at least three mite generations, pools of 300 females of each host plant were collected in three biological replicates. Samples were flash frozen and stored at −80 °C until use.
CiLV-C acquisition. Viruliferous samples were obtained by rearing the non-viruliferous false spider mites onto sweet oranges fruits with citrus leprosis symptoms for several generations. Three biological replicates with pools of the 300 viruliferous mites were collected and flash frozen stored at −80 °C. The CiLV-C acquisition was confirmed by RT-PCR using specific primers to detected a region within the viral movement protein (MP) gene, according to methodology described by Locali et al. 57 .
RNA extraction and cDNA synthesis. Total RNA was extracted from each development stage using RNAqueous ® -Micro Kit (Ambion, Life Technologies). The RNA concentration was measured with a NanoDrop ND-8000 spectrophotometer (Thermo Scientific) and RNA integrity was checked in 1% agarose gel. RNA samples with an A260/A280 ratio ranging from 1.8 to 2.2 were used for cDNA synthesis. First-strand complementary DNA was synthesized from 200 ng of total RNA with Platus Transcriber RNase H-cDNA First Strand kit (Sinapse Inc, catalog number S1402) following the manufacturer´s instructions and stored at −20 °C until use. The cDNA was diluted ten-fold for the subsequent RT-qPCR analysis.
Reference gene selection and primer design. Seven commonly reference genes used in arthropod studies were selected as candidates ( Table 2). Gene sequences were obtained from the B. yothersi database genome 15  Quantitative reverse transcription pCR. RT-qPCR was performed on 7500 Fast Real-Time PCR System (Thermo Scientific). Gene-specific primers ( Table 2) were used in reactions (12uL), containing 6.5 μL of GoTaq ® qPCR Master Mix (Promega), 3 μL of diluted cDNA template, 1.5 μL of ddH 2 O, and 0.5 μL of each specific primer (10 μM). The RT-qPCR program included an initial denaturation for 3 min at 95 °C followed by 40 cycles of 95 °C for 10 s, 60 °C for 30 s, and extension for 30 s at 72 °C. For melting curve analysis, a dissociation step cycle (55 °C for 10 s and then a gradual increase of 0.5 °C for 10 s until 95 °C) was added. Quantitative PCR reactions were carried out in 96-well plates with three technical replicates and three biological replicates. Three negative controls (no template) were included. The Cq values and PCR efficiency of each reaction were determined using Real-Time PCR Miner software 58 , a tool based on the kinetics of individual RT-qPCR reactions without the need of a standard curve. The arithmetic mean of amplification efficiencies of each primer pair was used for data normalization.
Gene expression stability. The stability of the candidate reference genes was evaluated by algorithms GeNorm 19 , Normfinder 27 , BestKeeper 28 , and the comparative ΔCq method 26 . The seven genes analyzed were then compared and ranked the tested candidates based on a web-based analysis tool RefFinder (http://leonxie. esy.es/RefFinder/) that integrates the four methods. The algorithms GeNorm e NormFinder used in this study are included in the software GenEx version 6.0.1 (www.multid.se). The optimal number of reference genes for accurate normalization was calculated by pairwise variation using GeNorm in qBaseplus software (https://www. qbaseplus.com/). www.nature.com/scientificreports www.nature.com/scientificreports/ The comparative ΔCq method establish a ranking based on pairwise of the genes by mean ΔCq values within a particular treatment and the stability of the gene is inversely proportional to its average standard deviation value (SD), and the gene with the lowest value is considered the most stable of Cq values between samples 26 . The GeNorm uses the pairwise variation between a specific candidate compared to the others to calculate an expression stability value (M-value). The lower the M-value, the higher stability of the gene, and an M value less than 1.05 is recommended 19 . The NormFinder calculates gene stability for all samples in any number of groups based on intra and inter group variations and combines values to ranking based on gene expression variation. The most stable gene will be the one with the lowest stability value (SV) according to the intra-group and inter-group variability of each gene 27 . Finally, the BestKeeper determines the most stably genes based on Pearson´s correlation coefficient for each gene, using raw data (Cq values) and PCR efficiency (E), with values closer to 1 indicating highest stability 28 .
Validation of reference genes. Genes were evaluated comparatively as its performance in data normalization of RT-qPCR, using as target the expression levels of the V-ATPase in different development stages of B. yothersi. Five separate normalization factors (NFs) were calculated based on: i) a single reference gene with the lowest GEOMEAN value, ii) the geometric mean of the three most stable genes with the lowest GEOMEAN value (as determined by RefFinder) and iii) the geometric mean of two most unstable genes with the highest GEOMEAN value (determined by RefFinder). The relative quantification of the expression of V-ATPase normalized by each of these NFs was calculated using the derived 2 −ΔΔCq method 59 .  Table 2. Gene name and primer information for RT-qPCR assays. *Mean efficiencies were calculated using Miner.