RNA-Seq Analysis of Quinoa Varieties with Different β-Carotene Contents Reveals Genes Related to Carotenoid Biosynthesis

1 Introduction

Quinoa (Chenopodium quinoa Willd.), a nutrient-dense pseudocereal crop, has garnered widespread attention for its exceptional environmental adaptability and promising application potential in functional food development [1]. Most existing studies have focused on the nutritional composition and bioactive components of quinoa grains [2,3,4]. Beyond grains utilization, quinoa greens represent a high-quality leafy vegetable resource, characterized by abundant nutrients, a short growth cycle, and year-round cultivability [5,6]. However, investigations into the nutritional properties and phytochemical profiles of quinoa leaf tissues remain relatively limited.

β-carotene is a ubiquitous plant secondary metabolite and one of the most stable and prevalent natural pigments. As a typical lipophilic compound, it functions as the primary precursor of provitamin A [7,8]. Moreover, β-carotene participates in plant photosynthetic processes and antioxidant defense [9]. Accordingly, β-carotene content is widely recognized as a core indicator for evaluating the nutritional quality of leafy vegetables. Previous comparative studies have demonstrated that the quinoa seedlings accumulate higher β-carotene than some conventional vegetables, including yellow pumpkin, cabbage, cauliflower and spinach [10,11]. Despite these findings, current researches on quinoa carotenoids remains restricted to quantitative detection and phenotypic evaluation, whereas the molecular mechanisms governing carotenoid accumulation and the key regulatory genes involved remain poorly understood. Benefiting from the high genetic diversity of quinoa germplasm resources [12], carotenoid content varies dramatically across different accessions, which provides a solid material foundation for dissecting the regulatory mechanism of β-carotene biosynthesis in quinoa.

High-throughput RNA sequencing (RNA-seq) has emerged as one of the most efficient approaches for exploring the molecular regulatory networks underlying diverse biological processes. Furthermore, the release of a high-quality quinoa reference genome has greatly advanced transcriptome-based investigations into the regulatory mechanisms of quinoa growth and development [13]. To date, many transcriptome studies have been conducted in quinoa to elucidate mechanisms associated with stress response, flavonoid biosynthesis, nitrogen deficiency adaptation and low phosphorus response [14,15,16,17]. Nevertheless, transcriptomic insights into β-carotene biosynthesis in quinoa leaves remain scarce.

In this study, we performed large-scale screening of leaf β-carotene concentrations across a diverse collection of quinoa germplasms. Four accessions with significantly contrasting β-carotene levels were selected for transcriptome analysis. This study aimed to identify DEGs related to β-carotene biosynthesis, characterize the potential molecular basis underlying β-carotene accumulation in quinoa leaves, and screen key candidate genes involve in carotenoid metabolism. Our findings will provide molecular insights into the genes and mechanisms related to quinoa nutritional quality.

2 Plant Materials Culture

2.1 Plant Materials

All quinoa germplasm materials used in this study were preserved at the China (Shanghai) international quinoa innovation center. All accessions were grown in a greenhouse under standardized growth conditions, including a 16 h light/8 h dark photoperiod, constant temperature of 23 ± 2°C, 60% relative humidity, and a photosynthetic photon flux density of approximately 6.02214 × 1021 μmol·m−2·s−1. Leaf samples were collected at the seedling stage 30 days after sowing, when plants reached an approximate height of 30 cm for subsequent phytochemical measurement. All samples were instantly snap-frozen in liquid nitrogen and preserved at −80°C for downstream experiments. Three independent biological replicates were prepared for each germplasm accession.

2.2 Determination of β-Carotene Concertration

The β-carotene content in quinoa leaves was quantified via a spectrophotometric approach reported in a previous study [18]. In brief, 0.15 g of lyophilized leaf powder was homogenized with a mixed solution containing 4 mL acetone and 6 mL n-hexane. The homogenate was fully vortexed for 1 min, and then centrifuged at 10,000× g for 10 min to separate the supernatant. The absorbance values of the collected supernatant was detected at four specific wavelengths (453 nm, 505 nm, 645 nm, and 663 nm). The final β-carotene contents were calculated using the established formula: β-carotene (mg/g) = 0.216 × A663 − 1.220 × A645 − 0.304 × A505 + 0.452 × A453.

2.3 RNA Extraction, Sequencing and Bioinformatic Analysis

Total RNA was isolated from frozen quinoa leaf samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer’s instructions. RNA concentration, purity and integrity were subsequently evaluated using an Agilent 2100 Bioanalyzer system (Agilent Technologies, Palo Alto, CA, USA). Qualified RNA samples were used for cDNA libraries construction with the Illumina RNA-Seq Library Preparation Kit (Illumina, Inc., San Diego, CA, USA). High-throughput sequencing was performed on the Illumina HiSeq™ 4000 platform by Shanghai OEbiotech Co., Ltd. (Shanghai, China). Raw sequencing reads underwent quality filtering to obtain high-quality clean reads: adapter sequences were removed, and reads with a Phred quality scores < 30 or sequence length < 50 bp were excluded. The retained clean reads were aligned to the Chenopodium quinoa reference genome retrieved from the Phytozome database (https://phytozome-next.jgi.doe.gov/info/Cquinoa_v1_0).

Differential gene expression analysis was conducted using the R package DESeq2. Genes with a p-value < 0.05, and |log2(fold change)| ≥ 1 were defined as significantly DEGs. Venn diagrams were generated to visualize shared and unique DEGs across different comparison groups were plotted on the OEbiotech online bioinformatics platform. KEGG enrichment analysis was performed via the same online platform to annotate the functional categories and metabolic pathways of the identified DEGs.

2.4 Quantitative Real-Time Analysis

qRT-PCR was carried out to validate the reliability of the transcriptome data. In brief, 1 μg total RNA from each sample was reverse-transcribed into first strand cDNA synthesis using a commercial reverse transcription kit (Toyobo, Osaka, Japan) following the manufacturer’s protocols. Gene-specific primer pairs for target genes were designed using the online tool (https://www.primer3plus.com/). Sequences of target genes and corresponding primers are listed in Table S1.

qRT-PCR amplification was performed on the ABI 7500 Fast Real-Time PCR system (Applied Biosystems, CA, USA) with SYBR qPCR Mix (Toyobo). The Cq-tubulin gene (gene ID: AUR62018378) was selected as the internal reference for normalization of gene expression leavels. The relative expression abundance of target genes were quantified based on the 2−∆∆CT algorithm [19].

3 Results

3.1 Screening of Quinoa Germplasms with Contrasting Leaf β-Carotene Contents

In this study, we determined leaf β-carotene concentrations in 196 quinoa germplasm accessions cultivated in a controlled growth chamber (Table S2). Based on the initial data in β-carotene accumulation, ten accessions with high and low β-carotene levels were selected for further validation of phenotypic stability (Fig. 1). Finally, two accessions with the highest β-carotene accumulation (Cq39 and Cq138) and two with the lowest β-carotene contents (Cq35 and Cq92) were chosen as experimental materials for subsequent transcriptomic analysis.

Figure 1: β-carotene contents in leaves of 10 selected quinoa germplasm accessions. Different lowercase letters indicate statistically significant differences at p < 0.05 level.

3.2 Transcriptome Sequencing and Data Quality Assessment

To explore the potential molecular regulatory mechanism of β-carotene biosynthesis in quinoa leaves, total RNA was isolated from leaf tissues of the four screened germplasms, and 12 digital gene expression libraries were constructed for high-throughput transcriptome sequencing. As summarized in Table S3, each library generated an average of 24.46 million raw reads. After rigorous quality filtering, and the clean read count per sample ranged from 20.82 to 24.47 million. All clean reads were mapped to the quinoa reference genome, yielding a high mapping rate of 96.53% to 98.73% across individual libraries. Moreover, Phred value > 30 (Q30) of all sequencing libraries varied between 93.52% and 95.04%, indicating high sequencing accuracy and reliability.

Principal component analysis (PCA) was performed to assess the reproducibility and consistency of biological replicates. PCA results showed that three replicates of each germplasm clustered closely with low intra-group variability (Fig. S1), confirming the high quality of transcriptome datasets for subsequent bioinformatic analysis.

To further validate the accuracy of gene expression patterns derived from RNA-seq, seven candidate genes were randomly selected for qRT-PCR detection. The results showed that the expression trends of these selected genes were highly consistent between transcriptome sequencing and qRT-PCR results, with a Pearson correlation coefficient (R2) of 0.837 (Fig. S2). The high correlation confirmed that the gene expression profiles generated by RNA-seq were reliable for subsequent differential expression analysis.

3.3 Identification of Differentially Expressed Genes

To characterize the genome-wide transcriptional features of the four contrasting quinoa germplasms, global gene expression profiles were analyzed. The overall gene expression abundance varied among the four accessions (Fig. 2A). A hierarchical clustering heatmap was generated based on the expression levels of all transcripts, which revealed distinct global expression patterns among the four germplasms (Fig. 2B). However, similar expression patterns were observed between the two high-carotenoid quinoa germplasms and between the two low-carotenoid quinoa germplasms, respectively. Specifically, genes in Cluster I and Cluster II exhibited higher transcript abundance in Cq39 and Cq138, whereas genes in Cluster III and Cluster IV presented lower expression levels in these two high-carotene accessions.

Figure 2: Comparative analysis of DEGs in four quinoa germplasms. (A) Genome-wide gene expression distribution in four quinoa germplasms. The y-axis represents gene number. Blue boxes indicate genes with fragments per kilobase per million reads (FPKM) values of 1–10, and orange boxes indicate genes with FPKM ≥ 10. Numbers above bars represent the total gene count. (B) Hierarchical clustering heatmap of genome-wide gene expression across four germplasms. The bottom color bar represents the log10 of FPKM for each gene, ranging from blue (−1.5) to orange (1.5). A deeper color indicates more transcript accumulation. (C) Number of DEGs in pairwise comparisons between high-β-carotene and low-β-carotene accessions. The y-axis represents DEG count.

3.4 KEGG Pathway Enrichment Analysis of DEGs

To screen key transcripts associated with carotenoid metabolism, Venn diagrams were performed to identify overlapping DEGs among distinct comparison groups (Fig. 3). A total of 386 upregulated and 390 downregulated DEGs were commonly detected across all four germplasm comparison groups. Given the inherent transcriptional differences between two high-β-carotene germplasms, we further focused on consensus DEGs shared in pairwise comparisons between high- and low-β-carotene accessions. A total of 123 upregulated and 302 downregulated DEGs were commonly identified in both Cq39 vs. Cq35 and Cq138 vs. Cq35 comparisons. Additionally, 1278 upregulated and 946 downregulated DEGs were identified between the Cq39 vs. Cq92 and Cq138 vs. Cq92 comparison sets.

Figure 3: Venn diagram of differentially expressed genes in pairwise comparisons between high-β-carotene and low-β-carotene accessions. Red and blue numbers represent upregulated and downregulated DEGs, respectively.

KEGG pathway enrichment analysis was carried out to annotate the biological functions of overlapping intersecting DEGs, pathways with p < 0.05 were defined as significantly enriched functional categories.

For the 776 DEGs identified across all four comparison groups, six pathways including sesquiterpenoid and triterpenoid biosynthesis were significantly enriched among the 386 upregulated DEGs (Fig. 4A,B). none of these enriched pathways were directly associated with β-carotene biosynthesis. Seven pathways were significantly enriched for the 390 downregulated DEGs, among which terpenoid backbone biosynthesis was the one of the most representative categories.

Figure 4: KEGG pathway enrichment analysis of 386 upregulated and 390 downregulated common DEGs across all comparisons. (A) Enrichment results of 386 upregulated DEGs. (B) Enrichment results of 390 downregulated DEGs.

For the 425 common DEGs in the high-carotene vs. Cq35 comparison set, β-alanine metabolism and anthocyanin biosynthesis were significantly enriched in the 123 upregulated DEGs (Fig. 5A,B). In the 302 downregulated DEGs, three pathways including protein export, caffeine metabolism and carotenoid biosynthesis were identified as significantly enriched categories.

Figure 5: KEGG pathway enrichment analysis of 123 upregulated and 302 downregulated common DEGs in high-carotene vs. Cq35 comparisons. (A) Enrichment results of 123 upregulated DEGs. (B) Enrichment results of 302 downregulated DEGs.

For the 2224 DEGs specifically shared in the high-β-carotene vs. Cq92 comparison set, multiple pathways were significantly enriched (Fig. 6A,B), yet none of these enriched pathways were associated with β-carotene biosynthesis.

Figure 6: KEGG pathway enrichment analysis of 1278 upregulated and 946 downregulated common DEGs in high-carotene vs. Cq92 comparisons. (A) Enrichment results of 1278 upregulated DEGs. (B) Enrichment results of 946 downregulated DEGs.

3.5 Expression Profiles of Core Genes in the Carotenoid Metabolic Pathway

We further analyzed expression patterns of core genes involved in carotenoid biosynthesis, degradation, and abscisic acid (ABA) metabolism, and constructed a heatmap based on their log10-transformed FPKM values (Fig. 7). The results revealed clear divergence in transcriptional profiles of some functional genes between high- and low-β-carotene accessions.

Most structural genes in the carotenoid biosynthetic pathway, including geranylgeranyl pyrophosphate synthase (GGPS), phytoene synthase (PSY), phytoene desaturase (PDS), ζ-carotene desaturase (ZDS), lycopene β-cyclase (LCYB), β-carotene hydroxylase 1 (BCH1), zeaxanthin epoxidase (ZEP), and violaxanthin de-epoxidase (VDE), exhibited comparable transcript abundance across all four quinoa germplasms.

For genes responsible for carotenoid degradation and ABA metabolism, overall expression levels showed no global difference among germplasms. Nevertheless, one homolog of carotenoid cleavage dioxygenase 4 (CCD4, gene ID: AUR62007366) showed significantly lower transcript abundance in high-β-carotene lines. Similarly, two homologs of 9-cis-epoxycarotenoid dioxygenase (NCED, gene IDs: AUR62037951 and AUR62043087) were also downregulated in Cq39 and Cq138. In contrast, CCD1 expression remained stable across all tested samples.

Collectively, these results suggested that transcriptional suppression of CCD4 and partial NCED homologs, rather than widespread upregulation of biosynthetic genes, may contribute to the accumulation of β-carotene in high-content quinoa germplasms.

Figure 7: Expression heatmap of core genes related to β-carotene metabolism across four quinoa germplasms. Red arrows indicate genes with differential expression patterns between high- and low-β-carotene accessions.

3.6 Validation of Key Genes Associated with β-Carotene Content

Combining KEGG enrichment results and carotenoid pathway gene expression profiles, three key genes closely related to leaf β-carotene content were selected for qRT-PCR validation, including one CCD4 homolog (AUR62007366), and two NCED homologs (AUR62037951 and AUR62043087). The qRT-PCR results showed that all five genes exhibited significantly lower transcript levels in Cq39 and Cq138 compared with Cq35 and Cq92 (Fig. 8), which was consistent with the transcriptome data.

Figure 8: qRT-PCR analysis of key genes related to β-carotene genes in four quinoa germplasms. * indicates statistically significant differences (p < 0.05).

4 Discussion

β-carotene is a vital plant secondary metabolite and the primary dietary precursor of provitamin A, whose content directly governs the nutritional quality of leafy vegetables including quinoa greens. Previous research on quinoa carotenoids has largely been focused on phenotypic quantification and component profiling [4], while the transcriptional regulatory mechanisms driving natural variation in β-carotene accumulation across quinoa germplasms remain largely unexplored. In this study, we performed large-scale screening of leaf β-carotene content across 196 quinoa germplasm resources, and selected two high-β-carotene (Cq39 and Cq138) and two low-β-carotene (Cq35 and Cq92) representative lines for comparative transcriptomic analysis.

Genome-wide transcriptome profiling revealed substantial variation in gene expression patterns across the four quinoa germplasms. Specifically, comparisons between the high-carotene accession Cq39 and the two low-carotene lines identified more upregulated than downregulated DEGs, whereas the opposite trend was observed in comparisons involving the other high-carotene accession Cq138. Quinoa has evolved extensive genetic diversity following centuries of fragmented and localized cultivation in its native Andean region [20], which likely explains the divergent transcriptomic landscapes observed among the four germplasms. This finding also implies that different quinoa accessions may achieve high β-carotene accumulation through distinct transcriptional regulatory networks. Nevertheless, hierarchical clustering of all expressed genes showed that the two high-β-carotene accessions shared consistent expression trends for a subset of genes relative to low-carotene lines, suggesting that these concordantly regulated genes are potentially involved in the regulatory network governing β-carotene accumulation in quinoa leaves.

KEGG pathway enrichment analysis of the overlapping DEGs further revealed germplasm-specific functional enrichment patterns. Among DEGs shared by all four comparison sets, upregulated genes were primarily enriched in sesquiterpenoid and triterpenoid biosynthesis, with no significant enrichment of carotenoid metabolism-related pathways. Conversely, downregulated common DEGs were significantly enriched in the terpenoid backbone biosynthesis pathway. Terpenoids represent one of the largest and most chemically diverse families of natural products, with tens of thousands of structurally distinct derivatives identified across all kingdoms of life [21]. As typical tetraterpenoid pigments, carotenoids are synthesized from precursors derived from the terpenoid backbone pathway [22]. Combined with our observation that core carotenoid biosynthetic genes maintained stably expressed (Fig. 7), the global downregulation of terpenoid backbone genes is not a limiting factor for carotenoid biosynthesis.

Expression profiling of core genes in the carotenoid metabolic pathway further clarified the regulatory mechanism underlying β-carotene enrichment in quinoa leaves. All detected structural genes involved in carotenoid biosynthesis, including GGPS, PSY, PDS, ZDS, LCYB, BCH1, ZEP, and VDE [23], displayed comparable transcript levels between high- and low-β-carotene germplasms. This expression pattern indicates that the elevated β-carotene accumulation in Cq39 and Cq138 leaves is not driven by the transcriptional upregulation of core carotenoid biosynthetic genes. In this study, one CCD4 homolog and two NCED homologs exhibited significantly lower transcript abundance in high-β-carotene germplasms. CCDs are key enzymes that catalyze the irreversible oxidative cleavage of mature carotenoid backbones, leading to pigment degradation [24]. CCD4 has been reported negative regulate carotenoids in lots of plants, such as Rosa damascene, Malus domestica, Chrysanthemum morifolium, rice and soybean [25,26,27,28]. NCEDs are rate-limiting enzymes in abscisic acid (ABA) biosynthesis, which cleave 9-cis-violaxanthin and 9-cis-neoxanthin to produce xanthoxin, the direct precursor of ABA [29,30]. Suppression of NCED1 increases both lycopene and β-carotene contents in tomato fruit [31]. These suggest that the enhanced β-carotene accumulation in Cq39 and Cq138 may be mediated the downregulated of CCD4 and NCEDs expression.

In conclusion, this study identified two quinoa germplasms with high leaf β-carotene content through large-scale germplasm screening, and systematically dissected the transcriptional regulatory mechanism underlying natural variation in β-carotene accumulation in quinoa leaves. Our findings revealed that three genes, including one CCD4 homolog and two NCED homologs, may be associated with β-carotene enrichment in high-content quinoa accessions. Further functional characterization of these genes is required to clarify their exact biological roles and precise regulatory mechanisms in quinoa carotenoid metabolism in future studies.