Carex muskingumensis and Osmotic Stress: Identification of Reference Genes for Transcriptional Profiling by RT-qPCR

Carex muskingumensis is a highly valued perennial ornamental grass cultivated worldwide. However, there is limited genetic data regarding this species. Selection of proper reference genes (RGs) for reverse transcription quantitative PCR (RT-qPCR) data normalization has become an essential step in gene expression analysis. In this study, we aimed to examine expression stability of nine candidate RGs in C. muskingumensis plants, subjected to osmotic stress, generated either by salinity or PEG treatment. The identification of genes exhibiting high expression stability was performed by four algorithms (geNorm, NormFinder, BestKeeper and deltaCt method). The results showed that the combination of two genes would be sufficient for reliable expression data normalization. ADP (ADP-ribosylation factor) and TBP (TATA-box-binding protein) were identified as the most stably expressed under salinity treatment, while eIF4A (eukaryotic initiation factor 4A) and TBP were found to show the highest stability under PEG-induced drought. A set of three genes (ADP, eIF4A and TBP) displayed the highest expression stability across all experimental samples tested in this study. To our best knowledge, this is the first report regarding RGs selection in C. muskingumensis. It will provide valuable starting point information for conducting further analyses in this and related species concerning their responses to water shortage and salinity stress.


Introduction
Carex L. is the largest and the most widespread genus in the Cyperaceae family [1]. Carex muskingumensis Schwein., also known as Muskingum sedge or Palm sedge, is a highly valued perennial ornamental grass cultivated worldwide. It is used as a ground cover or planted in flowerbeds or containers. The species is rather easy to grow in an average, medium to wet soil, in full sun or part shade [2].
During its life cycle plants are exposed to various environmental stresses. Water deficits, as well as soil salinization are currently becoming increasingly a problem in various regions of the world. Both drought and high salinity disrupt osmotic balance in plant cells and in consequence limit plant growth and development [3,4]. Therefore, knowledge of plants molecular response to such conditions is of great importance. In order to better understand plants tolerance adaptation mechanisms, many studies have focused on monitoring changes in genes expression under PEG-induced drought [5][6][7] or salinity treatments [8][9][10].
Reverse transcription quantitative PCR (RT-qPCR) constitutes commonly used approach in gene expression studies. However, in order to avoid erroneous results and consequently incorrect

Experiment Design
The plant material were C. muskingumensis tufts obtained from a stable tissue culture grown on the solidified Murashige and Skoog (MS) [30] medium supplemented with kinetin (KIN) in concentration of 2.5 mg·dm −3 and indole-3-acetic acid (IAA) in concentration of 0.25 mg·dm −3 . The pH was adjusted to 5.7. Explants for the experiment were 3-4 cm long single shoots with 4-5 leaves. The explants were placed on the MS medium supplemented with 1.5% or 3% of NaCl for salinity treatment and 2% or 8% of PEG (polyethylene glycol) for drought treatment. The medium without NaCl or PEG (0%) was treated as a control. Each treatments consisted of 4 flasks with 5 shoots. The explants were grown in the temperature of 22 • C day/20 • C night and 16 h photoperiod with 35 µmol·m −2 ·s −1 light intensity.

The RNA Extraction and Reverse Transcription
The RNA extraction was conducted after 7 days of cultivation, in four biological replicates for each treatment. Immediately after harvesting leaves from the treated and control plants, the tissue (100 mg) was homogenized in liquid nitrogen using mortar and pestle. Total RNA isolation was conducted with the use of TRIzol reagent (invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer's instructions. In order to prevent RNA degradation ribonuclease inhibitor (EURx, Gdansk, Poland) was added to all samples to final concentration of 1 U/µL. The RNA concentration was assessed with NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The integrity of RNA samples was analysed by the means of electrophoresis in 2% agarose gel stained with ethidium bromide.
Prior to cDNA synthesis, genomic DNA was removed by DNase I (EURx) treatment. The reverse transcription was performed with NG dART RT kit (EURx) according to the supplier's recommendations.
The components of reaction mixture were as following: 0.5 µg RNA, 4 µL of 5 × NG cDNA Buffer, 1 µL of 50 µM Oligo(dT) 20 primer, 1 µL of NG dART RT mix and RNase-free water to final volume of 20 µL.

RT-qPCR Reactions
Nine genes (ACT7, ADP, EF1a, eIF4A, GAPDH, PEPKR1, SAND, TUBa, TBP) that are more or less frequently used as internal controls for RT-qPCR data normalization were chosen as candidates for RGs selection. Since no sequences for C. muskingumensis were available, Carex rigescens sequences deposited in GeneBank or SRA database (Sequence Read Archive accession No. SRX2755644) [27] were used for primer designing. Arabidopsis thaliana and Zea mays nucleotide sequences were used as query sequences for performed BLASTN search. Primers for RT-qPCR were designed with PrimerBLAST tool [31] (Table S1).
The RT-qPCR reactions were based on SYBR Green chemistry. The reaction mixture (final volume of 25 µL) contained 10 ng of cDNA, 1 × SG/ROX qPCR Master Mix (EURx), 400 nM of each primer and 0.25 U uracil-N-glycosylase. The reactions were performed according to the following cycling program: 2 min at 50 • C, 10 min at 95 • C, 40 cycles of 15 s at 94 • C, 30 s at 60 • C and 30 s at 72 • C. All RT-qPCR reactions were conducted on a QuantStudio 3 apparatus (Applied Biosystems, Carlsbad, CA USA). In order to verify the specificity of PCR products melting curve analysis was performed after each run with continuous data collection from 60 to 95 • C. The reactions were performed in four biological replicates with two technical replicates for each sample. No-template controls were included in analyses. Standard curves were generated from serial dilution of pooled cDNA. The PCR efficiencies for each primer pair were determined according to the equation [10 (1/−S) − 1] × 100%, where S represents the standard curve slope value [25].

Analysis of Gene Expression Stability
The expression stability of tested genes was estimated using following algorithms: geNorm [23], NormFinder [24], BestKeeper [25] and deltaCt method [26]. According to the requirements, performed analyses were based, either on untransformed Cq values (for BestKeeper and deltaCt method) or relative quantities (for geNorm and NormFinder).
For the validation of selected reference genes, the expression level of target gene ascorbate peroxidase (APX) was analyzed under tested experimental conditions. The expression data of APX was normalized using most stable and least stable RGs according to the obtained results. The RT-qPCR amplification conditions were the same as those described above. The relative expression level of the target gene was calculated using the 2 −∆∆Ct method.

Growth of Shoots
C. muskingumensis shoots cultivated on the medium supplemented with PEG or NaCl had visible symptoms of the negative influence of the substances used. Explants grown on media, supplemented with PEG, were smaller than the control ones. It was noted that the higher concentration of PEG the smaller plants were obtained ( Figure 1). Explants treated with NaCl were definitely damaged after the first week of treatment. It was observed that leaves started to turn pale green after the first week of cultivation, changing into brown after the second and the third ones. On the medium supplemented with 8% NaCl most explants died after four weeks of cultivation ( Figure 2).

Candidate Reference Genes-Efficiency and Specificity of Amplification
Nine genes (ACT7, ADP, EF1a, eIF4A, GAPDH, PEPKR1, SAND, TUBa, TBP) were chosen as candidate RGs for normalization of gene expression data in C. muskingumensis plants subjected to NaCl and PEG treatments. Confirmation of designed primer pairs suitability for RT-qPCR analysis was based on the analysis of standard curves performed on pooled cDNA and melting curves. Determined reaction efficiency, slope and regression coefficient (R 2 ) are shown in Table S2. Amplification efficiency of each primer pair was above 90% and it ranged from 95.6% (ADP) to 111.2% (PEPKR1). Regression coefficients were all above 0.99 except for the EF1a where it reached 0.986. Specificity of each designed primer pair was confirmed by a single peak on dissociation curves ( Figure S1).

Reference Gene Selection
Nine candidate RGs were screened among all experimental samples in order to find the most stable ones in tested material. Tested genes displayed variation in their expression across all of the samples with Cq value ranging from 21.83 to 32.82 (Table S3). Out of all candidate RGs ACT7 had the highest expression with mean Cq value of 23.11, while PEPKR1 showed the lowest expression level with mean Cq value of 30.77. The widest range of expression (Cq difference of almost 5 cycles) suggesting low stability level was observed for TUBa gene.
The expression stability of selected candidate RGs was investigated with four different algorithms-geNorm, NormFinder, BestKeeper and deltaCt method. The analyses were conducted for three different datasets, that was data obtained from NaCl experiment, data obtained from PEG experiment and combined data of two abovementioned experiments.

Candidate Reference Genes-Efficiency and Specificity of Amplification
Nine genes (ACT7, ADP, EF1a, eIF4A, GAPDH, PEPKR1, SAND, TUBa, TBP) were chosen as candidate RGs for normalization of gene expression data in C. muskingumensis plants subjected to NaCl and PEG treatments. Confirmation of designed primer pairs suitability for RT-qPCR analysis was based on the analysis of standard curves performed on pooled cDNA and melting curves. Determined reaction efficiency, slope and regression coefficient (R 2 ) are shown in Table S2. Amplification efficiency of each primer pair was above 90% and it ranged from 95.6% (ADP) to 111.2% (PEPKR1). Regression coefficients were all above 0.99 except for the EF1a where it reached 0.986. Specificity of each designed primer pair was confirmed by a single peak on dissociation curves ( Figure S1).

Reference Gene Selection
Nine candidate RGs were screened among all experimental samples in order to find the most stable ones in tested material. Tested genes displayed variation in their expression across all of the samples with Cq value ranging from 21.83 to 32.82 (Table S3). Out of all candidate RGs ACT7 had the highest expression with mean Cq value of 23.11, while PEPKR1 showed the lowest expression level with mean Cq value of 30.77. The widest range of expression (Cq difference of almost 5 cycles) suggesting low stability level was observed for TUBa gene.
The expression stability of selected candidate RGs was investigated with four different algorithms-geNorm, NormFinder, BestKeeper and deltaCt method. The analyses were conducted for three different datasets, that was data obtained from NaCl experiment, data obtained from PEG experiment and combined data of two abovementioned experiments.

Candidate Reference Genes-Efficiency and Specificity of Amplification
Nine genes (ACT7, ADP, EF1a, eIF4A, GAPDH, PEPKR1, SAND, TUBa, TBP) were chosen as candidate RGs for normalization of gene expression data in C. muskingumensis plants subjected to NaCl and PEG treatments. Confirmation of designed primer pairs suitability for RT-qPCR analysis was based on the analysis of standard curves performed on pooled cDNA and melting curves. Determined reaction efficiency, slope and regression coefficient (R 2 ) are shown in Table S2. Amplification efficiency of each primer pair was above 90% and it ranged from 95.6% (ADP) to 111.2% (PEPKR1). Regression coefficients were all above 0.99 except for the EF1a where it reached 0.986. Specificity of each designed primer pair was confirmed by a single peak on dissociation curves ( Figure S1).

Reference Gene Selection
Nine candidate RGs were screened among all experimental samples in order to find the most stable ones in tested material. Tested genes displayed variation in their expression across all of the samples with Cq value ranging from 21.83 to 32.82 (Table S3). Out of all candidate RGs ACT7 had the highest expression with mean Cq value of 23.11, while PEPKR1 showed the lowest expression level with mean Cq value of 30.77. The widest range of expression (Cq difference of almost 5 cycles) suggesting low stability level was observed for TUBa gene.
The expression stability of selected candidate RGs was investigated with four different algorithms-geNorm, NormFinder, BestKeeper and deltaCt method. The analyses were conducted for three different datasets, that was data obtained from NaCl experiment, data obtained from PEG experiment and combined data of two abovementioned experiments.

geNorm Analysis
The geNorm algorithm ranks the candidate RGs according to the average expression stability value (M value). The lower the M value the more stable is the gene's expression. The M value threshold for RG selection was established at 1.5 by Vandesompele et al. [23]. All tested candidate genes in all analyzed datasets (Figure 3) had the M value lower than the 1.5 cutoff. The same three genes (eIF4A, ADP and TBP) displayed the most stable expression in combined dataset (PEG plus NaCl samples) as in NaCl experiment alone. The most stable RGs, found in PEG-treated samples, were eIF4A, TBP and PEPKR1. TUBa exhibited the least stable expression in all experimental conditions (regardless of the NaCl and PEG samples being analyzed separately or together).

geNorm Analysis
The geNorm algorithm ranks the candidate RGs according to the average expression stability value (M value). The lower the M value the more stable is the gene′s expression. The M value threshold for RG selection was established at 1.5 by Vandesompele et al. [23]. All tested candidate genes in all analyzed datasets (Figure 3) had the M value lower than the 1.5 cutoff. The same three genes (eIF4A, ADP and TBP) displayed the most stable expression in combined dataset (PEG plus NaCl samples) as in NaCl experiment alone. The most stable RGs, found in PEG-treated samples, were eIF4A, TBP and PEPKR1. TUBa exhibited the least stable expression in all experimental conditions (regardless of the NaCl and PEG samples being analyzed separately or together).

NormFinder Analysis
NormFinder algorithm generates a stability value (SV) for each candidate RG including in its estimations both intra-group and inter-group variation in the sample set. The lower the SV, the more stable is the gene's expression [24]. The best-scoring genes for samples subjected to NaCl stress were ADP, TBP and ACT7. In PEG treated samples the highest stability was displayed by PEPKR1, eIF4A and TBP (Table 1). When all of the sample sets were analyzed together the NormFinder ranked TBP, ADP and eIF4A as the most stable RGs for data normalization. Noteworthy, TBP gene was among three most stable RGs regardless of what data set was used for the calculations. On the other hand, the TUBa gene displayed the lowest expression stability in all of the data sets analyzed.

NormFinder Analysis
NormFinder algorithm generates a stability value (SV) for each candidate RG including in its estimations both intra-group and inter-group variation in the sample set. The lower the SV, the more stable is the gene's expression [24]. The best-scoring genes for samples subjected to NaCl stress were ADP, TBP and ACT7. In PEG treated samples the highest stability was displayed by PEPKR1, eIF4A and TBP (Table 1). When all of the sample sets were analyzed together the NormFinder ranked TBP, ADP and eIF4A as the most stable RGs for data normalization. Noteworthy, TBP gene was among three most stable RGs regardless of what data set was used for the calculations. On the other hand, the TUBa gene displayed the lowest expression stability in all of the data sets analyzed.

BestKeeper Analysis
BestKeeper software calculates its statistics based on Cq values. The best RGs are identified by assessing correlation coefficient (r) of each individual gene with the BestKeeper index (the geometric mean of all candidate genes). The most stably expressed genes are those exhibiting the highest coefficient of correlation [25]. Out of all tested candidate RGs, TBP was found to be the most stable, regardless of the treatment (Figure 4). High expression stability was also observed for ADP, eIF4A and EF1a, however, with interchanging ranking positions. Highest variation, on the other hand, was consistently displayed by PEPKR1 and TUBa genes.
Genes 2020, 11, x FOR PEER REVIEW 7 of 14 mean of all candidate genes). The most stably expressed genes are those exhibiting the highest 187 coefficient of correlation [25]. Out of all tested candidate RGs, TBP was found to be the most stable, 188 regardless of the treatment (Figure 4). High expression stability was also observed for ADP, eIF4A 189 and EF1a, however, with interchanging ranking positions. Highest variation, on the other hand, was 190 consistently displayed by PEPKR1 and TUBa genes.

deltaCt Method
The deltaCt method is an approach that identifies the most stable RGs by comparing relative expression of pairs of genes within each sample. Low average standard deviation (mean SD) reflects the low level of variability [26]. Based on the results generated from deltaCt method, highest stability of expression in all datasets was displayed by ADP, which was followed by either TBP or eIF4A (Table 2). Both TUBa and PEPKR1 could be considered as least stable in tested material, as they were characterized by the highest mean SD.

Determination of Reference Genes Optimal Number
The geNorm algorithm, additionally to ranking candidate RGs from the most to the least stable, allows for the determination of their optimal number for data normalization. The optimal number of RGs can be calculated by pairwise variation (V n /V n+1 ). If pairwise variation is below 0.15, the inclusion of additional RG does not make any significant contribution to the data analysis. In this study, V 2/3 was lower than the threshold value of 0.15 ( Figure 5) thus indicating that using just two best-performing RGs would be sufficient for expression data normalization. Nevertheless, Vandesompele et al. [23] and Pfaffl et al. [25] recommend the minimal use of at least three most stable RGs for reliable normalization. In the present study, when all samples were analyzed together, the same three genes were identified as the most stably expressed by all used algorithms, that is ADP, eIF4A, and TBP. However, their ranking positions slightly differed (Table 3). The ADP and TBP were among top three ranked genes in salinity treatment regardless of the calculation method. Both NormFinder and deltaCt method pointed out ACT7 as third most stable gene in NaCl treatment, whereas geNorm and BestKeeper indicated eIF4A, or EF1a, respectively.
According to all algorithms eIF4A and TBP were among the most stably expressed genes in drought treated plants. Moreover, PEPKR1 was indicated as stably expressed during PEG-induced drought by both geNorm and NormFinder. However, according to the results generated by BestKeeper and deltaCt method, the PEPKR1 gene showed high expression variability in this treatment. Instead, ADP was proposed as a gene showing relatively stable expression in this conditions. All algorithms indicated that both SAND and GAPDH showed medium to poor stability of expression. ACT7 and EF1a mostly displayed intermediate level of stability. The PEPKR1 gene, interchangeably with TUBa, were ranked by BestKeeper and deltaCt algorithms as two least stable genes in all sample sets. The geNorm and NormFinder invariably ranked TUBa at the last position, suggesting its high expression instability in tested conditions. Additionally, the expression of target gene ascorbate peroxidase (APX) was analyzed, using most stable and least stable RGs, according to obtained results. Based on the comprehensive evaluation of RGs stability, we normalized APX expression in NaCl treated samples with ADP, TBP, combination of ADP + TBP and TUBa. PEG treated samples were normalized with eIF4A, TBP, combination of eIF4A + TBP and TUBa. The expression profile of APX under salinity stress was consistent when normalized with RGs showing high stability ( Figure S2). However, a strong increase in APX expression could be observed when TUBa was used for data normalization. The expression of APX under drought treatment was at steady level with eIF4A, TBP or combination of eIF4A + TBP as a reference, whereas normalization with TUBa resulted in underestimation of APX expression ( Figure S2). This confirms that selection of reliable RGs is crucial for proper target gene expression analysis.

Discussion
The influence of specific treatments on candidate reference genes (RGs) expression has been well-presented in studies testing various stressing factors [16][17][18]32]. Depending on the conditions, applied candidate gene might exhibit, either stable or varied expression, thus, it might be useful or useless for data normalization. Since using unstable RGs may lead to erroneous results and incorrect conclusions, appropriate RT-qPCR data normalization has become necessity [13,22,33]. In this study, we investigated nine potential RGs (ACT7, ADP, EF1a, eIF4A, GAPDH, PEPKR1, SAND, TUBa, TBP) in terms of their stability in C. muskingumensis plants exposed to salinity and PEG-induced drought.
The results of our experiment showed that the combination of ADP and TBP genes was the most suitable for normalization of data obtained from C. muskingumensis plants subjected to NaCl stress. Only one of the algorithms, used in this study, presented partially different results, selecting ADP with eIF4A as the most stable pair of genes. Still, TBP performed well and was identified as third in the ranking. More inconsistent results were found in C. muskingumensis plants exposed to PEG-induced drought. Here different algorithms selected different pairs of genes. Yet, each pair included eIF4A, thus, confirming its high stability under water deficit. Among other genes identified as highly stable were TBP (selected by two algorithms), ADP and PEPKR1. Partially inconsistent results, generated through various programs, are due to distinct statistical algorithms they use in stability calculations. Such differences in ranking positions were also reported in other studies [15,22,34].
More consistency in current study was observed in rankings obtained when both drought and salinity sample sets were considered together. Although with different ranking positions, the same three genes (ADP, eIF4A, TBP) were identified by all used algorithms as the most stable across tested conditions. Consequently, we believe that they constitute strong candidates for RT-qPCR data normalization. We suggest that these genes are included in the gene expression analyses of C. muskingumensis subjected to osmotic stress.
The ADP, eIF4A and TBP have been previously identified as stable RGs in several plant species subjected to drought and salinity stresses. ADP performed well or very well in wild barley (Hordeum brevisubulatum) under different treatments that generated osmotic stress in cells [22]. TBP was found to be highly stable in drought treated celery (Apium graveolens) and moderately stable under salt stress [34]. The eIF4A was first in comprehensive rankings generated for perennial ryegrass (Lolium perenne) subjected to drought and high salinity stresses [35].
Study performed on close relative, C. rigescens, recommended eIF4A together with PEPKR1 and GADPH to be preferentially selected for RT-qPCR analyses of that species [32]. To some extent, this is in concordance with our results, as eIF4A and PEPKR1 ranked well or very well in the results generated by geNorm and NormFinder algorithms, particularly for PEG-treated samples. Both eIF4A and PEPKR1 were among top three most stable genes in these rankings. Nevertheless, GAPDH did not perform as good as it might have been expected in the conditions of our experiment.
Both GAPDH and ACT have been commonly used as RGs for RT-qPCR data normalization in plants [36][37][38][39][40]. However, it has been noticed that among many studies testing larger panels of potential RGs, few found GAPDH or ACT to show high stability of expression [13]. ACT and GAPDH exhibited variable expression in leaves of bermudagrass (Cynodon dactylon) under NaCl and PEG stress [18]. Likewise, both genes performed poorly in seashore paspalum (Paspalum vaginatum) exposed to drought and salinity [19]. Furthermore, GAPDH was reported to be least stable in annual ryegrass (Lolium multiflorum) under abovementioned stresses [16], while ACT was least stable in durum wheat (Triticum durum) subjected to drought [33].
In our study, ACT7 performed slightly better than GAPDH in terms of expression stability in NaCl experiment, contrary to PEG-induced drought experiment, where the situation was reversed. Both genes, however, are considered to be far from optimal for C. muskingumensis data normalization, as they showed rather medium to poor expression stability in most rankings. Moreover, according to most of the algorithms TUBa, which is another 'traditional' RG, proved to be the least stable of all tested RGs in conditions of our study. High variability of TUB expression in response to drought and salt treatment was also reported in leaves of creeping bentgrass (Agrostis stolonifera) [17], as well as leaves of bermudagrass (C. dactylon) [18].
It has been agreed that RGs should be selected individually for each experiment [12,15,16]. However, prior to performing RGs selection, the type of candidate genes tested should be carefully considered. Apart from using a few 'traditional' RGs (such as GAPDH, ACT, TUB, EF1a), it is worth analyzing less frequently used genes or even novel genes, as they may turn out to be more stable under given conditions. For instance, Dudziak et al. [41] analyzed expression stability of five traditional RGs together with five novel candidate genes in common wheat (Triticum aestivum) lines subjected to drought stress. Both geNorm and NormFinder algorithms found novel gene, CJ705892, to be the most stable. Similarly, Liu et al. [16] tested six common RGs together with four novel candidate genes in annual ryegrass (L. multiflorum). They found that under drought stress one of the 'traditional' RGs (eIF4A) was the most stable. However, under saline-alkali stress highest stability was exhibited by one of the novel genes (Unigene14912). Therefore, while satisfactory results might be obtained by using conventional RGs (as it was done within this study), considering some alternatives might bring about valuable outcomes.
The fact that gene expression analysis using RT-qPCR technique should be preceded by identification of most stable RGs within its sample set is nonnegotiable [12,13]. However, the results of other studies testing panels of potential RGs in the same, or closely related, species can provide valuable starting point information. Such studies might indicate which candidate RGs are worth being tested for expression stability. This is especially true for non-model organisms that are not yet well characterized [13].
Our experiment identified most stable RGs in C. muskingumensis plants subjected to osmotic stress. The data suggests that ADP and TBP would be the most suitable to use in salinity studies, while eIF4A in combination with TBP could be used in drought experiments. Across all tested conditions highest expression stability was displayed by ADP, eIF4A and TBP. To our best knowledge this is the first report regarding RGs selection in this species. The identification of stable RGs in C. muskingumensis will be helpful for accurate RT-qPCR data normalization and will facilitate further investigations of molecular mechanisms involved in its response to high salinity and drought.