Skip to main content
Download PDF

Genome-wide association study, combined with bulk segregant analysis, identify plant receptors and defense related genes as candidate genes for downy mildew resistance in quinoa

BMC Plant Biology volume 24, Article number: 594 (2024) Cite this article

Abstract

Background

Downy mildew is the most relevant disease of quinoa and the most widespread. Though, little is known about the genetics of resistance to this disease. The objective of this study was to identify the genomic regions controlling downy mildew resistance in quinoa and candidate genes for this trait. With this aim we carried out a GWAS analysis in a collection formed by 211 quinoa accessions from different origins. This approach was combined with inheritance studies and Bulk Segregant Analysis (BSA) in a segregating population.

Results

GWAS analysis identified 26 genomic regions associated with the trait. Inheritance studies in a F2 population segregating for resistance revealed the existence of a major single dominant gene controlling downy mildew complete resistance in quinoa accession PI614911. Through BSA, this gene was found to be located in chromosome 4, in a region also identified by GWAS. Furthermore, several plant receptors and resistance genes were found to be located into the genomic regions identified by GWAS and are postulated as candidate genes for resistance.

Conclusions

Until now, little was known about the genetic control of downy mildew resistance in quinoa. A previous inheritance study suggested that resistance to this disease was a quantitative polygenic trait and previous GWAS analyses were unable to identify accurate markers for this disease. In our study we demonstrate the existence of, at least, one major gene conferring resistance to this disease, identify the genomic regions involved in the trait and provide plausible candidate genes involved in defense. Therefore, this study significantly increases our knowledge about the genetics of downy mildew resistance and provides relevant information for breeding for this important trait.

Peer Review reports

Background

Quinoa (Chenopodium quinoa Willd.) is a member of the Amaranthaceae family. Quinoa is a predominantly autogamous (self-pollinated) species with varying rates of natural hybridization (10–17%) [1]. It is an allotetraploid (2n = 4x = 36) but shows disomic inheritance for most qualitative traits [2]. Quinoa was initially domesticated by the indigenous civilizations of Bolivian and Peruvian Altiplano [3] and, subsequently, the crop has expanded to western South America. Two germplasm pools have been reported in quinoa: Andean highland quinoa, which is the primary center of diversity, and central and southern Chilean quinoa, the second center of diversity [4]. Through a process of selection and diversification, the species is now divided in five major ecotypes: Altiplano, Salar, Yunga, Valley and Lowland [5].

The exceptional nutritional characteristics of quinoa, coupled with its intrinsic tolerance to drought, salinity and frost has attracted worldwide attention to quinoa cultivation [6]. Quinoa provides all the essential amino acids required for humans [7,8,9,10,11,12], being also rich in minerals, vitamins, dietary fiber, linolenate, and natural antioxidants. For all these nutritional qualities quinoa is considered a “superfood” and its consumption has increased in the last years.

Quinoa remains an important food crop in South America, but, these desirable characteristics of quinoa have led to its cultivation expanding to numerous countries, being currently grown in more than 95 countries [5]. In Spain, quinoa cultivation started around 10 years ago, being now an emerging crop with about 6000 ha planted.

For a sustainable cultivation of quinoa, controlling quinoa diseases through an environment-friendly method, as genetic resistance, is desirable. This is especially relevant in this crop because consumers demand mainly organic quinoa. The main disease affecting quinoa worldwide is downy mildew, caused by the biotrophic oomycete Peronospora variabilis Gäum. Therefore, resistance to this pathogen is a key breeding target. P. variabilis infects the leaves of the plant. The initial symptoms are small, isolated chlorotic spots on the upper face of the leaves that later grow into irregular chlorotic spots, that finally become necrotic. On the underside of the leaves, the sporulation of the pathogen in the lesions produces the appearance of a greyish or purplish layer. In cases of severe infection, defoliation occurs. This disease can cause up to 99% yield losses in susceptible cultivars [13]. In Spain, this disease affects severely this crop, affecting up to 90% of the plant area and causing defoliation in susceptible cultivars under conditions especially favourable for the disease.

Downy mildew resistance in quinoa ranges from complete resistance to high susceptibility [13,14,15,16]. These observations suggest that resistance to P. variabilis on quinoa could be controlled by both major and minor genes, depending on the accession, but little is known about the inheritance of the trait. Benlhabib et al. [17] evaluated several traits, including, resistance to downy mildew, in a F2:6 quinoa population derived from a cross between the slightly susceptible accession NL-6 and the resistant accession 0654. Their results suggest that resistance in these lines is a polygenic trait, as around 50% of F2:6 families were classified between the two parents and transgressive segregation for resistance was observed, indicating that resistance could be controlled by different genes in the parental lines. The genetics of downy mildew resistance in quinoa has also been analysed in germplasm collections in two studies using Genome-Wide Association Studies (GWAS) [18, 19]. However, these studies were unable to identify markers associated with the trait, or the markers identified were not consistent. Therefore, further studies are needed to unravel the genetic structure of downy mildew resistance in quinoa and to identify the genes controlling the trait. The identification of these genes, and molecular markers linked to them, would facilitate their introduction and combination into susceptible varieties.

The aim of the present study was to unravel the genetics of resistance to this important disease in quinoa and identify candidate genes for the trait. With this aim, in this study we combined a GWAS, carried out in a germplasm collection formed by 211 quinoa accessions, with inheritance studies and a Bulk Segregant Analysis performed using a cross segregating for resistance. GWAS identified several genomic regions associated with downy mildew resistance in quinoa and a set of plant receptors and defense related genes, located into these regions, are postulated as candidate genes for this trait. Furthermore, resistance in accession PI614911 was found to be controlled by a single dominant gene that is located in chromosome 4, into a region also identified by GWAS.

Methods

Evaluation of the response to downy mildew in a quinoa germplasm collection

The response to downy mildew was scored under field conditions in a collection formed by 211 quinoa accessions with different geographical origins and covering the two main germplasm pools described in quinoa. Cutivars F16 and Kancolla were also included in the assays. The code and origin of these accessions is included in Supplementary Table S1. Seeds were obtained from the USDA North Central Regional Plant Introduction Station of the US National Plant Germplasm SystemUSDA (EEUU) and IPK Gatersleben (Germany) genebanks, excepting cultivar F16, that was provided by Algosur S.A. company. To ensure homogeneity of accessions, before performing the sequencing and field experiments, one plant from each accession was selected and selfed at least twice. An initial set of 138 accessions were screened in 2019 at experimental plots located in Córdoba (37°53′4.226″N 4°46′46.443″W) (Spain) and the whole collection was screened during 2021 and 2022 in Córdoba and Guadajira (38°51’07’’ N, 6°40’49’’ O) (Spain). Accessions were sown according to a completely randomised block design with three blocks. In each block each accession was represented by a 1 m row (10 plants per row) separated 0.7 m from the other rows. Basal fertilization (400 kg/ha of 8:15:15 N: P:K fertilizer plus 87 Kg of urea/ha) was applied before sowing and top dressing (130 kg urea/ha) at flowering.

The severity of the disease, estimated as the percentage of the plant’s leaf area with symptoms, was scored using the “three-leaf screening method”, which considers the average percentage of leaf area in each plant that is infected by the pathogen in three leaves randomly selected: one from each of the lower, middle and upper part of the plant [20]. Disease severity was evaluated once a week, from the time the first symptoms appeared until the senescence of the plant made it difficult to distinguish the symptoms caused by downy mildew from those caused by senescence. Disease severity in the last assessment was considered as final disease severity and used in the analyses.

The correlation between the severities obtained in the different environments was calculated using Pearson’s correlation coefficient.

Genome sequencing and identification of genomic variations

For each quinoa accession forming the collection, DNA was isolated from frozen young leaf tissue obtained from plants grown in a greenhouse using “NucleoSpin Plant II” (Macherey-Nagel GmbH, Germany) kit. After checking its purity and quality by agarose gel electrophoresis, DNA concentration was determined using a Qubit instrument and optimum DNA samples were sent to Diversity Array Technology Pty Ltd (Camberra, Australia) for sequencing and genotyping as described in [21]. Briefly, DNA samples were processed as follows: PstI and MseI compatible adaptors with two different restriction enzyme overhangs were added [22]. The PstI and MseI compatible adaptors were designed to include the Illumina flowcell attachment sequence, sequencing primer sequence and “staggered”, varying length barcode region, similar to the sequence reported by [23]. The reverse adaptor contained the flowcell attachment region and MseI compatible overhang sequence. Only “mixed fragments” were effectively amplified in 30 rounds of PCR using the following reaction conditions: 94 °C for 1 min; 30 cycles of: 94 °C for 20 s, 58 °C for 30 s, 72 °C for 45 s; 72 °C for 7 min. After PCR, equimolar amounts of amplified product were bulked and subjected to 100 cycles of sequencing (single reads) on the Illumina Illumina NovaSeq sequencer. Sequences generated from each lane were processed using proprietary DArT analytical pipelines. In the initial pipeline, poor quality sequences were removed, with more stringent filtering parameters applied to the barcode region compared to the rest of the sequence, ensuring the assignments of the sequences to specific samples (based on the “barcode split”) was reliable. Filtering was performed on the raw sequences using the following parameters: Barcode region minimum Phred score 30, minimum percentage 75; whole read minimun Phred score 10, minimum percentage 50. Approximately 340,412 unique sequences per sample were used in marker calling. Identical sequences were collapsed into “fastqcoll files” which were “groomed” using DArT PL’s proprietary algorithm which corrects low quality base from singleton tag into a correct base using collapsed tags with multiple members as a template. The “groomed” fastqcol files were used in the secondary pipeline for DArT PL’s proprietary SNP and SilicoDArT (presence/absence of restriction fragments in representation) calling algorithms (DArTsoft14). For SNP calling, all tags from all libraries included in the DArTsoft14 analysis are clustered using DArT PL’s C + + algorithm at the threshold distance of 3, followed by parsing of the clusters into separate SNP loci using a range of parameters, including the balance of read counts for the allelic pairs. Additional selection criteria were added to the algorithm based on analysis of approximately 1,000 controlled cross populations. Testing for Mendelian distribution of alleles in these populations facilitated selection of technical parameters discriminating true allelic variants from paralogous sequences. In addition, multiple samples were processed from DNA to allelic calls as technical replicates and scoring consistency was used as the main selection criteria for high quality/low error rate markers. Markers identified (SNPs and SilicoDArT) were assigned to chromosomes using version one of quinoa reference genome (CoGe id33827).

Population structure analysis

The software fastSTRUCTURE [24] was used to estimate the number of populations (K) represented in the data. The input used was a reduced set of SNPs with linkage disequilibrium r² < 0.2 computed on a window of fifty markers that shifts five at the end of each step. PLINK 1.9 [25] was used for this filtering.

In addition, population structure was also inferred by Principal Component Analysis using GAPIT 3.1.0.

To test whether there were significant differences in ‘final disease severity’ between the different populations predicted by PCA and fastSTRUCTURE two analyses of variance (one for each software) were performed. In addition, in the case of the populations predicted by fastSTRUCTURE, as there were more than two populations, comparisons of mean values were performed by least significant difference (LSD) test. These analyses were carried out using the Statistix 8.0 package (Analytical Software, Tallahassee, FL, USA).

Linkage disequilibrium (LD) analysis

LD and squared correlation coefficients (r²) between SNPs within a sliding window of fifty SNPs were computed using TASSEL 5.0 [26]. LD decay and LD half decay distance were estimated using Hill and Weir [27, 28] formula in R [29].

GWAS analysis

First, in order to obtain accurate results, markers that could not be assigned to any chromosome and those that showed more than 20% missing values were excluded. GWAS analyses were carried out by GAPIT 3.1.0 and TASSEL 5.0 software. To take population structure and kinship into account, TASSEL 5.0. was used to obtain a Q matrix (using the multidimensional scaling (MDS) method) and a kinship matrix, that were subsequently used in GAPIT. Several models: MLM (Mixed Linear Model), GLM (General Linear Model), MLMM (Multiple Loci Mixed Model), FarmCPU and BLINK were used. MLM model was analysed using Tassel software while the rest of models were analysed using GAPIT software. “Model Selection” tool, as implemented by GAPIT, was used to determine the optimal number of PCs (covariates) to include for each phenotype. Multiple testing was corrected using Benjamini and Hochberg [30] false discovery rate (FDR) [31] (q = 0.1). Furthermore, quantile-quantile plots (QQplots) were obtained for each model and only when the QQplots showed that the data fit the model, the resulting marker-trait associations (MTAs) were considered accurate. In QQplots, the observed - log10 (P) for each marker are plotted against expected - log10 (P) values under the null hypothesis (no association of the markers with the trait). It is expected that only a few markers would be associated with the trait that is being evaluated. Therefore, if a model is suitable for analysing the data, in the QQplot most of the markers should be on or near the middle line between the x-axis and the y-axis and only a few (those associated with the trait) will be far from this middle line. In addition, MTAs were only considered reliable when their allele frequency was > 5%.

An analysis of variance was performed to check the effect of the factors ‘accession’, ‘location’ and ‘year’ on the variable ‘disease severity’. This analysis was performed using the software IBM SPSS statistic (version 26). As this analysis indicated that all factors, and their interactions, were significant (Suppl. Table S2), “disease severity” scored in each combination ‘location * year’ was considered a different trait and analysed separately in the GWAS analyses.

Identification of candidate genes controlling resistance to downy mildew

To identify candidate genes controlling resistance to downy mildew in the quinoa collection, we searched in the version one of quinoa reference genome of cultivar QQ74 (NCBI code ASM168347v1) in a range of 250 kb down and upstream the MTAs identified by GWAS, using the browse tool available at GeGo website (https://genomevolution.org/CoGe/SearchResults.pl?s=quinoa&p=genome; CoGe id33827). This threshold was selected because it was the maximum range for LD half decay distance calculated across the different chromosomes in our study.

Inheritance studies

The inheritance of resistance to P. variabilis was studied in an F2 population derived from the cross between the resistance accession PI614911 and the susceptible breeding accession Q122. The cross was made according to the method described by [32]. Q122 was used a female parent and PI614911 as male parent. In order to confirm that the seeds obtained were real F1, and not the result of self-pollination, a set of RAPD (Random Amplified Polymorphic DNA) markers were surveyed in the parental lines. DNA extraction and RAPD analyses were performed as reported in [33]. A primer OPC16 (Operon Technologies, Alameda California), showing polymorphism between the parents was next tested individually on DNA from the different F1 plants obtained. One F1 plant, showing bands from both parents, being, therefore the result of a real cross, was selfed in a greenhouse to obtain the F2 population. F2 plants were sown in the field and selfed to obtain the F3 families.

Resistance to P. variabilis in the Q122 x PI614911 F2 population and parental lines was evaluated in 2019 in an experimental plot located in Córdoba (Spain). Parental lines were sown in three replicates, having each 1 m row of each parent with 10 plants per row. To evaluate the F2 population, two hundred F2 seeds were sown in a row. F3 families were sown in rows having each ninety seeds of the corresponding family and evaluated during 2020 season. Plants were evaluated several times as described above and classified as resistant or susceptible according to their disease severity.

Goodness of fit to expected segregations was checked using chi-square tests.

Bulk segregant analysis

According to the results obtained in the evaluation of downy mildew resistance in the F2 and F3 generations derived from cross Q122 x PI614911, seven F2 plants homozygous for resistance and ten plants homozygous for susceptibility were selected. Their DNA was extracted, as described above for the GWAS panel, and two pools, one formed by the resistance plants an another formed by the susceptible plants, were created mixing equal amount of DNA from each of the plants forming the pool. These pooled DNA samples were sent to Diversity Array Technology Pty Ltd (Camberra, Australia) for sequencing and SNP calling, as described above for the GWAS collection, but, in this case, each sample was sequenced twice. Those SNPs that could not be assigned to chromosomes were omitted and read depth information was used for BSA analysis using BSAvis software (https://github.com/FadyMohareb/BSAvis_GP_2020/tree/main/BSAvis). Briefly, a SNP index was calculated across the different chromosomes [34], those SNPs having a SNP index in both pools < 0.3 or equal to 1 were excluded, and an average SNP index was calculated using the sliding window method (window size of 1 Mb and a step size of 10 kb). ΔSNP index was further calculated as the difference between the SNP index of the two pools and a ΔSNP-index graphs was generated by plotting ΔSNP-index against the position in each chromosome.

The hypothesis is that, for a marker unlinked to the resistance gene, having 50% mutant and 50% wild-type sequence reads is expected, while the causal SNP, and closely linked SNPs, should show 100% mutant and 0% wild-type reads. SNPs loosely linked to the causal mutation should have > 50% mutant and < 50% wild-type reads. If we define the SNP index as the ratio between the number of reads of a mutant SNP and the total number of reads corresponding to the SNP, we expect that this index would be equal 1 near the causal gene and 0.5 for the unlinked loci [35] We further calculated the difference between the SNP index of the two pools to obtain the ΔSNP index. Δ SNP-index would be equal to 1 when the genome of bulked DNA is consistent with that of one of the parents, while Δ SNP-index will be − 1 when the genome of bulked DNA is consistent with that of the other parent, and Δ SNP-index = 0 when both parents had the same SNP-indices at the genomic regions. Thus, the Δ SNP-index value should be different from 0 if a genomic region harbours a target gene.

Results

Genotyping

The collection formed by 211 quinoa accessions was sequenced and genotyped using DArTseq technology. DArTseq is a genome complexity reduction-based sequencing technology (https://www.diversityarrays.com/services/dartseq/) that produces two types of markers: SNPs (Single Nucleotide Polymorphism) and SilicoDArT (presence/absence of the tags sequences). An average of 2,564,206 reads were obtained per DNA sample and used for marker calling. Marker calling quality was validated by high average read depth per locus (average across all markers was over 15.7 reads/locus). After eliminating those markers with low quality, those that could not be assigned to chromosomes and those having more than 20% missing values, a total of 12,397 SNPs and 12,720 SilicoDArT markers were selected and used in further analyses. The distribution of markers and the average distance between markers, per chromosome, is shown in Suppl. Table S3.

Linkage disequilibrium analysis

LD half decay distance was estimated for each chromosome. This parameter varied between chromosomes ranging between 67,209 bp (in chromosome 13) and 248,075 bp (in chromosome 6) and showed an average value of 126,448 bp (Suppl. Fig. S1).

Population structure

Population structure in the quinoa collection was examined using Principal Components Analysis (PCA) and fastSTRUCTURE software. PCA, using the first and second principal components (that explained 25 and 11% of the variation, respectively) divided the collection into two main groups, in agreement with the two main germplasm pools reported for quinoa (Fig. 1A). Thus, one group was formed mainly by highland quinoa accessions (accessions mainly from Peru and Bolivia) and the other formed mainly by lowland quinoas accessions (accessions mainly from Chile and USA). Accessions from USA had been previously reported to be closely related with accessions from Chile as the USDA germplasm had been collected at these geographical regions [19].

fastSTRUCTURE software divided these two groups in additional subpopulations. So that, according to fastSTRUCTURE and the “chooseK.py” script included in it, that provides the value of K that best fits the data, the collection could be divided in four populations (Fig. 1B). One of this population corresponded to the PCA group containing the lowland accessions, while highland accessions were subdivided in three groups, one formed mainly by accessions from Bolivia, other formed mainly by accessions from Peru and a third one containing accessions from both Bolivia and Peru. These results show, in agreement with previous studies [19], that highland quinoa accessions show a wider genetic diversity than lowland accessions.

The assignment of the accessions included in the GWAS panel to each predicted population group can be seen in Suppl. Table S1.

Fig. 1
figure 1

Population structure of the quinoa germplasm collection according to Principal Component Analysis. (A) and according to fastSTRUCTURE software (B)

Full size image

Response to downy mildew in a quinoa germplasm collection

Response to downy mildew was scored in a germplasm collection formed by 211 quinoa accessions during 2019, 2021 and 2022 in Córdoba (South of Spain) and during 2021 and 2022 in Guadajira (Central West Spain). The collection showed substantial phenotypic variation for this trait (Fig. 2). In all seasons disease severity showed a continuous distribution ranging from high resistance to high susceptibility, although disease was more severe in Córdoba than in Guadajira (Suppl. Fig. S2). So that, in Córdoba disease severity ranged from 0 to 82.5% of the plant area affected by the disease in 2019, from 2 to 73.3% in 2021 and from 0 to 73.3% in 2022. In Guadajira, disease severity ranged from 10 to 55% in 2021 and from 5 to 31.6% in 2022. These data show that downy mildew can severely affect the quinoa crop in Spain when the accession is susceptible and the conditions conductive for the development of the disease, but that, fortunately, genetic resistance to the P. variabilis isolates present in Spain is available within quinoa germplasm. Disease severity values were not correlated between locations but showed a significant correlation between the scorings carried out in the same location (Table 1).

Fig. 2
figure 2

Quinoa accessions resistant to downy mildew (A) and susceptible to downy mildew (B)

Full size image
Table 1 Pearson’s correlation coefficient between downy mildew disease severity values scored in a collection of 211 quinoa accessions in different environments
Full size table

GWAS and candidate genes analyses

GWAS analysis was performed using MLM, GLM, MLMM, FarmCPU and BLINK models. Resulting QQplots indicated that, in general, MLMM, FarmCPU and BLINK were a better fit to the data than MLM and GLM (Suppl. Fig. S3). A total of 58 MTAs, corresponding to 26 genomic regions, passes the quality criteria (good fit to the model, according to the QQplot, p-adjusted < 0.1 and MAF > 0.05) and, were, therefore, considered to be associated with resistance to downy mildew accurately in the quinoa germplasm collection. Several of these regions were identified in more than one environment, while others were specific for a certain environment. These regions are summarized in Table 2 and the exact location of each MTA can be found in Supplementary Table S4.

We then searched for candidate genes with a putative function in defense within 250 kb down and upstream the MTAs identified. Interestingly, most of the genomic regions surrounding the MTAs contained plant receptor or defense related genes, such as “receptor like proteins”, “disease resistance proteins”, “Wall-associated receptor kinases”, “pathogenesis like protein”, “Zinc finger BED domain-containing protein” and “L-type lectin-domain containing receptor kinase”, among others (Table 2; Suppl. Table S4). Remarkably, for ten of the MTAs, these candidate genes with a putative function in defense were located exactly in the same genomic position as the MTA (not in the surroundings 250 kb) (Suppl. Table S4).

Table 2 Genomic regions associated with downy mildew resistance in a quinoa germplasm collection identified by GWAS analysis. For each region, position on quinoa genome, the environment where it was identified, and candidate genes located into the region are shown
Full size table

Inheritance studies

In 2019, the parental line Q122 was susceptible to P. variabilis, showing at the end of the disease assessment period, on average, 28.3% of its area covered by the disease and sporulation. By contrast, plants of the resistant parent PI614911 were indeed highly resistant, showing no symptoms or, at maximum, a few scattered yellow spots caused by the disease and no sporulation. The F2 population Q122 x PI614911 segregated into 92 resistant and 30 susceptible, fitting perfectly the 3:1 ratio expected for a single dominant gene (χ2 = 0.01; p = 0.92) (Table 3). Differences between resistant and susceptible plants were evident as resistant plants showed as maximum a few scattered spots caused by the disease while susceptible plants were at least as susceptible as their susceptible parent. The phenotype of the F2 plants was confirmed evaluating their derived F3 families during 2020. All F2 resistant plants produced F3 families that were resistant or segregated for resistance, while all F3 families derived from F2 susceptible plants were susceptible. However, due to severe problems with the emergence of the seeds, for most families less than 10 F3 plants emerged. In order to obtain accurate conclusions, we only considered for the resistance segregation analysis those families having at least 10 plants. Of the 56 F3 families having at least 10 plants, 12 were resistant, 12 were susceptible and 32 segregated for resistance. These ratios also fit to the 1:2:1 ratio expected for a single dominant gene (χ2 = 1.14; p = 0.56). The ratio of resistant plants to susceptible plants in the 32 families that segregated for resistance was consistent with the hypothesis of the existence of a major dominant gene controlling resistance.

Bulk segregant analysis (BSA)

The segregating F2 population derived from the cross between the susceptible line Q122 and the resistant PI614911, described above, was used to perform a BSA. After eliminating those SNPs that could not be assigned to chromosomes, a total of 11,418 SNPs were included in the analysis. ΔSNP-index graph across the different chromosomes identified a region showing an increase in ΔSNP-index values in the region 10.3–16.4 Mb on chromosome 4 (version 1 of quinoa reference genome CoGe id33827). In this region average ΔSNP index was high, being > 0.7 for 13 SNPs. This profile was not observed in any other region of the genome. Taking into account that in BSA, for a F2 population, the ΔSNP-index threshold to consider an imbalance of allele frequencies is 0.67 [36], these results demonstrate that the dominant gene controlling resistance to downy mildew in accession PI614911 is located into this region. SNPs identified by DArTseq sequencing were further positioned on the version two of quinoa reference genome (CoGe id607169) and a ΔSNP-index graph was also created. In agreement with the previous results, region 38.64–42.51 Mb on chromosome Cq2A, that corresponds to the region identified as a carrier of the resistance gene in version one of quinoa reference genome, also showed a clear increase in ΔSNP index (Fig. 3). Furthermore, this region is included into a region identified in the GWAS analysis performed in our study, confirming both, GWAS and BSA studies.

Table 3 Segregation of resistance to P. variabilis in Q122* PI614911 cross
Full size table
Fig. 3
figure 3

Genomic region containing the major gene conferring resistance to downy mildew in accession PI614911 identified by (A) BSA using version one of quinoa reference genome CoGe id33827 (B) BSA using version two of quinoa reference genome (CoGe id60716) (C) GWAS analysis performed using GAPIT software, SilicoDart markers and BLINK model. Manhattan plot for chromosome 4 (version one of quinoa reference genome) is shown. Chromosome 4 in version one of the quinoa reference genome corresponds to chromosome Cq2A in version two of the quinoa reference genome, and the region highlighted in (A) corresponds to the same region highlighted in (B)

Full size image

Discussion

Despite the relevance of downy mildew disease in quinoa little is known about the genetics of downy mildew resistance in quinoa. Knowledge of the genetic control of this trait would be useful to plan the best strategy to incorporate this trait into elite cultivars. The only study analysing the inheritance of P. variabilis resistance was performed by [17]. The authors evaluated the response to this disease in a recombinant inbred line population derived from a cross between a slightly susceptible accession and a resistant accession. In this population the trait behaved as a quantitative trait and transgressive segregation was observed, suggesting that resistance was polygenic and that the parents harboured different resistant genes. By contrast, in our study we have identified complete qualitative downy mildew resistance and demonstrated that this resistance is controlled by a single dominant gene in quinoa accession PI614911. Single gene-controlled resistance is easy to incorporate into susceptible material through backcrossing. Therefore, the identification of a major gene conferring complete downy mildew resistance in quinoa is a milestone that will greatly facilitate the development of resistant cultivars. So that, we have already used accession PI614911 to successfully incorporate resistance to downy mildew in some of our more interesting advanced breeding material.

Furthermore, through a BSA analysis, we have identified the genomic region containing this major gene. The resolution of BSA is expected to be lower than GWAS because the number of generations over which the population is interbred is limited [37]. However, BSA was useful to corroborate GWAS results and discern which of the different genomic regions identified by GWAS was responsible for resistance in accession PI614911. Our results reveal that the gene conferring resistance in accession PI614911 is located in the region 10.3–16.4 Mb on chromosome 4 (version one of quinoa reference genome). This region corresponds to the region 38.9–42.8 Mb on Chromosome Cq2A in version 2 of the quinoa reference genome (CoGe id60716) (Fig. 3), a region that was also postulated to be associated with resistance to downy mildew in the GWAS analysis performed by [19]. These results suggest that this resistance gene is not a rare gene but, rather, a gene that may be present in several quinoa accessions. Exciting, a gene predicted to be similar to “Disease resistance protein RGA2” is located exactly in the same positions as the MTA identified at 16,082,246 pb on Chromosome 4. This gene is an excellent candidate for this major resistance gene. Further sequencing of this candidate gene and gene expression studies, in PI614911and susceptible lines, and mapping in segregating populations are planned to confirm this hypothesis. Other interesting candidate genes located in this region are different types of “Receptor-like protein kinases” and a “Zinc finger BED domain-containing protein” (Suppl. Table S4).

The identification of both, quantitative and qualitative resistance to downy mildew in our quinoa collection suggest that both, major and minor genes may be involved in resistance to this pathogen depending on the accession. In agreement, in addition to the region on chromosome 4 commented above, our GWAS analysis identified several other regions associated with downy mildew resistance. Therefore, additionally to the major gene present in accession PI614911, there are, probably, other genes conferring resistance to this important disease in quinoa germplasm. In a previous study [16], we demonstrated the presence and high relevance of hypersensitive response as a defense mechanism against P. variabilis in quinoa. Hypersensitive response, which is a pathogen-induced cell death process at the site of infection that limits pathogen growth, is a common mechanism of resistance against biothrophic pathogens, as downy mildew. HR is the result of the recognition of pathogen effectors by the plant, unleashing effector-triggered immunity (ETI), and is activated by R-genes. In agreement with that, many of the regions associated with resistance to P. variabilis, according to our GWAS analysis, harbour plant receptor genes or resistance genes (Table 2). Especially exciting is the presence of nine of these genes located exactly in the same position as the MTAs identified by GWAS. These genes include two genes annotated as “disease resistance RPP13-like protein”, the gene annotated as RGA2 mentioned above as candidate gene for the major gene conferring resistance to downy mildew in accession PI614911, one gene annotated as RGA1, other gene annotated as RGA3, another annotated as “serine/threonine-protein kinase”, a “F-box/LRR-repeat protein”, a “LRR receptor-like serine/threonine-protein kinase” and a “wall-associated-receptor kinase” encoding genes (Suppl. Table S4). Plant resistance gene analogs (RGAs) act as intracellular receptors that perceive the presence of pathogen effectors by direct binding of the pathogen effector proteins, or by monitoring the modification of host proteins after associating with the pathogen, to activate multiple defense signal transductions to restrict pathogen growth [38]. RGAs include nucleotide binding site leucine rich repeats, receptor like kinases, receptor like proteins, pentatricopeptide repeats and apoplastic peroxidases [39]. The presence of genes similar to RPP13 in two MTAs is especially attractive, since RPP13 is a resistance gene that confers resistance to downy mildew in Arabidopsis [40]. Downy mildew in Arabidopsis is caused by Peronospora parasitica, a member of the same genus as the pathogen causing downy mildew in quinoa (Peronospora variabilis).

Another remarkable outcome is the presence of genes encoding “Zinc finger BED domain-containing protein” in the genomic regions associated with resistance identified by GWAS. BED domains have been found to be integrated into plant resistance genes from different plant species [41]. This kind of genes are frequent in quinoa genome, however, the presence of this kind of genes in more than one candidate region associated with resistance suggest that this type of genes can also play a relevant role in resistance against downy mildew in quinoa.

To validate the candidate genes identified in our GWAS analysis, the same approach as mentioned for the major gene conferring downy mildew resistance in accession PI614911 could be followed. That is, sequencing these genes and performing gene expression studies in accessions showing contrasting profiles for the MTAs associated.

Disease severity values were correlated between scorings performed in the same location. A Pearson’s correlation coefficient as high as 0.85 was found between scorings carried out in Córdoba in 2019 and 2022, demonstrating the accuracy of the method used to evaluate the response to the disease. However, disease severity between different locations was not correlated. Furthermore, the analysis of variance performed also indicated that disease severity was influenced by the location (Suppl. Table S2). In agreement with that, according to GWAS, some genomic regions were found to be associated with resistance to downy mildew only in one location. This differential response to downy depending on the location suggest the presence of P. variabilis races/isolates with different virulence in Córdoba than in Guadajira. The presence of races in the pathosysthem P. variabilis-Chenopodium quinoa was already suggested by Ochoa et al. (1999). Reinforcing this hypothesis, in a previous article [16], we already commented different reactions to P. variabilis (complete resistance vs. high susceptibility) of some quinoa accessions in the screenings performed in Córdoba compared with screenings carried out in other countries. The presence of R-genes into the genomic regions associated with resistance (that are typically race-specific) and the existence of HR in P. variabilis-C. quinoa (a mechanism present frequently in race-specific resistance) reinforce this hypothesis. Therefore, the genomic regions associated with resistance to downy mildew only in one location may harbour genes providing race-specific resistance to the P. variabilis races present in one region but, that can be overcome by other races with different virulent pattern present in other regions. PCA and fastSTRUCTURE divided the quinoa collection in two and four populations, respectively. Interestingly, there was a correlation between the response of the quinoa accessions to P. variabilis in each location and the belonging population (Suppl. Table S5), supporting the hypothesis of a differential pattern of resistance genes depending on the genetic population. So that, in general, quinoa accessions belonging to population 2 according to PCA (corresponding to highland quinoas), were more susceptible to downy mildew than those belonging to population 1 (corresponding to lowland quinoas) in Córdoba (Suppl. Table S5). Thus, the difference in average disease severity between lowland and highland accessions was statistically significant for the years 2021 and 2022. The opposite trend was observed in the case of Guadajira, where lowland quinoas were, as average, more susceptible than highland ones. In 2019, the same trend was observed, although differences were not statistically significant. In agreement, regarding the four populations identified by fastSTRUCTRE, lowland accessions, belonging to population 4, were, on average, the most susceptible accessions in Guadajira. Highland accessions were divided in 3 subpopulations in the analysis performed by fastSTRUCTRE and, similarly, population 2, corresponding to highland accessions mainly from Bolivia, was the population that showed higher disease severity, on average, in Córdoba (Suppl. Table S5).

Conclusion

Despite the relevance of downy mildew disease in quinoa, little was known about the genetic control of resistance to this disease. We here identified a set of genomic regions associated with this trait and provide plausible candidate genes located within these regions. The enrichment of these regions in plant receptors and resistance genes point out to a high relevance of gene-by-gene interactions controlling resistance/susceptibility to P. variabilis in quinoa. What more, one of these regions identified by GWAS was confirmed by BSA, and found to harbour a single dominant gene conferring complete resistance to downy mildew. All these outcomes markedly increased our current knowledge about the genetics of resistance to downy mildew in quinoa, providing valuable information for breeding for resistance to this important disease.

Data availability

The raw sequencing data that support the findings of this study have been deposited in the NCBI Sequence Read Archive (SRA) under the BioProject PRJNA1040439.

References

  1. Zurita-Silva A, Fuentes F, Zamora P, Jacobsen SE, Schwember AR. Breeding quinoa (Chenopodium quinoa Willd.): potential and perspectives. Mol Breed. 2014;34:13–30.

    Article  Google Scholar 

  2. Christensen SA, Pratt DB, Pratt C, Nelson PT, Stevens MR, Jellen EN, Coleman CE, Fairbanks DJ, Bonifacio A, Maughan PJ. Assessment of genetic diversity in the USDA and CIP-FAO international nursery collections of quinoa (Chenopodium quinoa Willd.) Using microsatellite markers. Plant Genet Resour. 2007;5(2):82–95.

    Article  CAS  Google Scholar 

  3. Risi C, Galwey NW. The Chenopodium grains of the Andes: Inca crops for modern agriculture. Adv Appl Biol. 1984;10:145–216.

    Google Scholar 

  4. Anuradha KM, Zinta G, Chauhan R, Kumar A, Singh S, Singh S. Genetic resources and breeding approaches for improvement of amaranth (Amaranthus spp.) and quinoa (Chenopodium quinoa). Front Nutr. 2023;10:1129723.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Bazile D, Jacobsen SE, Verniau A. The global expansion of quinoa: Trends and limits. Front Plant Sci. 2016;7:622.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Angeli V, Silva PM, Massuela DC, Khan MV, Hamar A, Khajehei F, Grae-Hönninger S, Piatti C. Quinoa (Chenopodium quinoa Willd.): an overview of the potentials of the Golden Grain and socio-economic and environmental aspects of its cultivation and marketization. Foods. 2020;9:216.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Repo-Carrasco R, Espinoza C, Jacobsen SE. Nutritional value and use of the Andean crops quinoa (Chenopodium quinoa) and kañiwa (Chenopodium pallidicaule). Food Rev Int. 2003;19:179–89.

    Article  Google Scholar 

  8. Abugoch L, Castro E, Tapia C, Anon MC, Gajardo P, Villarroel A. Stability of quinoa flour proteins (Chenopodium quinoa Willd.) During storage. Int J Food Sci Technol. 2009;44:2013–20.

    Article  CAS  Google Scholar 

  9. Gonzalez JA, Konishi Y, Bruno M, Valoy M, Prado FE. Interrelationships among seed yield, total protein and amino acid composition of ten quinoa (Chenopodium quinoa) cultivars from two different agroecological regions. J Sci Food Agr. 2012;92:1222–9.

    Article  CAS  Google Scholar 

  10. Aloisi I, Parrotta L, Ruiz KB, Landi C, Bini L, Cai G, Biondi S, Del Duca S. New insight into quinoa seed quality under salinity: changes in proteomic and amino acid profiles, phenolic content, and antioxidant activity of protein extracts. Front Plant Sci. 2016;7:656.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Escuredo O, Martin MIG, Moncada GW, Fischer S, Hierro JMH. Amino acid profile of the quinoa (Chenopodium quinoa Willd.) Using near infrared spectroscopy and chemometric techniques. J Cereal Sci. 2014;60:67–74.

    Article  CAS  Google Scholar 

  12. Filho AM, Pirozi MR, Borges JT, Pinheiro Sant’Ana HM, Chaves JB, Coimbra JS. Quinoa: Nutritional, functional, and antinutritional aspects. Crit Rev Food Sci Nutr. 2017;57:1618–30.

    Article  PubMed  Google Scholar 

  13. Danielsen S, Jacobsen SE, Echegaray J, Ames T. Impact of downy mildew on the yield of quinoa. In: CIP, editors. CIP Program Report 1999–2000. Peru: Lima; 2000. pp. 397–401.

    Google Scholar 

  14. Ochoa J, Frinking HD, Jacobs T. Postulation of virulence groups and resistance factors in the quinoa/downy mildew pathosystem using material from Ecuador. Plant Pathol. 1999;48:425–30.

    Article  Google Scholar 

  15. Mhada M, Ezzahiriand B, Benlhabib O. Assessment of downy mildew resistance. (Peronospora farinosa) in a quinoa (Chenopodium quinoa Willd.) Germplasm. Inter J Agr Biol Eng. 2014;8:277–80.

    Google Scholar 

  16. Calderón-Gónzález A, Matías J, Cruz V, Molinero-Ruiz L, Fondevilla S. Identification and characterization of sources of resistance to Peronospora Variabilis in quinoa. Agronomy. 2023;13:284.

    Article  Google Scholar 

  17. Benlhabib O, Boujartani N, Maughan PJ, Jacobsen SE, Jellen EN. Elevated genetic diversity in an F2:6 population of quinoa (Chenopodium quinoa) developed through an inter-ecotype cross. Front Plant Sci. 2016;7:1222.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Colque-Little C, Abondano MC, Lund OS, Amby DB, Piepho HP, Andreasen C, Schmöckel S, Schmid K. Genetic variation for tolerance to the downy mildew pathogen Peronospora Variabilis in genetic resources of quinoa (Chenopodium quinoa). BMC Plant Biol. 2021;21:41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Patiranage DSR, Rey E, Emrani N, Wellman G, Schmid K, Schmöckel SM, Tester M, Jung C. Genome-wide association study in quinoa reveals selection pattern typical for crops with a short breeding history. Elife. 2022;8:11.

    Google Scholar 

  20. Danielsen S, Munk L. Evaluation of disease assessment methods in quinoa for their ability to predict yield loss caused by downy mildew. Crop Prot. 2004;23:219–28.

    Article  Google Scholar 

  21. Premachandra HKA, la Cruz FLD, Takeuchi Y, et al. Genomic DNA variation confirmed Seriola lalandi comprises three different populations in the Pacific, but with recent divergence. Sci Rep. 2017;7:9386.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Sansaloni C, Petroli C, Jaccoud D, Carling J, Detering F, Grattapaglia D, Kilian A. Diversity arrays technology (DArT) and next generation sequencing combined: genome-wide, high throughput, highly informative genotyping for molecular breeding of Eucalyptus. BMC Proc. 2011;5(Suppl 7):54.

    Article  Google Scholar 

  23. Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE. 2011;6(5):e19379.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Raj A, Stephens M, Pritchard JK, FastSTRUCTURE. Variational inference of population structure in large SNP data sets. Genetics. 2014;197(2):573–89.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Chang CC, Chow CC, Tellier LCAM, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Giga Sci. 2015; 4(1), s13742-015-0047–0048.

  26. Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.

    Article  CAS  PubMed  Google Scholar 

  27. Hill WG, Weir BS. Variances and covariances of squared linkage disequilibria in finite populations. Theor Popul Biol. 1988;33(1):54–78.

    Article  CAS  PubMed  Google Scholar 

  28. Remington DL, Thornsberry JM, Matsuoka Y, Wilson LM, Whitt SR, Doebley J, Kresovich S, Goodman MM, Buckler IV. Structure of linkage disequilibrium and phenotypic associations in the maize genome. PNAS. 2001;98(20):11479–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. R Core Team. R: A Language and Environment for Statistical Computing. 2023. https://www.R-project.org/.

  30. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multipletesting. J R Statist Soc B. 1995;57:289–300.

    Article  Google Scholar 

  31. Newell MA, Asoro FG, Scott MP, White PJ, Beavis WD, Jannink JL. Genome-wide association study for oat (Avena sativa L.) beta-glucan concentration using germplasm of worldwide origin. Theor Appl Genet. 2012;125:1687–96.

    Article  CAS  PubMed  Google Scholar 

  32. Peterson A, Jacobsen SV, Bonifacio A, Murphy K. A crossing method for quinoa. Sustainability. 2015;7:3230–43.

    Article  Google Scholar 

  33. Fondevilla S, Rubiales D, Moreno MT, Torres AM. Identification and validation of coupling and repulsion phase RAPD and SCAR markers linked to the gene Er3 conferring resistance to Erysiphe pisi DC in pea. Mol Breed. 2008;22:193–200.

    Article  CAS  Google Scholar 

  34. Li B, Zhao Y, Zhu Q, Zhang Z, Fan C, Sikandar A, Gao P, Luan F. Mapping of powdery mildew resistance genes in melon (Cucumis melo L.) by bulked segregant analysis. Sci Hort. 2017;220:160–7.

    Article  CAS  Google Scholar 

  35. Abe A, Kosugi S, Yoshida K, Natsume S, Takagi H, Kanzaki H, Matsumura H, Yoshida K, Mitsuoka C, Tamiru M. Genome sequencing reveals agronomically important loci in rice using MutMap. Nat Biotechnol. 2012;30:174–8.

    Article  CAS  PubMed  Google Scholar 

  36. Xu X, Luo W, Guo J, Chen H, Akram W, Xie D, Li G. Fine mapping and candidate gene analysis of the yellow petal gene ckpc in Chinese kale (Brassica oleracea L. Var. Alboglabra Bailey) by whole-genome resequencing. Mol Breed. 2019;39:96.

    Article  Google Scholar 

  37. Shen R, Messer PW. Predicting the genomic resolution of bulk segregant analysis. G3. 2022;12(3):jkac012.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Zhou Z, Bar I, Sambasivam PT, Ford R. Determination of the Key Resistance Gene analogs involved in Ascochyta Rabiei Recognition in Chickpea. Front Plant Sci. 2019;10:644.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Sekhwal MK, Li P, Lam I, Wang X, Cloutier S, You FM. Disease Resistance Gene analogs (RGAs) in plants. Int J Mol Sci. 2015;16(8):19248–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Bittner-Eddy PD, Crute IR, Holub EB, Beynon JL. RPP13 is a simple locus in Arabidopsis thaliana for alleles that specify downy mildew resistance to different avirulence determinants in Peronospora Parasitica. Plant J. 2000;21(2):177–88.

    Article  CAS  PubMed  Google Scholar 

  41. Zuluaga AP, Bidzinski P, Chanclud E, Ducasse A, Cayrol B, Gomez Selvaraj M, Ishitani M, Jauneau A, Deslandes L, Kroj T, Michel C, Szurek B, Koebnik R, Morel J-B. The Rice DNA-Binding protein ZBED controls stress regulators and maintains Disease Resistance after a mild Drought. Front Plant Sci. 2020;11:1265.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Susana Vilariño (Algosur S.L, Spain) for providing seeds of quinoa cultivar F16.

Funding

Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. This research was funded by MCIU/AEI/https://doi.org/10.13039/501100011033/FEDER, UE, projects PID2019-103978RB-I00 and PID2022-140238OB-I00, and by Conserjería de Transformación Económica, Industria, Conocimiento y Universidades, Junta de Andalucía, and FEDER, projects P18-TP-592 and QUAL21_023 IAS. In addition, the authors gratefully acknowledge the financial support received from the MCIN/AEI NutriCrop Spanish Network RED2022-134382-T.

Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Author information

Authors and Affiliations

  1. Institute for Sustainable Agriculture, Consejo Superior de Investigaciones Científicas (CSIC), Avd. Menéndez Pidal s/n, Córdoba, 14004, Spain

    Sara Fondevilla, Álvaro Calderón-González & Borja Rojas-Panadero

  2. Agrarian Research Institute “La Orden-Valdesequera” of Extremadura, Centro de Investigaciones Científicas y Tecnológicas de Extremadura (CICYTEX), Autovía A-5, km 372 - 06187, Guadajira, Badajoz, 06187, Spain

    Verónica Cruz & Javier Matías

Authors
  1. Sara Fondevilla
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Álvaro Calderón-González
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Borja Rojas-Panadero
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Verónica Cruz
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Javier Matías
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

S.F. and J.M conceptualized the study. S.F., B.R-P., A.C-G., J.M. and V.C. analysed the data. S.F., J.M. and V.C. carried out the field screenings, A.C-G. extracted the DNA. S.F. drafted the article. All authors reviewed and approved the final manuscript.

Corresponding author

Correspondence to Sara Fondevilla.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fondevilla, S., Calderón-González, Á., Rojas-Panadero, B. et al. Genome-wide association study, combined with bulk segregant analysis, identify plant receptors and defense related genes as candidate genes for downy mildew resistance in quinoa. BMC Plant Biol 24, 594 (2024). https://doi.org/10.1186/s12870-024-05302-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-024-05302-2

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us