The quinoa variety NL-6 and Longli-4 were obtained from Prof. Heng Zhang and Prof. Pengshan Zhao, respectively. Quinoa seeds were directly sown in soil and kept at a growth chamber at 28°C under a light cycle of 12 h light/12 h dark. The phenotype of plants was observed and photographed after growth for indicated days. Arabidopsis Columbia (Col-0) ecotype was used for the transformation of 35S::CqPORB-GFP. Arabidopsis seeds were sterilized and sown on half Murashige and Skoog (MS) solid medium and kept at 4°C for 3 days before being transferred to an incubator at 22°C with a long-day cycle (16 h light/8h dark).

For EMS treatment, a set of healthy seeds of quinoa NL-6 variety were presoaked in distilled water at 4°C for 8 h, and then the seeds were immersed in 1.0% EMS solution (v/v) (Sigma-Aldrich, M0880) and incubated at 4°C for 12 h. After treatment, a volume of 10% sodium thiosulphate (w/v) solution was added to the EMS solution, and the quinoa seeds were gently stirred and stand for 10 min before they were washed five times with distilled water to neutralize and remove EMS. The mutated seeds were sown in field for the collection of progeny seeds.

For Chl and carotenoid measurement, the leaves of four quinoa seedlings grown in the conditions of 28°C/26°C day/night temperatures and a 12 h/12 h light (216 µmol m⁻²·s⁻¹)/dark cycle were sampled. Total pigments were extracted using 95% ethanol. The contents of Chl a, Chl b and carotenoid were measured using a NanoDrop™ 2000C UV-vis spectrophotometer (Thermo Scientific) as described previously (Lichtenthaler, 1987).

Total Pchlide and Chlide were calculated following the method described previously (Kwon et al., 2017). 20 seedlings were grown in darkness for 5 days or further exposed to white light for 10 s, 30 s, and 5 min. The seedlings were homogenized in 500 μL of ice-cold 80% acetone and incubated for 8 h in darkness (4°C). After spinning for 5 min, 150 μL of the supernatant was transferred to a 96-well plate (black, LABSELECT, 33113) that is made of polypropylene. Fluorescence was detected at an excitation wavelength of 440 nm, and emissions were collected from top between 600 nm and 720 nm at room temperature using a Thermo Scientific Varioskan Flash Plate Reader.

The quinoa seeds were sown directly in soil and grown in a light incubator at 22°C (light)/20°C (dark) in long-day conditions (16 h light/8 h dark). The first or second pair of the true leaves of quinoa seedlings after growth for 15-20 days were collected using adhesive tape and dissociated in a dissociation solution (1.5% Cellulase R10, 0.4% Macerozyme R10, 0.5 M Mannitol, 20 mM KCl, 20 mM MES (pH 5.7), 10 mM CaCl₂, and 0.1% BSA). The dissociation was carried out on a shaker with a speed of 50 rpm for 2.5 h under darkness at 22°C. The protoplasts were filtered through a Filter Mesh 100 and the filtered cells were centrifuged at 100 g for 3 min. The supernatant was removed and the cells were washed with 5-10 mL of W5 solution (154 mM NaCl, 125 mM CaCl₂, 5 mM KCl, 2 mM MES (pH 5.7), and 0.1% (w/v) glucose) for three time with gentle shake. After washing, appropriate amount of MMg solution (0.4 M Mannitol, 15 mM MgCl₂, and 4 mM MES (pH 5.7)) was added to adjust the cell concentration to about 2 × 10⁵ cells/mL. Subsequently, 100 μL protoplasts, 5 μL CqPORB-GFP plasmid (800-1000 ng/μL), and 105 μL 40% (w/v) PEG (0.2 M mannitol and 100 mM Ca(NO₃)₂) were slowly mixed, and the mixture was incubated at room temperature for 15 min. The transformation was terminated by adding 500 μL W5 solution. After centrifugation at 100 g for 3 min, the supernatant was discarded and protoplasts were washed twice with W5 solution. The protoplasts were transferred to cell culture plate and incubated at 22°C under low light conditions (5 µmol·m⁻²·s⁻¹) for 36 h. Then protoplasts were collected and total Chl was extracted under light conditions, and Chl was measured according to the protocol described above.

The first true leaves of NL-6 and nl6-35 harvested from 10-day-old seedlings grown in low-light conditions were used for the detection of granum in chloroplast. The samples were immersed in TEM fixative (G1102, Servicebio, China) containing 2.5% glutaraldehyde overnight at 4°C, and rinsed with 0.1 M sodium phosphate buffer (pH 7.4), and post-fixed with a 1% (w/v) aqueous solution of OsO₄ at room temperature for 7 h. The materials were then dehydrated in a gradient ethanol series, and embedded in EMBed 812 resin (90529-77-4, SPI, USA). Tissues were processed into ultrathin sections, and stained with a 2% (w/v) uranium acetate and subsequently with 2.6% (w/v) lead citrate solution. The chloroplast was observed using a TEM device (HT7800, HITACHI, Japan).

For BSA-seq analysis, nl6-35 was backcrossed with its parent NL-6. In F₂ population, 80 plants with faded green leaf were collected and pooled together for total DNA isolation using the DNeasy Plant Maxi kit (Qiagen). The library construction, sequencing, and bioinformatic analysis were all conducted by Personalbio (Shanghai, China). 150 bp pair-end sequencing libraries were generated and sequenced with an Illumina NovaSeq 6000 platform (Illumina lnc., San Diego, CA, USA). The clean reads that were filtered from raw reads were aligned to the reference genome of NL-6 using BWA (Li and Durbin, 2009). Picard 1.107 tools were used to convert alignment files to bam files and remove duplicated reads. GATK was applied to conduct local realignment around Indels and single-nucleotide polymorphism (SNP) calling (Zhu et al., 2014). Annovar program (Wang et al., 2010) was used to annotate all the SNPs based on the GFF3 files of the reference genome of NL-6. The SNP-index (Takagi et al., 2013) was calculated for all SNPs to identify candidate regions associated with the mutant trait, and the SNPs with a read depth ≥ 8 and SNP-index ≥ 0.3 were retained, and the average SNP-index across a 1-Mb genomic interval was measured individually using a 100-kb sliding-window approach and then plotted against all the chromosomes of quinoa genome. Finally, the SNPs with SNP-index = 1 were chosen for the manual analysis of candidate genes. To verify the SNP in CqPORB gene in nl6-35 mutant, genomic sequence of CqPORB gene was amplified and sequenced using traditional Sanger sequencing. The gene-specific primers used for sequencing are listed in Supplementary Table 1 .

The identified SNP in CqPORB gene was used to design Kompetitive Allele Specific PCR (KASP) marker. The developed marker was used for the linkage analysis of the F₂ population plants derived from a hybridization between nl6-35 and its parent NL-6. Totally 93 plants with mutant phenotype and 95 plants with wild type phenotype were selected for KASP analysis. Besides, the mutant-type plants from the F₂ progenies of Longli-4 × nl6-35 crossing population were also selected for linkage analysis using the gene-specific KASP marker in CqPORB gene.

Homologs of CqPORB were obtained by a Blastp search of GenomeNet (www.genome.jp), Plant GDB (www.plantgdb.org), and The Rice Annotation Project Database. Multiple protein sequences were aligned by ClustalW and conserved regions were visualized by CLC sequence viewer 8. The phylogenetic tree was constructed by iqtree with the best-fit model and 1,000 replicates of ultrafast bootstrap (Minh et al., 2020). The resulting tree files were submitted to iTOL (Letunic and Bork, 2021) for modification.

To analyze the expression pattern of CqPORs in response to light, surface sterilized seeds were sowed on ½ MS solid medium and incubated at 22°C for 3 days, and then etiolated seedlings were collected after exposure to white light for 0, 3, and 6 h. For the analysis of CqPOR gene expression in different developmental stages, the first true leaves of 10-day-old seedlings, leaves of early flowering stage (30 d), and late flowering stage (50 d) were sampled. Three independent replicates were performed. Total RNAs of all these samples were extracted by using EastepTM Super Total RNA Extraction Kit (Promega, United States). cDNAs were synthesized using HiScript III RT SuperMix for qPCR (+gDNA wiper) according to the manufacturer’s instructions (Vazyme, Nanjing, China). qRT-PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and Bio-Rad CFX connect real time detection system (BIO-RAD, United States). Primers used for qRT-PCR are listed in Supplementary Table 1 . The transcript levels of all three CqPOR genes in four different tissues were obtained from previous RNA-seq data (Jiang et al., 2022). Gene expression patterns of CqPOR genes in different tissues were visualized in heat map using TBtools software package (Minh et al., 2020).

CqPORB was amplified using the cDNA template generated from the leaves of NL-6 using gene-specific primers ( Supplementary Table 1 ). CqPORB was constructed into plant expression binary vector pCAMBIA1305, which carries a GFP reporter gene driven by the cauliflower mosaic virus (CaMV) 35S promoter. The recombinant plasmid CqPORB-GFP was introduced into the competent cells of Agrobacterium tumefaciens GV3101 strain, which was further used for transformation into Nicotiana benthamiana and Arabidopsis thaliana (Columbia). For transient expression in N. benthamiana, GV3101 strain harboring CqPORB-GFP was infiltrated into leaves according to the protocol described previously (Liu et al., 2010). Transformation of Arabidopsis plants was performed by using floral dip method (Clough and Bent, 1998). To detect the subcellular localization of CqPORB, GFP fluorescence was detected by Leica confocal laser scanning microscope SP8 with 488 nm excitation light and 500-550 nm emission light.

To investigate the regulatory mechanism of Chl biosynthesis in quinoa, we performed EMS-mediated mutagenesis for quinoa NL-6 variety and screened for mutants that displayed a reduced Chl content. NL-6 is a coastal germplasm that originates from lowland Chilean and exhibits relatively stable agronomic traits (Benlhabib et al., 2016). Here we report one of the mutants, designed nl6-35, that showed faded green leaf when grown in both greenhouse and field ( Figure 1A , Supplementary Figures 1A, B ). Specifically, we noted that at cotyledon stage, there was no significant difference between the wild type NL-6 and nl6-35 ( Supplementary Figure 1C ), but the leaves of nl6-35 became faded green after growth for approximately 14 days under white light ( Figure 1A ). Accompanying the faded green leaf, the phenotype of curled leaf was also observed in the nl6-35 mutant, but not in the wild type NL-6 ( Figure 1B ). Genetic analysis indicated that the phenotypes of faded green leaf and curled leaf were tightly associated, suggesting that they are likely caused by a same mutation in the nl6-35 mutant. Quantification of Chl content revealed that Chl a and Chl b, and total Chl contents were significantly reduced in the nl6-35 mutant ( Figure 1C ), but the ratio of Chl a/Chl b was slightly increased in the nl6-35 mutant ( Figure 1D ). In addition, the carotenoid content was also reduced by approximately 15.6% in the nl6-35 ( Figure 1E ). These results suggested that the faded green leaf in nl6-35 was likely caused by the reduced Chl and carotenoid contents.

To identify the mutated gene in nl6-35 that is responsible for the phenotype of faded green leaf, nl6-35 was backcrossed with its parent NL-6. The F₁ generation showed a wild type phenotype and F₂ population exhibited a segregation of 3:1 ( Table 1 ), indicating that nl6-35 is a recessive mutation controlled by a single gene. To perform BSA-seq analysis, 80 plants of F₂ population with faded green leaf were collected and mixed for whole genome sequencing. After filtering the raw reads, a total of 324,781,844 high quality reads were obtained, and 99.94% of the reads were mapped to NL-6 reference genome (1.39 GB), with an average of 34× reads coverage. Totally 11,960 SNPs were identified after removing the ones with low coverage (<8 reads). Because EMS mutagenesis most frequently causes C to T or G to A change in the genome, so we searched for these two types of nucleotide acid changes and finally obtained 3,013 SNPs with C to T or G to A change ( Figure 2A , Supplementary Figure 2 , Supplementary Table 2 ). Among these SNPs, 35 were identified with a mutation rate of 100%, and 15 of them were located on chromosome (chr)17 ( Supplementary Table 3 ), suggesting that the objective gene is most likely located on chr17. Analysis of these 15 SNPs showed that 14 SNPs were located either in intronic or intergenic regions ( Supplementary Table 3 ), while one SNP (G2258A) was identified in the fourth exon region of CqNL-6_047471, resulting in an amino acid change of Leu279 to Phe ( Figure 2B ). CqNL-6_047471 gene was annotated as a light-dependent protochlorophyllide reductase, which is homologous to POR genes reported in other plant species. It is worth noting that the faded green leaf of the nl6-35 mutant is similar to the phenotype reported in OsporB mutant in rice (Kwon et al., 2017), suggesting that CqNL-6_047471 gene mutation is responsible for the phenotype of nl6-35 mutant in quinoa. Because CqNL-6_047471 gene mutation resulted in a similar phenotype as the OsporB mutant in rice, we named CqNL-6_047471 as CqPORB.

Currently genetic transformation in quinoa is still technically limited, so we explored other strategies to verify the association of CqPORB gene mutation with the faded green leaf of nl6-35. In the F₂ hybridization population of nl6-35 with NL-6, we selected 93 plants that showed faded green leaf phenotype and 95 plants that showed wild type phenotype for the analysis of G2258A mutation in CqPORB gene. To perform genotyping in a high-throughput way, KASP (Kompetitive Allele Specific PCR) molecular marker that targeted this mutation site was designed. KASP analysis showed that in all 93 plants with mutant phenotype, the nucleotide acid G2258 was mutated to A, while in all 95 plants with wild type phenotype, G2258 site was either wild type or heterozygous, with a ratio of 27:68 (close to 1:2), verifying that CqPORB gene mutation is genetically associated with the mutant phenotype of nl6-35 mutant.

In addition, we also crossed nl6-35 mutant with another quinoa variety Longli-4 that originates from Bolivia. In F₂ generation, the segregation of the faded green leaf phenotype is approximately 3:1 ( Table 1 ), corroborating that the mutant phenotype in nl6-35 was controlled by a recessive mutation. Similarly, 180 plants with faded green phenotype were selected for KASP analysis. The result showed that all these plants with mutant-type phenotype carried CqPORB allele that were derived from nl6-35 but not from Longli-4 variety, which further supported that CqPORB gene mutation is associated with the faded green leaf of nl6-35 mutant.

To further demonstrate that the reduced chlorophyll content in nl6-35 mutant was caused by the CqPORB gene mutation, we generated protoplasts from both wild type NL-6 and nl6-35 and transformed wild type CqPORB to the protoplasts of nl6-35. In protoplast, total Chl content in the nl6-35 mutant was much lower than that of the NL-6. However, transformation of CqPORB gene, but not the empty vector, to the protoplasts of nl6-35 significantly increased Chl content ( Supplementary Figure 3 ), suggesting that CqPORB gene mutation is responsible for the reduced Chl content in the nl6-35 mutant.

Studies have shown that disruption of POR genes affects the development of thylakoids in both Arabidopsis and rice (Paddock et al., 2010; Paddock et al., 2012; Sakuraba et al., 2013). Here we took a close view of chloroplast in nl6-35 by using transmission electron microscopy. In the wild type NL-6, numerous grana stacks were observed in the chloroplasts of light-grown plants. In the nl6-35 mutant, however, only unstacked thylakoids were observed ( Figure 3A ), indicating that CqPORB is crucial for grana development in quinoa.

In Chl biosynthesis pathway, POR is an essential enzyme that catalyzes the reduction of Pchlide to Chlide in a light-dependent way (Tanaka and Tanaka, 2007). POR-bound Pchlide is termed ‘photoactive’ and be able to converted to Chlide instantaneously in the presence of light, whereas Pchlide not bound by POR is termed ‘non-photoactive’ and has been indicated as a factor causing photooxidative damage (Paddock et al., 2010). Here, we measured the content of total Pchlide and Chlide in NL-6 and nl6-35 before and after light exposure for 10 s, 30s, and 5 min. Under dark conditions, the total Pchlide content was comparable between wild type NL-6 and nl6-35. However, after illumination, the photoactive Pchlide was reduced in a higher speed in NL-6 than in nl6-35 mutant, while Chlide was increased in a higher rate in NL-6 than in nl6-35 mutant, suggesting that CqPORB is important for the conversion of Pchlide to Chlide ( Figure 3B ).

The genomic length of CqPORB annotated in NL-6 genome was 2780 nucleotide acids, and this gene harbors five exons and encodes a protein of 395 amino acids. CqPORB shares high sequence similarities with its homologs in other plants species, including Arabidopsis thaliana At5g54190/PORA (77.5%), At4g27440/PORB (76.6%) and At1g03630/PORC (76.6%); Oryzae sativa Os04g0678700/PORA (70.3%) and Os10g0496900/PORB (72.8%); Physcomitrella patens Ppls146_112V6.1 (78.2%) and Ppls108_171V6.1 (68.2%); Selaginella moellendorffii SMpep_g142334 (71.0%); and Chlamydomonas reinhardtii Cre01.g015350.t1.1 (65.8%) ( Supplementary Figure 4 ). This result suggested that POR proteins are conserved from algae to higher plants.

Alignment of CqPORB against the protein database of NL-6 genome revealed that there are two additional isoforms, CqNL-6_007045 and CqNL-6_022507, which are located on Chr4 and Chr12, respectively, based on the annotation of NL-6 genome. Quinoa is an allotetraploid plant that consists of A sub-genome and B sub-genome that originated from diploid Chenopodium pallidicaule and diploid Chenopodium suecicum, respectively (Jarvis et al., 2017). Based on synteny analysis of quinoa genome, it has been shown that Chr4 and Chr12 are the colinear chromosomes in quinoa that originated from the A and B sub-genomes of its progenitors (Jarvis et al., 2017). These clues suggest that CqNL-6_007045 and CqNL-6_022507 are evolutionarily the same isoform, and thereby they were designated as CqPORA (CqNL-6_007045) and CqPORA-like (CqNL-6_022507), respectively.

Interestingly, we only identified one copy of CqPORB on Chr17 in NL-6 variety that originated from A sub-genome, indicating that the copy of CqPORB might be lost in the B sub-genome. To understand whether this is also the case in other quinoa varieties, we aligned CqPORB protein sequence with the protein database of the sequenced QQ74 variety (Jarvis et al., 2017). Unexpectedly, there are four CqPOR isoforms in QQ74 variety, and the two copies of CqPORA and CqPORB were identified in the A sub-genome and B sub-genome, respectively ( Supplementary Table 4 ). These results indicated that the CqPORB copy was only lost in the NL-6, but not in QQ74. To understand whether CqPORB-like gene was truly lost in the genome of NL-6 or it was caused by the incomplete genome sequencing. The collinearity of the CqPORB-like gene and its flanking genes between QQ74 and NL-6 varieties was analyzed. Totally 20 genes, including CqPORB-like gene, were selected in the QQ74 genome, and among them 9 genes were also identified in the collinear region of NL-6 genome, but CqPORB-like gene and its neighboring genes was completely lost in the NL-6 genome ( Figure 4A ). To further exclude the possibility that the loss of CqPORB-like gene in NL-6 was caused by incomplete sequencing, primers that matched to both CqPORB and CqPORB-like genes were used for gene amplification using the genomic DNAs isolated from NL-6, Longli-4 and Real varieties, and PCR products were sequenced by sanger sequencing. The result showed that only one isoform was amplified in the NL-6 variety, but two isoforms were amplified in both Longli-4 and Real varieties ( Figure 4B ), verifying that CqPORB-like gene was lost in the NL-6 genome.

In a previous study, we have performed RNA-seq analysis for four different tissues, including leaf, root, stamen, and pistil in NL-6 variety (Jiang et al., 2022). Using these transcriptomic data, we found that all the three CqPOR genes were highly expressed in the leaves, but with a much less level in pistil, and they were almost not detected in root and stamen ( Figure 5A ). This gene expression pattern is consistent with the role of POR genes in the regulation of Chl biosynthesis in photosynthetic tissues. Moreover, the transcriptomic data showed that CqPORA and CqPORA-like genes rather than CqPORB gene exhibited a more similar gene expression pattern ( Figure 5A ).

It has been known that the expression of POR genes is altered upon illumination (Armstrong et al., 1995). Here, we tested the expression of these three CqPOR genes in NL-6 etiolated seedlings during the transit from dark to white light. qRT-PCR analysis showed that both CqPORA and CqPORA-like genes were down-regulated, while CqPORB gene was up-regulated after exposure to light ( Figure 5B ). The similar expression pattern of CqPORA and CqPORA-like genes in response to light further supported that they are evolutionarily the same isoform. The different expression pattern between CqPORA and CqPORB genes indicated that these two genes may exhibit a distinct function in the regulation of Chl synthesis in response to ambient light signal. We also analyzed the expression of these three POR genes in different developmental stages in quinoa, and found that after growth on soils for 10, 30, and 50 d, the expression levels of all these three POR genes were gradually decreased, but the CqPORB gene was reduced in a slower rate than the CqPORA genes ( Figure 5C ). At later developmental stages, we noted that the transcript level of CqPORB was much higher than the CqPORA and CqPORA-like genes ( Figure 5D ), indicating that CqPORB plays a predominant role in Chl biosynthesis and chloroplast development in mature plants.

To study the subcellular localization of CqPORB in quinoa, we infiltrated Agrobacterium strain harboring 35S::CqPORB-GFP plasmid into quinoa leaves according to the protocol described previously (Xiao et al., 2022). However, after transformation for 48 h, no GFP fluorescence signal was detected in quinoa leaves, suggesting that the efficiency of transient gene expression in quinoa leaves is still very low. We then transformed 35S::CqPORB-GFP plasmid to tobacco leaves, and found that CqPORB was exclusively localized in chloroplast ( Figure 6A ). Besides, 35S::CqPORB-GFP was transformed into Arabidopsis to obtain stable transgenic plants, and in these transgenic plants, GFP fluorescence was also clearly observed in chloroplast ( Figure 6B ).