Next Article in Journal
Critical Issues in the Study of Magnesium Transport Systems and Magnesium Deficiency Symptoms in Plants
Next Article in Special Issue
Regulatory Proteolysis in Arabidopsis-Pathogen Interactions
Previous Article in Journal
Early Pregnancy Biomarkers in Pre-Eclampsia: A Systematic Review and Meta-Analysis
Previous Article in Special Issue
Physcomitrella patens Activates Defense Responses against the Pathogen Colletotrichum gloeosporioides
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Reference Genes in the Pathosystem Phakopsora pachyrhizi/ Soybean Suitable for Normalization in Transcript Profiling

Department of Phytopathology, Institute of Phytomedicine, Faculty of Agricultural Sciences, University of Hohenheim, Otto-Sander-Straße 5, 70599 Stuttgart, Germany
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Int. J. Mol. Sci. 2015, 16(9), 23057-23075; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms160923057
Submission received: 19 August 2015 / Revised: 14 September 2015 / Accepted: 18 September 2015 / Published: 23 September 2015
(This article belongs to the Special Issue Plant Microbe Interaction)

Abstract

:
Phakopsora pachyrhizi is a devastating pathogen on soybean, endangering soybean production worldwide. Use of Host Induced Gene Silencing (HIGS) and the study of effector proteins could provide novel strategies for pathogen control. For both approaches quantification of transcript abundance by RT-qPCR is essential. Suitable stable reference genes for normalization are indispensable to obtain accurate RT-qPCR results. According to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines and using algorithms geNorm and NormFinder we tested candidate reference genes from P. pachyrhizi and Glycine max for their suitability in normalization of transcript levels throughout the infection process. For P. pachyrhizi we recommend a combination of CytB and PDK or GAPDH for in planta experiments. Gene expression during in vitro stages and over the whole infection process was found to be highly unstable. Here, RPS14 and UbcE2 are ranked best by geNorm and NormFinder. Alternatively CytB that has the smallest Cq range (Cq: quantification cycle) could be used. We recommend specification of gene expression relative to the germ tube stage rather than to the resting urediospore stage. For studies omitting the resting spore and the appressorium stages a combination of Elf3 and RPS9, or PKD and GAPDH should be used. For normalization of soybean genes during rust infection Ukn2 and cons7 are recommended.

1. Introduction

Phakopsora pachyrhizi, the causative agent of Asian soybean rust, is one of the most devastating pathogens of soybean. The pathogen originates from East Asia, but by way of Hawaii and Australia it spread to South America where it now also causes enormous damages, especially in Brazil. In severely infected fields the pathogen can cause yield losses of up to 80%. Since 2004 the pathogen is also present in the continental United States, the biggest soybean producer worldwide, threatening production there [1].
P. pachyrhizi belongs to the Pucciniales, the rust fungi. It is an obligate biotrophic pathogen, which means that it is entirely dependent on its host plant for nutrient supply. Rust fungi also only sporulate growing on their respective host plant [2]. In contrast to other rust fungi that produce up to five different spore forms, for P. pachyrhizi infection only from urediospores has been observed. On the host leaf the urediospore attaches to the surface, germinates, forms an appressorium, and then the fungus penetrates directly through the epidermis—a feature unique to P. pachyrhizi. The epidermal cell that is penetrated dies in the process [1]. Afterwards the fungus grows intercellularly; the only intracellular structures formed are haustoria, the structures responsible for nutrient uptake. Seven to nine days after infection sporulation begins. Of the whole infection process only germination and appressorium formation can be induced in vitro. Under laboratory conditions spores can be germinated in watery suspension or on the water surface [3]. On artificial surfaces, for example polyethylene (PE) membranes, also the formation of appressoria can be induced. Since later structures cannot be observed in vitro it is assumed that at this point a major metabolic switch occurs and that only after this point the biotrophic interaction with the plant begins, which is established after formation of haustoria. To identify genes involved in the biotrophic interaction it is therefore common to compare gene expression in germinated spores or in vitro formed appressoria with gene expression in the infected leaf or in fungal structures isolated from infected leaves, for example haustoria [4,5,6]. An important practical reason to use the in vitro approach is that very early in the infection process the fungal biomass on an inoculated leaf is so small [7], that fungal gene expression cannot be properly quantified.
To elucidate the biotrophic interaction in more detail, genes responsible for nutrient uptake and involved in primary metabolism have been investigated. Results from these experiments support the idea that a major switch in the fungal metabolism occurs after the plant is penetrated [8,9]. A current focus in research of plant pathogens is on the so-called ‘effectors’, or their action. Effectors are defined as proteins secreted by the pathogen and transferred to the host plant to influence host metabolism. Functions of effectors can be suppression of host defense responses or to induce the release of nutrients by the host. If the host is able to recognize such effector proteins by special receptors—resistance proteins—effectors can also have the unintended effect of inducing resistance. In this case they are called avirulence proteins [10,11]. An important part in characterizing effector genes is to study their expression during the infection process.
Since rust fungi cannot be cultivated in vitro, it is almost impossible to genetically modify them. The system for transforming Melampsora lini that depends on selection for host resistance [12] so far is the exception to the rule. In consequence, other ways to influence the pathogen have come to the fore. One way to influence gene expression of plant pathogens is Host Induced Gene Silencing (HIGS). Here, the host plant is made to produce double stranded RNA (dsRNA), complementary to a pathogen gene (the messenger RNA sequence). Either the dsRNA or the downstream siRNA (small interfering RNA) that results from cleavage of the dsRNA by the RNA interference (RNAi) machinery of the host plant, notably the Dicer and Argonaute protein complexes, is transferred to the pathogen and there leads to silencing of the targeted gene. Apart from producing transgenic plants carrying constructs to produce inhibitory RNA, plant viruses can be used to introduce dsRNA [13] into the host plant. Thus, Virus Induced Gene Silencing (VIGS) has become a common tool in plant research. Since in this special case VIGS is used to perform HIGS, we would suggest the term HIGS by VIGS (HbV).
Transcript profiling to measure gene expression is an important part of characterizing pathogen genes, effectors as well as genes involved in metabolism or signaling. The effectiveness of HIGS is also tested by measuring transcript reduction of the target gene. Today the gold standard for measuring transcript levels is RT-qPCR (reverse transcript quantitative PCR).
Transcript levels are normally given relative to a starting sample and they need to be normalized by biomass, total RNA amount, or—most common, time efficient, and most exact—stably expressed genes, so-called reference genes [14]. It is common to use housekeeping genes—genes involved in primary metabolic processes—for normalization, assuming that these important genes need to be expressed unchangingly during every stage of development or in every tissue. However, this assumption is often wrong; many housekeeping genes have been proven to be regulated [15,16]. It is therefore necessary to experimentally determine stably expressed genes that can be used as reference genes. Even more so, since genes that are stable under one set of conditions can be unstable under a different set of conditions it is necessary to determine suitable reference genes for every experiment or every set of samples to be tested [14]. Since stability of expression can only be measured relative to other genes, reference genes are determined by choosing from a set of candidate reference genes comprising housekeeping genes and/or genes that were proven stable in other experiments.
Here, we determined P. pachyrhizi reference genes suitable for normalization in transcript profiling during the infection process including in vitro stages and in measuring transcript levels in HbV experiments. We also determined Glycine max reference genes for normalizing transcript levels of soybean genes during infection with P. pachyrhizi.

2. Results and Discussion

2.1. Choice of Samples

In vitro structures of P. pachyrhizi normally tested in gene expression studies are resting urediospores (s), germ tubes (gt), and appressoria (ap). Therefore, these stages were also included as samples in this study. We also included a wide range of time points of infected leaves. Non-infected leaves (ni), leaves 6 h post inoculation (hpi), 12, 18, 24, 48, 72, 120, 168, and 336 hpi. Time points 6–72 hpi cover the early stages of infection since during this time strong changes in plant gene expression have been observed [17]. Later time points (168 and 336 hpi) represent the sporulation stage and also take into account leaf senescence that is caused by strong P. pachyrhizi infection. Silencing efficiency in our own work using HbV is usually checked at 168 hpi. In order to ensure that our reference genes are not affected by virus infection, we included samples of plants infected with Bean Pod Mottle Virus (BPMV) (21 days post infection (dpi)) and plants infected with both BPMV and P. pachyrhizi (168 hpi).
All samples from leaves were used to determine G. max reference genes. The amount of fungal biomass and therefore the fraction of fungal RNA in infected leaves during early infection, however, is very low. As a consequence Cq values for any fungal gene at these time points are very high and reliable determination of fungal gene expression at these time points is almost impossible. Therefore, the first of the infected leaf samples used for determining P. pachyrhizi reference genes was 72 hpi.

2.2. Choice of Candidate Reference Genes

Housekeeping genes are generally considered good candidates for reference genes. We kept to this rule and to further enhance our chances for choosing stable genes we selected genes that were identified as stable in earlier studies.
For P. pachyrhizi to stay close to the system, we chose from studies dealing with rust fungi. Schmitz et al. [18] worked with soybean rust and mentioned determining reference genes. The other studies deal with Hemileia vastatrix [19], Melampsora larici-populina [20], and Puccinia triticina [21]. In addition five more housekeeping genes were chosen. The actual sequences used were taken from data published by Link et al. [9]. The candidate reference genes are presented in Table 1.
Table 1. Candidate reference genes for P. pachyrhizi.
Table 1. Candidate reference genes for P. pachyrhizi.
Gene DesignationDescriptionGenBank Accession/NameReference
ActactinGACM01001510/TSA:GACM01:Pp_contig05019[18]
ASUBATP synthase β subunitPpGI_Contig363rc a-
aTubα tubulinGACM01003080/TSA:GACM01:Pp_contig00787[20]
CytBcytochrome bGQ332420[18,19]
Elf1aelongation factor 1 αGACM01002155/TSA:GACM01:Pp_contig07365-
Elf3elongation factor 3GACM01003209/TSA:GACM01:Pp_contig00583-
GAPDHglyceraldehyde-3-phosphate dehydrogenaseGACM01002868/TSA:GACM01:Pp_contig00008[20]
PDHpyruvate dehydrogenaseGACM01002893/TSA:GACM01:Pp_contig00239-
PDKpyruvate dehydrogenase kinaseGACM01003388/TSA:GACM01:Pp_contig02726-
RPS940S ribosomal protein S9GACM01002779/TSA:GACM01:Pp_contig00421[19]
RPS1140S ribosomal protein S11GACM01003017/TSA:GACM01:Pp_contig00682[19]
RPS1440S ribosomal protein S14GACM01002939/TSA:GACM01:Pp_contig00153[19]
SDHsuccinate dehydrogenaseGACM01002923/TSA:GACM01:Pp_contig00241[21]
UbcubiquitinGACM01002338/TSA:GACM01:Pp_contig00942[21]
UbcE2ubiquitin conjugated enzymeGACM01000888/TSA:GACM01:Pp_contig06751[21]
a No NCBI accession. See [9] supplementary information for sequence.
For soybean, several studies to determine reference genes have been performed [22,23,24]. One study used G. max homologs of housekeeping genes of Arabidopsis thaliana and validated and ranked them based on stable expression in different organs, under different light regimes, and in different cultivars [22]. Another study [23] searched for new reference genes, using a large dataset of different microarrays to find promising candidates and then validated and ranked these with qRT-PCR data. These authors even included soybean rust infected plant material, however just one sample. A third study [24] compared the classical housekeeping genes [22] to the new ones [23], using an even bigger set of samples. Apart from providing a good source to choose candidate genes (Table 2), these studies also found that no reference gene was stable over all sample sets analyzed and that different genes were the most stable in different sample sets. This once more highlights the necessity for validating reference genes for any given experiment.
Table 2. Candidate reference genes for G. max.
Table 2. Candidate reference genes for G. max.
Gene Designation aDescriptionGene IDReference
cons7metalloproteaseAW310136 b[23]
cons15CDPK-related protein kinaseAW396185 b[23]
CYP2cyclophilinTC224926 c[22,24]
ELF1Belongation factor 1 βTC203623 c[22,24]
SKIP16F-box proteinCD397253 b[23,24]
TIP41TAP42 interacting protein-signalingEV263725 b[24]
Ukn1hypothetical protein/possibly ABC transporterBU578186 b[23,24]
Ukn2hypothetical proteinBE330043 b[24]
a abbreviations as used in the main text; b NCBI accession numbers; c soybase accession numbers.

2.3. Reference Genes for P. pachyrhizi: All Tested Genes Vary in Expression for in Vitro Structures but Are Stable in in planta Samples

Using the method described by Pfaffl [25], we first determined amplification efficiencies of the selected primers. The calculated efficiencies can be found in Table 3.
We then performed RT-qPCR for all samples—in vitro structures and infected leaves starting at 72 hpi. Considering primer efficiencies we first plotted corrected Cq values over time (Figure 1) and then used the algorithms geNorm [26] and NormFinder [27] to determine the most stable genes (Figure 2). To be able to detect differences between different stages in rust development we analyzed the samples from in vitro structures separately from the samples from infected leaf stages and also did an analysis that combined all samples.
In our plots of the Cq values there is a strong variation between samples. Most genes are highly expressed in the germ tube stage and much less in the appressorium stage and in all the in planta stages (Figure 1a). There seems to be considerable co-regulation or rather a generally higher transcriptional activity in the germ tube than in all the other stages.
Figure 1. Variation of Cq values for P. pachyrhizi candidate reference genes between samples (vertical axis: Cq). (a) Average and standard deviation of Cq values for every sample; (bd) box plots representing the range of Cq values; (b) in vitro samples; (c) in planta samples; (d) all samples. For better representation of RNA abundance the Cq values were corrected for primer efficiency: Cqcorr = lg2 ((2 × e)Cq).
Figure 1. Variation of Cq values for P. pachyrhizi candidate reference genes between samples (vertical axis: Cq). (a) Average and standard deviation of Cq values for every sample; (bd) box plots representing the range of Cq values; (b) in vitro samples; (c) in planta samples; (d) all samples. For better representation of RNA abundance the Cq values were corrected for primer efficiency: Cqcorr = lg2 ((2 × e)Cq).
Ijms 16 23057 g001
Table 3. Primers used in qRT-PCR for P. pachyrhizi genes.
Table 3. Primers used in qRT-PCR for P. pachyrhizi genes.
NameGeneSequenceReferenceAmplicon (bp)Efficiency (%)
ActinDis 1 f (KES 1461 fw)Actacagtttcaccacaaccgcc[18]14493
ActinDis 1 r (KES 1462 rv)tgaccgtcgggaagttcg
qPpGI363fwASUBtgtcccaaccctttgctgttthis study14686
qPpGI363rvggcaccgaccatgtaaaacg
AtubDis 1 f (KES 1463 fw)aTubctgcgaacaactatgctcgtc[18]11690
AtubDis 1 r (KES 1464 fw)cacgaagaagccttggagtcc
CytB 1 f (KES 1688 fw)CytBtcaagacgcatccaattctaggtc[28]13885
CytB 1 r (KES 1689 rv)gtgttacacccgtgataatctgaatgat
Elf1a 1 fElf1agtgagcgtggtatcaccatcthis study14392
Elf1a 1 rcagaatggcgcaatcagc
q00583fwElf3aatgcgtggtctctctggtgthis study8995
q00583rvgctcgtccaagatcaccaca
GAPDH 1 f (KES 1465 fw)GAPDHggtatggctttccgagttcca[18]15790
GAPDH 1 r (KES 1466 rv)tcagttgataccaaatcatcctcag
q00239fwPDHgcggaaggaaaggataaggggthis study13695
q00239rvtccgatccttagtctggcct
q02726fwPDKacctcccgttcagctagtctthis study13894
q02726rvaattcatcagagtcggcccc
RibPro 3 fRPS9gtgaatgggagaccaatctcagthis study12894
RibPro 3 rttgcctcctccatgagtcag
q00682fwRPS11ggactgggcttcaagactccthis study9786
q00682rvgaatcctgcccctgatcgag
q00153fwRPS14agttgctcgagtgactggtgthis study8689
q00153rvcatcttgggcagccaacatg
q00241fwSDHcaatcgcctgaggaccgtaathis study8088
q00241rvctggggcaacttgtagagca
Ubc 1 fUbccggaccagtacccttacaaatcthis study12893
Ubc 1 ratcaaacatcggcgaccag
UbcE2 3 fUbcE2gtcgaactgtgacgagtttgthis study11789
UbcE2 3 racggccttagtcttcgatg
That Cq values for the in planta samples are higher than the Cq values for the in vitro samples, especially the germinated spore, is mostly due to the fact, that here fungal RNA is mixed with plant RNA. At the later time points, Cq values in the plant samples decrease; this tendency reflects the increase of fungal biomass in infected leaves.
Since spores are resting structures with no metabolic activity, it is not surprising that transcript levels of most genes are relatively low in this sample. Nevertheless, the fact that there are large differences in transcript levels between spores and germ tubes, also for housekeeping genes, is a highly relevant finding, especially in light of the fact that many studies on gene expression in fungi give expression levels relative to those in spores [20]. What was surprising to us, are the generally high Cq values for the appressorium sample. There is no fundamental reason why expression levels should be lower in appressoria than in germ tubes. One explanation could be the way the appressorium structures are prepared. As described in Section 3.1.2 these structures are grown on polyethylene (PE) sheets for a relatively long time. After 16 h it must be expected that, while some specimens are still growing and forming appressoria, many others are no longer viable, and the RNA within is already degrading. At this stage also contaminations with bacteria or other fungi may be higher than in the germinated spore stage. This in turn would lead to higher Cq values despite the fact that RNA was quantified and tested on a gel—and considered good quality.
Judging from the ranges of Cq values, CytB and SDH are the most stable genes in the in vitro samples (Figure 1b). However, even the range for CytB with 6.5 is still high. For the in planta samples the ranges are generally lower and would be lower still if not for the changes caused by the increase of fungal biomass in the leaf. Here, RPS11 and RPS14 can be rated as most stable (Figure 1c). Over all samples (Figure 1d), again CytB has the smallest range (8.0), followed by Act and RPS14.
In analyses with the algorithms geNorm and NormFinder other genes are considered the most stable (Figure 2). geNorm analysis ranks best RPS14 and UbcE2 for in vitro samples, PDK and CytB for in planta samples, and PDK and Ubc for all samples. Results of the NormFinder analysis are relatively similar to those of the geNorm analysis. With NormFinder RPS14 and UbcE2 are ranked best for in vitro samples, GAPDH and CytB for in planta samples, and PDK and RPS14 for all samples.
However, there are big differences how the genes are ranked between the different combinations of samples. While PDK and RPS14 that are ranked best for all samples by NormFinder are also good in the in vitro and in planta subsets of samples, Ubc, ranked second for all samples by geNorm, is ranked much lower, especially in the in planta subset. Most striking, however, is the difference in ranking for CytB that is ranked among the best for the in planta samples, but worst for the in vitro samples and for the combination of all samples.
The rankings, especially for CytB, are also contrary to what could be deduced from the Cq ranges. While CytB is overall the most stable gene based on its Cq range, it is ranked worst by both algorithms for all samples and for the in vitro samples. Since both algorithms are based on comparing the difference between the Cq values in different samples, they will rank genes that have a similar regulation better than genes that have constant Cq values but differ in their regulation from other genes. This is obviously the case for CytB. CytB is highly expressed but its expression in the germ tube is not as much elevated as that of most of the other genes (Figure 1a). So it seems that the higher stability of CytB is responsible for its lower ranking by geNorm and NormFinder.
It is hard to give a recommendation for what genes should actually be used for normalization for in vitro samples and for all samples. On the one hand CytB, SDH, and Act are the most stable genes, while on the other hand RPS14 and UbcE2 go best with the tendency of all other genes. Therefore, normalizing with either CytB or RPS14 and UbcE2 could be recommended. CytB would be particularly useful, if one would like to get close to absolute quantification with normalization. Using RPS14 and UbcE2 would be more suitable if researchers have the intention to mask the elevated metabolic activity in the germ tube.
Figure 2. Ranking of P. pachyrhizi candidate reference genes for stability by the algorithms geNorm (a,c,e,g) and NormFinder (b,d,f,h). Analyses performed for (a,b) in vitro samples; (c,d) in planta samples; (e,f) all samples; (g,h) in planta samples together with germ tube (gt) sample; M-value: gene-stability measure calculated by geNorm—the arithmetic mean of all pairwise variations of the respective gene against all other genes between the different samples; standard deviation (here): stability measure calculated by NormFinder based on group wise comparisons.
Figure 2. Ranking of P. pachyrhizi candidate reference genes for stability by the algorithms geNorm (a,c,e,g) and NormFinder (b,d,f,h). Analyses performed for (a,b) in vitro samples; (c,d) in planta samples; (e,f) all samples; (g,h) in planta samples together with germ tube (gt) sample; M-value: gene-stability measure calculated by geNorm—the arithmetic mean of all pairwise variations of the respective gene against all other genes between the different samples; standard deviation (here): stability measure calculated by NormFinder based on group wise comparisons.
Ijms 16 23057 g002
The results for the in planta samples are better and much more straightforward. Here we recommend use of CytB in combination with either PDK or GAPDH for normalization.
Expression of potential reference genes is obviously not stable between the three in vitro samples tested; choosing the right reference gene(s) is difficult, and consequently analysis of gene expression during the in vitro stages will be precarious. Following the arguments detailed above the fault lies with the spore and the appressorium samples. For reliable gene expression analyses one should therefore consider to omit these samples. We strongly recommend not to use spores as the sample to which expression in other samples should be compared to, but rather use the germ tube stage as the reference sample. When spore and appressorium are included in expression studies we recommend that expression changes between the in vitro stages should only be considered relevant if they are greater than a thousand fold, since even for the most stable gene tested expression varies roughly a hundredfold.
To identify the most suitable gene(s) for normalizing in a sample set omitting spore and appressorium stages another analysis was performed. This combines in planta samples with germ tube as the only in vitro sample (Figure 2g,h). Here Elf3 and RPS9 or PKD and GAPDH were identified as the most suitable reference genes. Consequently, we recommend to include only in planta stages, and a germ tube sample, and to use a combination of the above genes for normalization in studies that do not necessarily need to include resting urediospores and appressorium stages.

2.4. Reference Genes for G. max: Ukn2 and cons7 Can Be Used for Normalizing RT-qPCR in Tissues Infected with Soybean Rust and Also for Superinfection of Soybean Rust on BPMV

During determination of primer efficiencies (Table 4) the primer pairs chosen for CYP2 (Gm CYP2 f and Gm CYP2 r) and ELF1B (Gm ELF1B f and Gm ELF1B r) gave rise to unspecific side products; consequently, these genes were excluded from further analyses.
After RT-qPCR for the other genes and all samples, the ranges of corrected Cq values suggested cons7, Ukn1, and Ukn2 as the most stable genes (Figure 3). Even though the genes differ strongly in their expression, all appear relatively stable over the tested samples.
Ranking of the genes by geNorm and NormFinder identified Ukn2 and cons7 as the most stable genes (Figure 4). In this case we find no discrepancies between our Cq range analysis and results of the algorithms. We therefore recommend using cons7 and Ukn2 for normalization in expression profiling of soybean genes during soybean rust infection and for quantification of either virus or fungus during HbV experiments.
More studies dealing with reference genes were only found, or were only published, after we finished our experiments [29,30,31,32]. In addition to genes selected from [24] these studies include other frequently used reference genes. Of special interest here is [30] who included infection with Soybean Mosaic Virus in their study and found ELF1B and Ukn2 to be most stable under these conditions. Closest to our own study is [32] which included BPMV infection and here found an ABC transporter and elongation factor 1 α most stably expressed. During infection with powdery mildew, another biotrophic fungus, they found CYP2 and tubulin α-4 most stable. So while this goes well together with our result for Ukn2, this also suggests that there are alternatives to our findings. Also it seems most regrettable in retrospect that CYP2 and ELF1B did not work in our hands.
Figure 3. Variation of Cq values for G. max candidate reference genes between samples (vertical axis: Cq). (a) Average and standard deviation of Cq values for every sample; (b) box plots representing the range of Cq values; For better representation of RNA abundance the Cq values were corrected for primer efficiency: Cqcorr = lg2 ((2 × e)Cq).
Figure 3. Variation of Cq values for G. max candidate reference genes between samples (vertical axis: Cq). (a) Average and standard deviation of Cq values for every sample; (b) box plots representing the range of Cq values; For better representation of RNA abundance the Cq values were corrected for primer efficiency: Cqcorr = lg2 ((2 × e)Cq).
Ijms 16 23057 g003
Figure 4. Ranking of G. max candidate reference genes for stability by the algorithms geNorm (a) and NormFinder (b).
Figure 4. Ranking of G. max candidate reference genes for stability by the algorithms geNorm (a) and NormFinder (b).
Ijms 16 23057 g004

2.5. Number of Reference Genes Needed for Efficient Normalization

Since generally no single gene can be considered truly and absolutely stable, it is also accepted practice to use more than one gene for normalization, expecting that expression changes in the different genes will balance each other. Besides determining the most stable gene, the NormFinder algorithm also calculates a value designated accumulated standard deviation. This value indicates the optimal number of genes to be used for normalization.
Figure 5 shows our results for these analyses in the different sample sets. The optimal numbers differ strongly for different sample sets or combinations of samples. In two cases we have a steep drop of the accumulated standard deviation up to the optimal gene number, in two cases the drop flattens out and the optimal number is reached much later. For the plant genes there is altogether only little alteration between different numbers of genes.
We observe that there are bigger changes between the first few genes whereas differences for higher numbers are less important. We therefore ignore the optimal numbers for in planta samples and all samples (Figure 5b,c), and assume that a lesser number of genes for normalization is also adequate. Of course this is also a question of cost and effort, and here using nine, ten, or five genes, respectively, for normalization is too much. Since in this case the plant genes all have good stability, using fewer genes is justified, and we therefore recommend lower numbers.
Above we are recommending gene combinations for normalization; these combinations always include what we consider the right number of genes.
Figure 5. Number of reference genes recommended by NormFinder (optimal number marked red); (a) for P. pachyrhizi in vitro structures; (b) for P. pachyrhizi in planta structures; (c) for all P. pachyrhizi samples; (d) for all P. pachyrhizi samples excepting spore and appressorium; (e) for G. max.
Figure 5. Number of reference genes recommended by NormFinder (optimal number marked red); (a) for P. pachyrhizi in vitro structures; (b) for P. pachyrhizi in planta structures; (c) for all P. pachyrhizi samples; (d) for all P. pachyrhizi samples excepting spore and appressorium; (e) for G. max.
Ijms 16 23057 g005
Overall the gene combinations we recommend provide a good option for normalization in the studied sample combinations if researchers heed the mentioned caveats, especially limited stability in the in vitro structures. Studies by Vieira et al. [19] and Hacquard et al. [20] suggest using DNA for additional normalization or as a correction for the normalization with reference genes. Even though our results for the in vitro structures seem to suggest that additional normalization is necessary, using DNA for these structures might not be appropriate, since during appressorium formation also nuclear division takes place, so the DNA content increases even though up to this point the biomass stays the same. If researchers need additional normalization we would rather recommend to measure the amount of total RNA used in reverse transcription and take this into account. Only for direct comparisons between resting spore and germ tube, where DNA content stays constant, DNA could be used for normalization. In this case however we would recommend using DNA only, without combining it with RNA references. On the other hand for the in planta structures we find that the tested genes are stable and are reflecting the increase in biomass, so for these samples additional normalization is clearly not necessary.
For the soybean genes, we found that all those tested have good stability, which could be expected, judging from the amount of samples that were already tested in the studies the candidates were chosen from [22,23,24]. Nevertheless, our study now provides evidence that these genes can be used to track expression changes during rust infection and in host-induced gene silencing.

3. Experimental Section

3.1. Plant Material, Fungal Isolate, Inoculation with Virus and Rust Fungus

3.1.1. Plant Material

Soybean (G. max) cultivar Thorne (Bayer CropScience AG, Lyon, France) was grown in the greenhouse with 16 h light/8 h dark at 22 °C. Plants were transplanted after formation of the primary leaves.

3.1.2. Fungal Isolate, Plant Inoculation and in Vitro Structures

Urediospores of isolate Thai 1 (laboratory collection, Phytopathology, University of Hohenheim, Stuttgart, Germany) were first hydrated by stirring vigorously for 45 min in a watery suspension with 0.2% urediospores, 0.2% milk powder, 0.01% Tween20. This suspension was then sprayed on 21 days old G. max plants (in case of superinfection of P. pachyrhizi on BPMV, plants 21 days after virus inoculation (see below)) using a chromatography vaporizer. Plants were then kept in the dark for 12 h at 95% humidity and 20 °C. After that plants were returned to the greenhouse until sampling or for collection of urediospores that appeared starting at 9 dpi.
Germ tubes (gt) were obtained as described by [3] by spreading urediospores on the water surface and incubation for 12 h at room temperature in the dark.
To induce in vitro formation of appressoria (ap), the same suspension used for plant inoculation was sprayed on PE sheets instead and incubated for 16 h at 95% humidity at room temperature in the dark.

3.1.3. BPMV Inoculation

Plasmids pBPMV-IA-R1M and pBPMV-IA-V1 [33] were introduced into soybean plants using particle bombardment as described [34]. 21 dpi symptomatic leaves were harvested, dried and frozen. This infected plant material was homogenized using mortar and pestle, mixed with Sörensen inoculation buffer (1 M NaH2PO4/Na2HPO4, pH 7.2), and with the help of carborundum the mix was rubbed onto freshly developed primary leaves of soybean plants. The carborundum was washed off again and plants were kept in the dark at 95% humidity for 12 to 16 h before they were returned to the greenhouse.

3.2. Sample Homogenization and RNA Isolation

Homogenization was performed using the FastPrep-24 (MP Biochemicals GmbH, Eschwege, Germany) for 20 s at 4.5 m/s.
RNA isolation was performed using the Plant RNA Isolation Mini Kit (Agilent Technologies, Santa Clara, CA, USA). We followed the procedure described in the kit protocol except for adding an additional centrifugation step (1 min, 16,000× g, room temperature) between homogenization and addition of the material to the prefiltration column.
Below are listed particulars for the different sample materials:
Fifty mg urediospores were used directly from storage at −80 °C. Together with two metal beads they were put in a homogenization tube, cooled in liquid nitrogen and homogenized deep frozen. Afterwards 600 μL Extraction Solution was added.
Germ tubes were removed from the water surface by filtering. Together with 600 μL Extraction Solution roughly 100 mg germ tubes were added to a homogenization tube with Lysing Matrix E.
Six hundreds μL Extraction Solution were dispersed on the PE membrane with appressorial structures and the structures scraped together with a rubber and transferred to a homogenization tube with Lysing Matrix E.
One hundred mg of leaf material were excised using a cork borer. Together with two metal beads they were put in a homogenization tube, cooled in liquid nitrogen and homogenized deep frozen. Afterwards 600 μL Extraction Solution were added.
Integrity of the RNA and absence of contaminations were checked by gel electrophoresis using 1% agarose gels with TAE buffer containing 1% bleach (6% NaClO). RNA concentrations were measured using a Qubit® 2.0 fluorometer (Life technologies, Darmstadt, Germany) with the Qubit® Broad Range Assay Kit (Life technologies). Storage of the RNA was as precipitate in 0.3 M sodium acetate (pH 5.3) and 70% ethanol at −80 °C.

3.3. RT-qPCR

3.3.1. Reverse Transcription

Before reverse transcription, DNA contaminations were removed by digestion with PerfeCTa® DNase I (Quanta Biosciences, Gaithersburg, MD, USA) using 5 μg total RNA in a 10 μL reaction.
For reverse transcription the Tetro cDNA Synthesis Kit (Bioline Reagents Ltd., London, UK) was used. The 10 μL from DNase I digestion were directly used in this reaction; we followed the recommendations of the manufacturer using random hexamer primers. cDNA was stored at −80 °C.

3.3.2. Machinery, Chemistry, and Reaction Conditions of Real Time PCR

Real Time PCR reactions were run on a CFX96™ Real-Time-PCR System (Bio-Rad Laboratories, Hercules, CA, USA).
We were using the SensiFAST™ No-ROX mix (Bioline Reagents Ltd., London, UK). Primers were used at a concentration of 0.3 μM; we used 2 μL of cDNA in 20 μL reactions. PCR cycling followed a two step protocol with an initial denaturation step of 95 °C for 5 min, and then 5 s 95 °C and 15 s at 60 °C. 39 cycles were completed; then a melt curve analysis was run: dissociation at 95 °C for 10 s, then the actual melt curve spanning 65 to 95 °C with 0.5 °C temperature increase per 5 s.

3.3.3. Data Analysis

Data logging and the determination of Cq, ΔCq, and ΔΔCq values was done using Bio-Rad CFX Manager 2.1 (Bio-Rad Laboratories). We used automatic threshold determination.
As proposed by [25] we considered primer efficiency in our analysis. Primer efficiencies were determined using rows of four 1:10 dilutions of cDNA derived from germ tubes and 7 dpi infected plant material for P. pachyrhizi primers and also 7 dpi plant material for G. max primers. The analysis was performed in biological triplicate and technical replicate.
Both for determination of primer efficiency and determination of the most stable candidate reference genes, we used the GenEx software package: GenEx 6.0.1.612 (MultiD Analyses AB, Göteborg, Sweden).

3.4. Primers

Many sequences for primers used in this study were taken from previous studies (see references in Table 3 and Table 4). Otherwise we used the webtool Primer3 [35,36] for primer design. We set the following parameters: primer size: minimum 18 bases, optimum 20 bases, maximum 25 bases, primer Tm: min 58 °C, opt 60 °C, max 62 °C; primer GC%: min 30, opt 50, max 70. The product size range was set to 75–150 bp.
The products of the PCR reactions for determination of primer efficiency were also loaded on TAE (tris, acetate, ethylenediaminetetraacetic acid) agarose gels and only when a clear single band was obtained the genes were analyzed further.
Table 4. Primers used in qRT-PCR for G. max genes.
Table 4. Primers used in qRT-PCR for G. max genes.
NameGeneSequenceReferenceAmplicon [bp]Efficiency [%]
Gm cons7 fcons7atgaatgacggttcccatgta[23]11481
Gm cons7 rggcattaaggcagctcactct
Gm cons15 fcons15taaagagcaccatgcctatcc[23]97106
Gm cons15 rtggttatgtgagcagatgcaa
Gm CYP2 fCYP2cgggaccagtgtgcttcttca[22]154na
Gm CYP2 rcccctccactacaaaggctcg
Gm ELF1B fELF1Bgttgaaaagccaggggaca[22]118na
Gm ELF1B rtcttaccccttgagcgtgg
Gm SKIP16 fSKIP16gagcccaagacattgcgagag[24]6087
Gm SKIP16 rcggaagcggagaactgaacc
Gm TIP41 fTip41aggatgaactcgctgataatgg[24]8895
Gm TIP41 rcagaaacgcaacagaagaaacc
Gm UKN1 fUkn1tggtgctgccgctatttactg[24]7490
Gm UKN1 rggtggaaggaactgctaacaatc
Gm UKN2 fUkn2gcctctggatacctgctcaag[24]7996
Gm UKN2 racctcctcctcaaactcctctg

4. Conclusions

Our results show that there are enormous changes in gene expression even among housekeeping genes during early development of soybean rust. There is only little gene expression in resting spores and there may be problems with viability of in vitro structures in the appressorium stage. We therefore recommend to specify gene expression data relative to the germ tube stage rather than the resting urediospore stage. We also suggest that gene expression changes in the appressorium should only be considered relevant when they are a thousand fold or higher. We could identify reference genes that are suitable for expression profiling during the infection process in planta and combined with the germ tube stage: CytB/PDK or CytB/GAPDH and Elf3/RPS9 or PKD/GAPDH are suitable reference gene combinations.
We could show that reference genes Ukn2 and cons7 are most stable in our experimental setup and are best suited for normalization of expression of soybean genes during soybean rust infection and in HbV experiments.

Acknowledgments

We thank the University of Hohenheim, Faculty of Agricultural Sciences for hosting this research. Manuel Müller was supported by a fellowship provided by the Faculty of Agricultural Sciences.

Author Contributions

Daniela Hirschburger, Manuel Müller, Ralf T. Voegele, and Tobias Link planned the research; Daniela Hirschburger and Manuel Müller performed experiments; Daniela Hirschburger, Manuel Müller, and Tobias Link analyzed the data; Tobias Link wrote the paper. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Goellner, K.; Loehrer, M.; Langenbach, C.; Conrath, U.; Koch, E.; Schaffrath, U. Phakopsora pachyrhizi, the causal agent of Asian soybean rust. Mol. Plant Pathol. 2010, 11, 169–177. [Google Scholar] [CrossRef] [PubMed]
  2. Voegele, R.T.; Hahn, M.; Mendgen, K. The Uredinales: Cytology, biochemistry, and molecular biology. In The Mycota V—Plant Relationships, 2nd ed.; Deising, H., Ed.; Springer: Berlin, Germany; Heidelberg, Germany, 2009; pp. 69–98. [Google Scholar]
  3. Posada-Buitrago, M.L.; Frederick, R.D. Expressed sequence tag analysis of the soybean rust pathogen Phakopsora pachyrhizi. Fungal Genet. Biol. 2005, 42, 949–962. [Google Scholar] [CrossRef] [PubMed]
  4. Hahn, M.; Mendgen, K. Characterization of in planta-induced rust genes isolated from a haustorium-specific cDNA library. Mol. Plant Microbe Interact. 1997, 10, 427–437. [Google Scholar] [CrossRef] [PubMed]
  5. Duplessis, S.; Cuomo, C.A.; Lin, Y.C.; Aerts, A.; Tisserant, E.; Veneault-Fourrey, C.; Joly, D.L.; Hacquard, S.; Amselem, J.; Cantarel, B.L.; et al. Obligate biotrophy features unraveled by the genomic analysis of rust fungi. Proc. Natl. Acad. Sci. USA 2011, 108, 9166–9171. [Google Scholar] [CrossRef] [PubMed]
  6. Link, T.I.; Voegele, R.T. Secreted proteins of Uromyces fabae: Similarities and stage specificity. Mol. Plant Pathol. 2008, 9, 59–66. [Google Scholar] [CrossRef] [PubMed]
  7. Voegele, R.T.; Schmid, A. RT real-time PCR-based quantification of Uromyces fabae in planta. FEMS Microbiol. Lett. 2011, 322, 131–137. [Google Scholar] [CrossRef] [PubMed]
  8. Voegele, R.T.; Mendgen, K.W. Nutrient uptake in rust fungi: How sweet is parasitic life? Euphytica 2011, 179, 41–55. [Google Scholar] [CrossRef]
  9. Link, T.I.; Lang, P.; Scheffler, B.E.; Duke, M.V.; Graham, M.A.; Cooper, B.; Tucker, M.L.; van de Mortel, M.; Voegele, R.T.; Mendgen, K.; et al. The haustorial transcriptomes of Uromyces appendiculatus and Phakopsora pachyrhizi and their candidate effector families. Mol. Plant Pathol. 2014, 15, 379–393. [Google Scholar] [CrossRef] [PubMed]
  10. Dodds, P.N.; Rafiqi, M.; Gan, P.H.; Hardham, A.R.; Jones, D.A.; Ellis, J.G. Effectors of biotrophic fungi and oomycetes: Pathogenicity factors and triggers of host resistance. New Phytol. 2009, 183, 993–1000. [Google Scholar] [CrossRef] [PubMed]
  11. Petre, B.; Joly, D.L.; Duplessis, S. Effector proteins of rust fungi. Front. Plant Sci. 2014, 5, 416. [Google Scholar] [CrossRef] [PubMed]
  12. Lawrence, G.J.; Dodds, P.N.; Ellis, J.G. Transformation of the flax rust fungus, Melampsora lini: Selection via silencing of an avirulence gene. Plant J. 2010, 61, 364–369. [Google Scholar] [CrossRef] [PubMed]
  13. Nowara, D.; Gay, A.; Lacomme, C.; Shaw, J.; Ridout, C.; Douchkov, D.; Hensel, G.; Kumlehn, J.; Schweizer, P. HIGS: Host-induced gene silencing in the obligate biotrophic fungal pathogen Blumeria graminis. Plant Cell 2010, 22, 3130–3141. [Google Scholar] [CrossRef] [PubMed]
  14. Bustin, S.A.; Benes, V.; Garson, J.A.; Hellemans, J.; Huggett, J.; Kubista, M.; Mueller, R.; Nolan, T.; Pfaffl, M.W.; Shipley, G.L.; et al. The MIQE guidelines: Minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 2009, 55, 611–622. [Google Scholar] [CrossRef] [PubMed]
  15. Ferguson, B.S.; Nam, H.; Hopkins, R.G.; Morrison, R.F. Impact of reference gene selection for target gene normalization on experimental outcome using real-time qRT-PCR in adipocytes. PLoS ONE 2010, 5, e15208. [Google Scholar] [CrossRef] [PubMed]
  16. Kwon, M.J.; Oh, E.; Lee, S.; Roh, M.R.; Kim, S.E.; Lee, Y.; Choi, Y.L.; In, Y.H.; Park, T.; Koh, S.S.; et al. Identification of novel reference genes using multiplatform expression data and their validation for quantitative gene expression analysis. PLoS ONE 2009, 4, e6162. [Google Scholar] [CrossRef] [PubMed]
  17. Van de Mortel, M.; Recknor, J.C.; Graham, M.A.; Nettleton, D.; Dittman, J.D.; Nelson, R.T.; Godoy, C.V.; Abdelnoor, R.V.; Almeida, A.M.; Baum, T.J.; et al. Distinct biphasic mRNA changes in response to asian soybean rust infection. Mol. Plant Microbe Interact. 2007, 20, 887–899. [Google Scholar] [CrossRef] [PubMed]
  18. Schmitz, H.K.; Medeiros, C.A.; Craig, I.R.; Stammler, G. Sensitivity of Phakopsora pachyrhizi towards quinone-outside-inhibitors and demethylation-inhibitors, and corresponding resistance mechanisms. Pest Manag. Sci. 2014, 70, 378–388. [Google Scholar] [CrossRef] [PubMed]
  19. Vieira, A.; Talhinhas, P.; Loureiro, A.; Duplessis, S.; Fernandez, D.; Silva Mdo, C.; Paulo, O.S.; Azinheira, H.G. Validation of RT-qPCR reference genes for in planta expression studies in Hemileia vastatrix, the causal agent of coffee leaf rust. Fungal Biol. 2011, 115, 891–901. [Google Scholar] [CrossRef] [PubMed]
  20. Hacquard, S.; Veneault-Fourrey, C.; Delaruelle, C.; Frey, P.; Martin, F.; Duplessis, S. Validation of Melampsora larici-populina reference genes for in planta RT-quantitative PCR expression profiling during time-course infection of poplar leaves. Physiol. Mol. Plant Pathol. 2011, 75, 106–112. [Google Scholar] [CrossRef]
  21. Song, X.; Rampitsch, C.; Soltani, B.; Mauthe, W.; Linning, R.; Banks, T.; McCallum, B.; Bakkeren, G. Proteome analysis of wheat leaf rust fungus, Puccinia triticina, infection structures enriched for haustoria. Proteomics 2011, 11, 944–963. [Google Scholar] [CrossRef] [PubMed]
  22. Jian, B.; Liu, B.; Bi, Y.; Hou, W.; Wu, C.; Han, T. Validation of internal control for gene expression study in soybean by quantitative real-time PCR. BMC Mol. Biol. 2008, 9, 59. [Google Scholar] [CrossRef] [PubMed]
  23. Libault, M.; Thibivilliers, S.; Bilgin, D.D.; Radwan, O.; Benitez, M.; Clough, S.J.; Stacey, G. Identification of four soybean reference genes for gene expression normalization. Plant Genome 2008, 1, 44–54. [Google Scholar] [CrossRef]
  24. Hu, R.; Fan, C.; Li, H.; Zhang, Q.; Fu, Y.F. Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR. BMC Mol. Biol. 2009, 10, 93. [Google Scholar] [CrossRef] [PubMed]
  25. Pfaffl, M.W. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acid Res. 2001, 29, e45. [Google Scholar] [CrossRef] [PubMed]
  26. Vandesompele, J.; de Preter, K.; Pattyn, F.; Poppe, B.; van Roy, N.; de Paepe, A.; Speleman, F. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3, research0034.1–research0034.11. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Andersen, C.L.; Jensen, J.L.; Orntoft, T.F. Normalization of real-time quantitative reverse transcription-PCR data: A model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64, 5245–5250. [Google Scholar] [CrossRef] [PubMed]
  28. Schmitz, H.K. In Vivo und molekularbiologische Untersuchungen zur Sensitivität von Phakopsora pachyrhizi gegenüber Demethylierungs-Inhibitoren und Qo-Inhibitoren; Universität Hohenheim: Stuttgart, Germany, 2013; p. 41. (In Germany) [Google Scholar]
  29. Le, D.T.; Aldrich, D.L.; Valliyodan, B.; Watanabe, Y.; Ha, C.V.; Nishiyama, R.; Guttikonda, S.K.; Quach, T.N.; Gutierrez-Gonzalez, J.J.; Tran, L.S.; et al. Evaluation of candidate reference genes for normalization of quantitative RT-PCR in soybean tissues under various abiotic stress conditions. PLoS ONE 2012, 7, e46487. [Google Scholar] [CrossRef] [PubMed]
  30. Ma, S.; Niu, H.; Liu, C.; Zhang, J.; Hou, C.; Wang, D. Expression stabilities of candidate reference genes for RT-qPCR under different stress conditions in soybean. PLoS ONE 2013, 8, e75271. [Google Scholar] [CrossRef] [PubMed]
  31. Miranda, V.D.J.; Coelho, R.R.; Viana, A.A.; de Oliveira Neto, O.B.; Carneiro, R.M.; Rocha, T.L.; de Sa, M.F.; Fragoso, R.R. Validation of reference genes aiming accurate normalization of qPCR data in soybean upon nematode parasitism and insect attack. BMC Res. Notes 2013, 6, 196. [Google Scholar] [CrossRef] [PubMed]
  32. Bansal, R.; Mittapelly, P.; Cassone, B.J.; Mamidala, P.; Redinbaugh, M.G.; Michel, A. Recommended reference genes for quantitative PCR analysis in soybean have variable stabilities during diverse biotic stresses. PLoS ONE 2015, 10, e0134890. [Google Scholar] [CrossRef] [PubMed]
  33. Zhang, C.; Bradshaw, J.D.; Whitham, S.A.; Hill, J.H. The development of an efficient multipurpose bean pod mottle virus viral vector set for foreign gene expression and RNA silencing. Plant Physiol. 2010, 153, 52–65. [Google Scholar] [CrossRef] [PubMed]
  34. Zhang, C.; Yang, C.; Whitham, S.A.; Hill, J.H. Development and use of an efficient DNA-based viral gene silencing vector for soybean. Mol. Plant Microbe Interact. 2009, 22, 123–131. [Google Scholar] [CrossRef] [PubMed]
  35. Untergasser, A.; Cutcutache, I.; Koressaar, T.; Ye, J.; Faircloth, B.C.; Remm, M.; Rozen, S.G. Primer3—New capabilities and interfaces. Nucleic Acids Res. 2012, 40, e115. [Google Scholar] [CrossRef] [PubMed]
  36. Koressaar, T.; Remm, M. Enhancements and modifications of primer design program Primer3. Bioinformatics 2007, 23, 1289–1291. [Google Scholar] [CrossRef] [PubMed]

Share and Cite

MDPI and ACS Style

Hirschburger, D.; Müller, M.; Voegele, R.T.; Link, T. Reference Genes in the Pathosystem Phakopsora pachyrhizi/ Soybean Suitable for Normalization in Transcript Profiling. Int. J. Mol. Sci. 2015, 16, 23057-23075. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms160923057

AMA Style

Hirschburger D, Müller M, Voegele RT, Link T. Reference Genes in the Pathosystem Phakopsora pachyrhizi/ Soybean Suitable for Normalization in Transcript Profiling. International Journal of Molecular Sciences. 2015; 16(9):23057-23075. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms160923057

Chicago/Turabian Style

Hirschburger, Daniela, Manuel Müller, Ralf T. Voegele, and Tobias Link. 2015. "Reference Genes in the Pathosystem Phakopsora pachyrhizi/ Soybean Suitable for Normalization in Transcript Profiling" International Journal of Molecular Sciences 16, no. 9: 23057-23075. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms160923057

Article Metrics

Back to TopTop