Genome-wide identification, phylogeny and expression analysis of the R2R3-MYB gene family in quinoa (Chenopodium quinoa) under abiotic stress

Introduction

Quinoa (Chenopodium quinoa; 2n = 4x = 36) is a nutritious crop well adapted under varying agro-ecosystems around the world. Quinoa was first domesticated approximately 7000 years ago during pre-Columbian era, and the seed was known as ‘mother grain’ during Inca Empire of 13th century (Jarvis et al. 2017). Quinoa possesses high genetic diversity and ability to be adapted under harsh conditions; i.e. marginal soils and resilient to abiotic stresses, e.g. low and high temperature, drought, salinity and large day-night duration and temperature variations (Jacobsen et al. 2005, 2007; Ruiz et al. 2014). Quinoa produces highly nutritious and gluten-free grains with an excellent balance between carbohydrates, essential amino acids, oils, minerals, vitamins and dietary fibre (Graf et al. 2015; Imamura et al. 2018). Quinoa improvement through breeding approaches has been challenging owing to multiple reasons. Therefore, genetic based interventions hold promising effect for crop improvement. Therefore, identification of candidate genes could be reliable approach can be useful technique utilised for crop improvement. The development of quinoa genome is considered useful resource for understanding functional genomics of integralagronomic attributes and the identification of candidate genes being utilised for quinoa crop improvement programmes (Jarvis et al. 2017).

Approximately all biological processes taking place in eukaryotic cells or organisms are controlled through gene expression (Pucker et al. 2020). Therefore, gene transcription regulation level is an important starting point for genetic understanding. The transcription factors (TFs) (sequence-specific DNA-bindings factor) are the proteins regulating rate of transcription of genetic information from DNA to messenger RNA through linkage with a specific DNA sequence (Sabir et al. 2022). TFs can be divided into various multigene families as per their DNA-binding domains. One of the largest TFs families present in plant kingdom has been identified as MYB gene family, present among all eukaryotes. The COLORED1 (C1) was first plant MYB gene identified in Zea mays, which was responsible for regulating anthocyanin biosynthesis both in aleurone and scutellum tissues of maize kernel (Arce-Rodríguez et al. 2021). The highly conservative MYB domain of MYB gene family members consisted of one to four neighbouring incomplete tandem repeats. Furthermore, each repeated unit is comprised of approximately 50–55 amino acid residues, with three spaced tryptophan residues, forming a helix-turn-helix (HTH) fold. The MYB TFs possess several regulatory roles, the amino acid sequence outside the conserved MYB domain are highly divergent (Zhang et al. 2020). These TFs are further divided into four classes based on the number of MYB domain repeats: 1R-MYB, R2R3-MYB, 3R-MYB and 4R-MYB (Zhao et al. 2020). The TFs have diverse functions and may play an indispensable role during plant growth and stress resistance; i.e. such as cell morphogenesis, secondary metabolism, organ morphogenesis, chloroplast development, leaf senescence, response to phosphate starvation and tolerance to cold, drought and salt stress (Park et al. 2018; Sun et al. 2018). The R2R3-MYB, TFs are considered largest subfamily among MYB family in plants whereas, the number of 4R-MYB and 3R-MYB genes are relatively smaller.

During the past decade, substantial huge amount of genome sequencing data, The R2R3-MYB gene family has been comprehensively studied for gene identification in several plant genomes; e.g. 126 genes in Arabidopsis thaliana (Stracke et al. 2001), 192 in Populus ciliata (Wilkins et al. 2009), 244 in Glycine max (Du et al. 2012) and 406 in Gossypium hirsutum (Wang et al. 2019). In addition, literature has displayed R2R3-MYB TFs role for controlling physiological and biochemical mechanisms; i.e. responses to ambient stimuli (He et al. 2020), specialised metabolism (Yan et al. 2020a), plant development (Sun et al. 2019) and cell differentiation (Yan et al. 2020b).

Previous literature has displayed involvement of R2R3-MYB gene family may play central role during plant growth and development and response of plant under salinity stress condition. However, the role of R2R3-MYB gene family toward salt stress for quinoa plant species is limited. In this study, we identified R2R3-MYB TFs their comprehensive phylogenetic analysis, chromosomal distribution, gene duplication events and synteny analysis. Based on available transcriptomic data, the tissue differential expression patterns and response of gene family under salinity stress will be explored for the characterisation of CqMYB genes and providing novel insights and valuable information for the further identification of this R2R3-MYBs gene family in quinoa.

Materials and methods

Results

Phylogenetic analysis of the MYB family

To explore the relationships between the R2R3-MYB proteins, we constructed the NJ-phylogenetic tree by using the amino acid sequences from both Arabidopsis and quinoa. As per classification in Arabidopsis, quinoa R2R3-MYB proteins were divided into 27 groups. Among these, group 20 was identified as largest subgroup and contained nine CqR2R3-MYB proteins (Fig. 1).

Fig. 1.

Phylogenetic tree of MYB protein quinoa, and Arabidopsis thaliana with each colour representing a subfamily (S1-C27) of MYB genes.

Chromosomal localisation of the CqR2R3MYB gene family

The chromosomal distribution of R2R3-MYBs gene family were analysed via employing quinoa genome annotation and TBtools was employed for their visualisation (Jarvis et al. 2017; Chen et al. 2020). The chromosomal location analysis revealed that 103 CqR2R3-MYB genes were unevenly distributed across all 18 quinoa chromosomes, whereas five CqR2R3-MYB genes belonged to unmapped scaffolds. Chromosome 1 had 14 of CqR2R3-MYB genes, followed by 12 on chromosome 6. In contrast, only a single CqR2R3-MYB genes was identified on chromosome 6 (Fig. 2).

Fig. 2.

Chromosome distribution of CqMYB genes. The chromosome number is displayed next to each chromosome.

Gene structure and protein motif analysis of the CqR2R3-MYBs

The conserved motifs of CqR2R3-MYBs were examined through employing MEME program. We identified 20 conserved motifs; however, the number of motifs (3–8) were different in identified proteins. Generally, the location of introns and exons within the genome may provide important information on evolutionary relationships of quinoa. Results also displayed that number of exons in CqR2R3-MYBs ranged from two to five; among these, mostly contained five exons (Fig. 3). This is the first time that the distribution of exons in CqR2R3-MYBs has been identified, providing valuable information on its evolutionary processes and expansion in quinoa and other species.

Fig. 3.

Phylogenetic relationships and gene structures (introns/exons) MYB genes of the quinoa. The grey line represents the intron(s), and yellow and green colours represent the exon(s) and UTR(s).

Gene duplication analysis of CqR2R3MYBs

Gene duplication, resulting from polyploidisation or during tandem and segmental duplication is associated with replication, is an integral component responsible for expansion of any gene family. Of the 103 genes, 89 were CqR2R3-MYBs, which accounts for 86.4% of the total genes identified that participated in gene replication events. Furthermore, 10 genes took part in tandem duplication events and 89 genes participated in WGD or segmental events, forming four tandem duplication gene pairs (Fig. 4).

Fig. 4.

Duplicated CqMYB genes pairs identified in quinoa. Quinoa chromosomes are shown in different colours. Duplicated gene pairs are linked by lines with red lines.

Expression of CqR2R3MYBs genes in different tissues

To understand the tissue specificity and genes controlling stress response quinoa, the expression pattern of 103 CqR2R3-MYBs genes were studied. The RNA-seq data of various tissues of quinoa; i.e. stems, leaves and roots were downloaded from the SRA database (https://www.ncbi.nlm.nih.gov/sra). The FPKM values were calculated by Hisat (ver. 2.0.4, http://ccb.jhu.edu/software/hisat2/index.shtml) and Tophat software (ver. 2.1.1, http://ccb.jhu.edu/software/tophat/index.shtml). Approximately half of the CqR2R3-MYBs genes shown clear and similar tissue-specific expression (Fig. 5). CqMYB98 shown lowest expression in leaves, whereas the strongest expressions of genes were detected in the roots. In addition, the strongest expressions of CqMYB79 were found in stem. However, they showed low expression in rest of the tissues. Among the CqR2R3-MYBs genes, the most obvious differences in tissue expression were observed for genes CqMYB98, CqMYB103, CqMYB56, CqMYB95, CqMYB79 and CqMYB47.

Fig. 5.

Relative expression level of 103 CqMYB genes in specific tissues. l, leaf; r, root; s, stem.

Expression of CqMYB genes under salt treatment

To understand the stress response of CqR2R3-MYBs, both salttolerant and salt-sensitive quinoa tissue were selected to study the expression levels of CqR2R3-MYBs under salt stress condition. The response of CqR2R3-MYBs in root were observed significantly higher compared to rest of the tissues. Approximately two thirds of genes displayed significant expression under salt stress condition, suggesting their role under salt stress condition. In roots, the expression level of CqMYB79, CqMYB95, CqMYB49, CqMYB58, CqMYB75, CqMYB92, CqMYB90, CqMYB17, CqMYB3, CqMYB2, CqMYB16, CqMYB41, CqMYB67, CqMYB1 and CqMYB15 were significantly higher for both salt-tolerant and salt-sensitive quinoa genotypes under control condition. In contrast, the expression level of CqMY-B98 was significantly lower compared to control condition for both investigated genotypes. However, in shoot, the expression levels of CqMYB78, CqMYB96, CqMYB56 and CqMYB103 were slightly higher than control condition for both salt-tolerant genotype Q68 and salt-sensitive genotype Q30 (Fig. 6).

Fig. 6.

Relative expression level of 103 CqMYB quinoa genes under salt treatment for salt-tolerant genotype Q68, and salt-sensitive genotype Q30.

Discussion

In this study, R2R3-MYB TPs of Arabidopsis were studied as reference genes. A total of 103 CqR2R3-MYB genes were identified in quinoa. The size of the R2R3-MYB gene family in quinoa were much higher compared to Prunus avium (69 members) (Zhao et al. 2020; Sabir et al. 2022), but lower than that Phyllostachys edulis (114 members), and Arabidopsis (125 members) (Stracke et al. 2001; Hou et al. 2018). The results published by various studies and considering phylogenetic analysis we divided the quinoa R2R3-MYB family into 27 groups. The chromosome distribution results displayed R2R3-MYB genes were unevenly distributed on 17 chromosomes. The abundance of TFs were found largely dependent on sequence duplications during genome evolution process (Yue et al. 2018). The relatively high number of CqR2R3-MYB genes in quinoa indicated that gene duplication events may have been occurred during genome evolution. It can also be hypothesised that presence of most R2R3-MYB genes in quinoa genome contain specific requirements and are involved in complicated mechanism of transcriptional regulation. These findings may provide with an insights on evolutionary relationship of R2R3-MYB family genes.

To understand the CqMYB gene family evolutionary relationships, the R2R3-MYB proteins were divided into S1-C27 groups based phylogenetic tree constructed with Arabidopsis. The CqR2R3-MYB genes grouped with Arabidopsis orthologues in several branches, suggesting R2R3-MYB function was extensively conserved across the species. The phylogenetic connection of MYB proteins between quinoa and Arabidopsis also revealed clades with varied numbers of AtMYB and CqMYB genes, which demonstrate that these two species share the similar evolutionary history. In phylogenetic tree, several CqMYB genes from quinoa were not grouped with the AtMYBs from Arabidopsis, displaying these genes were not retained between quinoa and Arabidopsis and may hold unique role. Furthermore, these results indicated the existence of divergence among primitive ancestor, species-specific MYBs genes were included in quinoa or eliminated from Arabidopsis lineages.

The genes fall in identical group contain similar positions of exon-intron, its structure; however, few dissimilarities were also observed which may be outcome of gain, loss, or skimming of introns during the evolution of MYB gene family. Similar observations were recorded for quinoa, suggesting similar evolutionary process, suggesting elimination of few genes during evolutionary process (Du et al. 2012). We showed that several genes contain only single intron, indicating the loss of intron during the evolution process. In other studies, it has been concluded that plants during evolution, gain/loss introns which is outcome of selection pressure (Wang et al. 2016; Zhao et al. 2020). The non-existence of intronic region in genes may helpful in acceleration of evolution of genes through process of gene copying (Lecharny et al. 2003). Moreover, the function of genes without an intron and those with single intron were found identical in their evolution and functional properties. We demonstrated that genes conserved through the evolutionary process contain highly functional similarities, and the majority of CqMYB genes have similar motif compositions in the identical subfamily. A sum of 20 diverse motifs were identified among MYB genes of sweet cherry (Prunus avium) with each gene containing at least two identical motifs. The arrangement and amount of the 12 distinct motif types in quinoa expressed that MYB members were functionally divergent. The presence of similar motifs and intron/exon structure among identified genes suggest their functional similarities.

A number of studies have shown that R2R3-MYB TFs play crucial roles for regulating plant growth, development and response under abiotic stresses. We investigated differential expression analysis among various quinoa tissue to salt stress, and found that five genes that were differentially expressed in any two tissues; e.g. leaf, stem and root and further divided into three groups based on difference in expression patterns. By comparing RNA-Seq data collected before and after salt stress, we identified 15 genes having significantly higher expression compared to control condition for both salt-tolerant and salt-sensitive quinoa genotypes. Furthermore, our results showed that the majority of genes were only expressed in roots, followed by four genes expressed in shoots.

Conclusions

We found 103 R2R3-MYB gene family members classified into 27 subgroups of quinoa genome. The transcriptome data showed that expression abundance of MYB family members in quinoa was significant, with CqMYB98, CqMYB103, CqMYB56, CqMYB95 and CqMYB79 genes having the high expression level and significant tissue specificity. The expression of genes under salt stress suggested members of quinoa MYB gene family are responsible for salinity stress response in quinoa. However, CqMYB90 maintained highly significant response under salt stress In conclusion, this study provided insight on functional roles of CqR2R3-MYBs, for plant growth and abiotic stress. These results will be useful for our enhancing understanding on molecular mechanisms of stress tolerance in crops, especially quinoa.

Supplementary material

Supplementary material is available online.