The analysis of reference genes expression stability in susceptible and resistant Apera spica-venti populations under herbicide treatment

Weed resistance to herbicides constitutes a serious problem to world crop production. One of the weeds that are significantly threatening the crops’ yield and quality is Apera spica-venti. The target-site resistance (TSR) mechanism of A. spica-venti has been widely studied, though, little is known about its non-target-site resistance (NTSR) mechanisms at the molecular level. Molecular examination of NTSR is, to a great extent, based on the expression profiles of selected genes, e.g. those participating in detoxification. However, to obtain reliable results of gene expression analysis, the use of a normalizer is required. The aim of this study was to select the best reference genes in A. spica-venti plants of both populations, susceptible and resistant to ALS inhibitor, under treatment with herbicide. Eleven housekeeping genes were chosen for their expression stability assessment. The efficiency correction of raw quantification cycles (Cq) was included in the gene expression stability analyses, which resulted in indicating the TATA-box binding protein (TBP), glyceraldehyde-3-phosphate dehydrogenase, cytosolic (GAPC), and peptidyl-prolyl cis–trans isomerase CYP28 (CYP28) genes as the most stably expressed reference genes. The obtained results are of vital importance for future studies on the expression of genes associated with the non-target-site resistance mechanisms in the A. spica-venti populations susceptible and resistant to herbicides.


Results
Plant populations selection for analysis of the expression stability of housekeeping genes. Eleven populations (two susceptible and nine resistant) were examined in order to determine ED50 (effective dose of active ingredient (ai) causing a 50% reduction in plant biomass). In the susceptible populations, ED50 was lower than 1/16N dose (N-the maximal recommended dose of the herbicide (120 g ha −1 , i.e. 9 g ha −1 of ai)) ( Table S1). For four out of the nine resistant populations the ED50 values were higher than 32N dose of the herbicide. In order to choose A. spica-venti populations to be subjected to HKGs stability analysis, plants from the populations resistant and susceptible to pyroxsulam were examined for the absence of mutations in the ALS sequence. ALS sequence analysis permitted selection of two susceptible and two resistant A. spica-venti populations without known mutations conferring the resistance. Single mutations were found, however, they were present in both susceptible and resistant plants, or the interdomain region of ALS.
RT-qPCR primer specificity and efficiency. In order to choose the most reliable reference genes in susceptible and resistant to ALS inhibitors plants of A. spica-venti exposed to treatment with herbicide, 11 candidate HKGs were selected for the expression stability analysis. These candidates were the primary metabolism genes frequently used as reference ones in other weeds species and were indicated as stable in other plants exposed to different abiotic stresses. The specificity of the primers used was confirmed by a single peak of melting curves ( Fig. 1) and agarose gel electrophoresis (Fig. 2). The sequences of RT-qPCR products were compared to the previously published A. spica-venti transcriptome data (NCBI Gene Expression Omnibus GSE86989) 16 showing high identity to the corresponding genes in the transcriptome. These sequences were deposited in the GenBank database (accession numbers provided in Table 1). To confirm the similarity of the obtained sequences to the  www.nature.com/scientificreports/ selected genes, the latter were blasted in the following databases: NCBI blastn and UniProt blastp, as well as, against model plant species Arabidopsis thaliana and Brachypodium distachyon. The results with the highest scores are presented in Table S2. The analysis of standard curves resulted in the amplification efficiency ranging from 88.2% (RPL23A) to 108.2% (EF1A), while the regression coefficient varied from 0.994 (ACT ) to 0.999 (GAPC, RPL23A, RCA , and TBP) ( Table 1). The distribution of Cq values revealed the highest expression of RCA with the mean Cq value of 18.82, while the lowest mean Cq value was observed for ARF1 (27.19) (Fig. 3). The lowest Cq range (5.11 cycles) was detected for the TBP, whereas the greatest discrepancy in Cq values (11.15) was observed for EF1A. .
To gain a comprehensive summary of the results, RefFinder software (https:// www. heart cure. com. au/ reffi nder) was used.
Most of the gene expression stability analysis softwares omit the fact that not every primer pair and reaction conditions ensure 100% efficiency in RT-qPCR. The use of qPCR efficiency corrected Cq values in statistical analyses of the most stably expressed genes impacts the final expression stability ranking, which may result in selection of different sets of reference genes for normalization. Therefore, the efficiency correction was applied to the raw Cq values.

Candidate reference genes expression stability analysis in all samples indicated that TBP, GAPC, and CYP28
show the most stable expression under tested conditions. The first analysis, performed with the geNorm software, using Bioconductor "NormqPCR" package in R software, indicated that the first three of the most stably expressed genes were TBP, ARF1, and GAPC, while the genes with the least stable expression were-EF1A, UBQ, and RPL23A ( Table 2). The last three genes also exhibited unstable expression according to BestKeeper and Nor-mFinder analyses. The genes found to be the most stably expressed by BestKeeper were TBP, ACT , and GAPC, whereas the ones CYP28, TBP, and GAPC were indicated as the most stably expressed by NormFinder. Similarly, the ΔCt method pinpointed CYP28, TBP, and GAPC to have the highest expression stability, contrary to RCA, EF1A, and RPL23A that were classified as the least stably expressed. The general ranking obtained from Ref-Finder analysis indicated TBP, GAPC, and CYP28 as the reference genes with the most stable expression in the plants susceptible and resistant to pyroxsulam, and EF1A, UBQ, and RPL23A as the least stably expressed ones.
Analysis of the reference genes expression stability separately for the pyroxsulam susceptible and resistant plants indicated that TBP was among the most stably expressed genes in both types of samples. The same methods were used in the analyses of HKGs expression stability, separately in the plants showing the susceptibility or resistance to sulfonamide herbicides.
In the resistant plants (Table 2), the genes classified to be the most stably expressed were as follows: TBP, ARF1, and ACT (according to geNorm); GAPC, ACT , and TBP (according to BestKeeper); GAPC, TBP, and CYP28 (according to NormFinder); CYP28, TBP, and GAPC (according to ΔCt method). The genes identified as the least stably expressed were: EF1A, RPL23A, and UBQ (according to geNorm); EF1A, RPL23A, and SPS (according to BestKeeper); EF1A, SPS, and RPL23A (according to NormFinder); RCA , EF1A, and ARF1 (according to ΔCt method). RefFinder analysis indicated GAPC, TBP, and ARF1 as the best reference genes with the most stable expression, while EF1A, RPL23A, and SPS showed the least stability in their expression. www.nature.com/scientificreports/ Statistical analyses of efficiency adjusted Cq values has changed the positions of the candidate reference genes in the expression stability ranking. The efficiency correction to the raw Cq values obtained from RT-qPCR was applied using the GenEx software. The values of efficiency used are given in Table 1. The gene expression stability analysis using the corrected Cq values was performed with all statistical methods used in the previous analyses namely: geNorm, BestKeeper, NormFinder, ΔCt, and RefFinder. This approach resulted in differences in the gene expression stability ranking (Table 3). Then, two genes expression stability rankings based on the data resulting from analyses of efficiency corrected Cq and non-corrected Cq were compared. No differences in the results obtained from geNorm analyses were observed. The results of the BestKeeper and NormFinder softwares revealed the changes in 39% (13/33) and 33% (11/33) of HKGs positions, respectively. The biggest discrepancies in the rankings were observed in the results obtained by the ΔCt method, with 73% (24/33) of changes in genes positions. The comprehensive ranking obtained by RefFinder showed that 45% (15/33) HKGs positions were altered. In all samples analyzed (Table 3), the three most stably expressed genes taking into account the efficiency correction were as follows: ARF1, TBP, and GAPC (according to geNorm); TBP, ACT , and ARF1 (according to BestKeeper); CYP28, TBP, and ARF1 (according to NormFinder); TBP, GAPC, and ACT (according to ΔCt method); TBP, CYP28, and GAPC (according to RefFinder). The ranking of the least stably expressed genes changed only according to the results of ΔCt method and comprised RCA , RPL3A, and EF1A.
In the samples collected from the susceptible plants, the three most stably expressed genes taking into account the efficiency correction were as follows: CYP28, SPS, and TBP (according to geNorm); TBP, ACT , and SPS (according to BestKeeper); CYP28, SPS, and TBP (according to NormFinder); TBP, ACT , and GAPC (according to ΔCt method); CYP28, SPS, and GAPC (according to RefFinder). Changes in the order of the unstably expressed genes were found in the results obtained from NormFinder (for genes EF1A, UBQ, and UBC), ΔCt (for RPL23A, RCA , and SPS) and RefFinder (for EF1A, UBQ, and RPL23A). www.nature.com/scientificreports/ The greatest differences in the initial three positions were found in the results obtained for the samples harvested from the resistant plants. The most stably expressed genes after efficiency correction ranked as follows: TBP, ARF1, and ACT (according to geNorm); ACT , TBP, and ARF1 (according to BestKeeper); GAPC, TBP, and CYP28 (according to NormFinder); TBP, GAPC, and UBC (according to ΔCt method); TBP, GAPC, and ACT (according to RefFinder). However, the positions of the most as well as the least stably expressed genes were unaltered in the NormFinder ranking. The unstably expressed genes in the rankings were changed accordingly to the following methods: BestKeeper (EF1A, SPS, and UBQ), ΔCt method (RCA , RPL23A, and EF1A), and RefFinder (EF1A, SPS, and RPL23A).
Two reference genes are required for the normalization of target genes expression results. Pairwise variation (V n/n+1 ) between the two sequential normalization factors (NF n and NF n+1 ) was calculated using geNorm. The threshold value of pairwise variation (V) below 0.15 indicates the minimum number of reference genes needed for the normalization of the target gene. In all samples collected from the plants susceptible and resistant to the herbicide, the target gene is required to be normalized using two reference genes. The V 2/3 value was lower than 0.15, regardless of the experimental variant (Fig. 4).

Validation of the selected reference genes revealed the constitutive ABCC10 and CYP89A2
up-regulation of their expression in plants resistant to pyroxsulam. ABC transporters as well as CYP450s are found to be involved in the metabolism-based NTSR to ALS inhibitors 21 . ABCC10 was highly expressed in Myosoton aquaticum plants resistant to ALS inhibitor 22 , while CYP89A2 was identified to confer metabolism-based diclofop resistance in Lolium rigidum 23 . To normalize the ABCC10 and CYP89A2 expressions, the selected pairs of reference genes were chosen according to the efficiency corrected results. The analysis was performed using the most stably (TBP and GAPC) and the least stably (EF1A and UBQ) expressed HKGs in www.nature.com/scientificreports/ all samples, to show the differences in target genes expression between the pyroxsulam resistant and susceptible plants.
The descriptive statistics of the results obtained with TBP and GAPC used as normalizers showed that the expression of ABCC10 was higher in the non-treated resistant plants than in susceptible plants (Fig. 5A). Although the change in the ABCC10 expression was higher in the pyroxsulam resistant plants 1 h after treatment with herbicide than in the susceptible ones, however, at 24 h time point, the ABCC10 expression rapidly rose in the susceptible plants. Similarly, the change in the CYP89A2 expression in the treated and non-treated resistant plants was higher than in the treated and non-treated susceptible ones, both at 1 h and 24 h (Fig. 5B). Moreover, the treatment with herbicide caused up-regulation of the CYP89A2 expression.
Finally, to show how important it is to choose the most stably expressed reference genes for accurate normalization, the ABCC10 and CYP89A2 expression analyses were performed using the least stably expressed HKGs as normalizers, namely EF1A and UBQ. ABCC10 was expressed at the same level in the majority of experimental conditions, except for the samples from the non-treated resistant plants harvested 24 h after treatment with herbicide, which exhibited lower ABCC10 expression than the remaining samples (Fig. 5A). The CYP89A2 expression analysis confirmed slight differences in the expressions between the resistant and susceptible plants, except for the samples from the susceptible plants harvested 24 h after treatment with herbicide with downregulation in CYP89A2 expression (Fig. 5B). Moreover, higher variance values calculated for each condition with the use of the most and the least stably expressed HKGs as normalizers, were observed for EF1A and UBQ, which indicated greater discrepancies of the obtained results within each subset of samples after normalization (Table S3). This analysis showed that the use of inappropriate normalizers changes the final results of the target gene expression analysis. Comparison of the results obtained using the most and the least stably expressed genes for normalization showed that there was not one direction of the changes in gene expression between samples. For example, the use of the least stably expressed genes for normalization led to the results indicating downregulation of the CYP89A2 expression in the non-treated susceptible plants harvested 24 h after treatment with herbicide, while the normalization performed with the most stably expressed genes indicated up-regulation of the CYP89A2 expression.

Discussion
Weed resistance to herbicides poses a serious threat to world crop production. Repetitive field applications of herbicides with the same mode of action resulted in the emergence of herbicide-resistant weed populations. The herbicide resistance may be initiated by mutations in the coding sequence of the target site of the herbicide active substance or over-production of the target enzyme. The second resistance mechanism is based, among others, on the increased metabolism and detoxication of xenobiotics. In the case of A. spica-venti, which was the subject of this study, TSR has been widely examined for the presence of mutations conferring resistance to ALS inhibitors 6,24,25 , while the experiments concerning NTSR mechanisms at the molecular level have been scarce. A. spica-venti transcriptome has been recently published 16 , however, to study the expression of the selected genes that may be associated with herbicide resistance, an internal control (reference gene) for normalization should be applied. Reference genes use is desired for the avoidance of the errors originating from the differences in RNA integrity, starting RNA quantity between samples, reverse transcription efficiency, and primer efficiency 18 . Therefore, in this study, A. spica-venti HKGs expression stability analysis was performed to choose the most stably expressed reference genes for the further target gene expression analyses.
According to the statistical analyses of the efficiency corrected Cq values, the best genes for expression data normalization were TBP, GAPC, and CYP28. The least stably expressed genes included EF1A, UBQ, and RPL23A. However, when the plants susceptible and resistant to ALS inhibitor were analyzed separately, the genes with the most stable expression were CYP28, SPS, and TBP, GAPC, respectively. TBP is the most stably expressed gene in Lolium multiflorum under salt stress and in A. fatua under herbicide stress 14,26 . Moreover, TBP was indicated as www.nature.com/scientificreports/ the most suitable reference gene under treatment with herbicide in Galium aparine 27 . GAPDH was also among the most stably expressed genes in herbicide-resistant grass species such as A. myosuroides, A. fatua, and Alopecurus japonicus 10,14,28 . CYP28, belonging to the cyclophilin (peptidyl-prolyl cis-trans isomerase) gene family was shown in our studies to be one of the most stably expressed reference genes. Cyclophilin was also used for normalization in Cucumis melo 29 , whereas, SPS was implemented as an endogenous reference gene in Oryza sativa 30 .
Various biotic and abiotic factors may influence housekeeping genes expression. Also, the expression stability of indicated experimentally genes might be affected in specific conditions including the plants' developmental stage, which should be taken into account when designing an experiment using reference genes for normalization.
The validation of the most and the least stably expressed reference genes showed that the selection of the inappropriate genes in normalization might lead to the incorrect results of target gene expression analysis in the www.nature.com/scientificreports/ test conditions. In our study, it was revealed that the ABCC10 expression in the pyroxsulam susceptible plants 24 h after herbicide application, varied substantially when normalization was carried out with the best and the worst reference genes. The use of the genes with unstable expression also led to the results indicating the opposite gene expression changes (from up-to down-regulation of the target gene expression). The analysis of ABCC10, as well as CYP89A2 expression, showed the elevated expression of these genes after the treatment with herbicide, which implies their involvement in the plant's reaction to the herbicide application. ABC transporters are known for their involvement in herbicide resistance through sequestration of the herbicides and their metabolites 21 .
ABCC10 was highly expressed in the plants resistant to tribenuron-methyl M. aquaticum and indicated to play an essential role in NTSR 22 . By contrast, the transcriptome analysis in diclofop resistant L. rigidum revealed the CYP89A expression up-regulation in the resistant, non-treated with the herbicide plants 23 . In our results, the ABCC10 expression in the untreated A. spica-venti resistant to pyroxsulam plants was higher than in susceptible plants. This implies that this gene is constitutively expressed in resistant plants at a higher level than in susceptible plants, which is consistent with previous studies of M. aquaticum. However, CYP89A2, belonging to the CYP450s gene family participating in phase I of herbicide detoxification, also exhibits the increased expression in the non-treated resistant plants at two time points, which is also in accordance with transcriptome analysis of diclofop-resistant L. rigidum. Usually, in analyses of the expression stability of HKGs by different algorithms, no efficiency corrected Cq values have been used as an input. In such situations, the assumption of 100% primers efficiency is taken, however, in practice, such a value of primer efficiency is not always achievable. Our results show that the implementation of the efficiency corrected input data changed the gene's positions in the expression stability rankings, therefore it is important to consider the efficiency correction of the raw input data. This issue has also been brought up by another study, which indicated the necessity of performing efficiency correction of the results obtained in RT-qPCR with primers exhibiting efficiencies differing by more than 10% from the optimal 100% 15 . RNA extraction and cDNA synthesis. Plant material was ground in a mortar using liquid nitrogen.

Materials and methods
Total RNA was extracted using 500 µl of TriReagent solution (Thermo Fisher Scientific) followed by RNA precipitation with 2-propanol. The precipitate was washed with 70% ethanol, air-dried, and suspended in nucleasefree water. The RNA concentration and purity (the 260 nm/230 nm and 260 nm/280 nm values) of each sample were estimated using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific).
Total RNA (2.5 µg) was reverse transcribed using Maxima First Strand cDNA Synthesis Kit for RT-qPCR with dsDNase (Thermo Fisher Scientific). The cDNA samples were diluted with 30 µl of nuclease-free water. www.nature.com/scientificreports/ Gene selection, primer design, and RT-qPCR. Eleven candidate reference genes were selected for the analysis: ACT , ARF1, CYP28, EF1A, GAPC, RPL23A, RCA , SPS, TBP, UBC, and UBQ. Moreover, 2 genes associated with NTSR were chosen for validation: encoding the ABC transporter C family member 10-like (ABCC10) and cytochrome P450 89A2 (CYP89A2). Primers sequences were designed based on the A. spica-venti published transcriptome (NCBI Gene Expression Omnibus GSE86989) 16 using Primer3 software (version 0.4.0) 32 . Primers sequences list is presented in Table 1. RT-qPCR was performed in 10 µl solution containing 1× iTaq Universal SYBR Green Supermix (Bio-Rad, Hercules, USA), 0.5 µM forward primer, 0.5 µM reverse primer, and 1 µl cDNA. The PCR program consisted of an initial incubation step at 95 °C for 3 min, followed by 40 cycles of 95 °C for 15 s, an annealing step for 30 s (temperatures listed in Table 1), and 72 °C for 30 s. The dissociation curves were generated from 60 to 95 °C. The standard curve for each gene was made on the basis of n-fold dilutions of cDNA. Reactions were performed in triplicate (technical replicates) in QuantStudio5 (Thermo Fisher Scientific). However, to optimize reactions' conditions for each pair of primers, firstly RT-qPCR with temperature gradient at the annealing step was performed. To confirm the estimated size of the amplicons, the reactions' products were resolved on a 2% agarose gel, purified from the gel with Wizard ® SV Gel and PCR Clean-Up System (Promega), ligated to pJET1.2 plasmid using CloneJET PCR Cloning Kit (Thermo Fisher Scientific), and cloned into DH10B Escherichia coli competent cells. The plasmids were isolated from E. coli cells using NucleoSpin ® Plasmid (Mecherey-Nagel), followed by the confirmation of insert presence in plasmids by the digestion with BglII. DNA inserts were sequenced by Genomed (Warsaw, Poland). Sequencing data were analyzed using the BioEdit Sequence Alignment Editor 7.5.5 31 . The obtained sequences were compared with sequences of the corresponding transcript in A. spica-venti transcriptome and submitted to the GenBank database. To confirm the similarity of the sequences of products amplified by the designed primers targeting HKGs encoding genes to the selected genes, they were blasted to the following databases: NCBI blastn and UniProt blastp. Moreover, the sequences were blasted against a model plant species A. thaliana TAIR database (Araport11 Proteins) using the website BLASTP search (version 2.9.0+, https:// www. arabi dopsis. org/ Blast/ index. jsp) and against the closest to A. spica-venti model plant species B. distachyon protein sequences (v3.2) 33 downloaded from Phytozome 34 using OmicsBox 35 .
Data analysis and validation. Distribution of candidate reference genes Cq, efficiency correction of Cq data as well as the normalization of the expression of CYP89A2 and ABCC10 using the most and the least stably expressed genes were performed with GenEx (version 6.1.1.550, Multid Analyses AB, Göteborg, Sweden). The expression of CYP89A2 and ABCC10 was calculated using descriptive statistics at the confidence level of 95%.
To determine the expression stability of the reference gene candidates for A. spica-venti, the analyses were carried out using four independent bioinformatics algorithms: geNorm 17 and BestKepeer 18 softwares, performed in the integrated development environment (IDE) of R (RStudio), NormFinder software 19 in GenEx (version 6.1.1.550), and ΔCt method 20 recalculated in Microsoft Excel application (Microsoft Office 2018). To gain a comprehensive list of genes expression stability by calculating the geometric mean, RefFinder software (https:// www. heart cure. com. au/ reffi nder) was implemented. The analyses were performed using pooled Cq data from susceptible and resistant samples together, pooled Cq values for the susceptible samples only, and separately for the resistant samples only.
To conduct genes expression stability analysis in the R software environment, the data from RT-qPCR analysis (non-corrected or efficiency corrected Cq values) were loaded to R software (version 3.6.2) 36 . To conduct expression stability analysis of the selected candidate genes, two different packages were installed: "NormqPCR" (version 1.7.1) 37 (geNorm) and "ctrlGene" (version 1.0.1) 38 (BestKeeper). The first one calculates the average expression stability values (M) of the analyzed genes during stepwise exclusion of the least stably expressed genes in each round until the most stably expressed candidates remain. The second one calculates descriptive statistics from Cq data and pairwise correlation between all analyzed genes. The results of gene expression stability analysis from BestKeeper software are calculated given the standard deviation (SD) of Cq values between all analyzed genes, where the lowest SD score represents the gene with the most stable expression in ranking index 37 . In order to determine the minimum required number of reference genes for the analysis, the algorithm in geNorm calculates the pairwise variation (V n/n+1 ) between the two sequential normalization factors (NF n and NF n+1 ) with the cut-off threshold of 0.15 38 .
The NormFinder software calculates a global average expression of all genes in the model studied, to which the individual genes are compared. Next, the SD is estimated for each candidate gene and separates the variation into intragroup and intergroup contribution 19 . The ΔCt method calculates the expression stability of individual genes by comparing the relative expression of two reference genes pairwise. At the first step, the delta Cq is estimated between each pair of tested genes with the SD values of the obtained results. Then, a fixed base index of SD values for each tested-gene with its all possible pair is generated. The expression stability of candidate genes is ranked according to the mean values of SD recalculated for each index 20 .
Consent for publication. All authors have consented to this publication.

Data availability
All data generated or analyzed during this study are included in this published article.