We used BLASTP (Johnson et al., 2008) to compare the Arabidopsis NADP-ME protein sequences with the cassava reference genome (AM560-2-JGI-v8.1). Genes with an E-value ≤ 1e−5 were selected, and their chromosomal locations were determined. The chromosomal map was created using the “Gene Location Visualize from GTF/GFF” tool in TBtools (Chen et al., 2020). The conserved protein domains of the MeNADP-ME gene family were analyzed using the MEME suite (https://meme-suite.org/meme/), setting the motif number to 15. Protein analysis, including amino acid length, isoelectric point, molecular weight, and hydrophilic coefficient, was done using the ExPASy server (https://web.expasy.org/protparam). The subcellular localization of MeNADP-ME proteins was predicted using the Plant-mPLoc server (http://www.csbio.sjtu.edu.cn/cgi-bin/PlantmPLoc.cgi). To study cis-regulatory elements, the 3,000-bp upstream sequence of each gene was submitted to the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

This study used cultivated cassava varieties KU50, TMS60444, SC205, and SC16 and wild varieties FLA4047 and W14. The cultivated varieties are widely planted and closely related to the reference genome AM560-2. TMS60444 is commonly used in tissue culture, while SC205 is frequently selected for research. KU50 exhibits superior resistance to abiotic stress compared to SC205, and SC16, developed by our team, shows a 40% higher yield than SC205. All plants were grown in Chengmai, Hainan, China (19.85°N, 110.08°E), with samples collected from storage roots, fibrous roots, leaves, and stems during the root expansion stage. Shading transcriptome analysis was performed using SC205 and SC16, chosen for their significant yield difference, hypothesizing divergent photosynthetic efficiencies. This approach helps eliminate variability in photosynthetic gene responses and construct a co-expression network to study NADP-ME’s role in cultivated cassava’s photosynthetic physiology. Shading treatment began at 8:00 a.m., with samples collected at 8:30 a.m. and 10:00 a.m., along with control samples. For circadian rhythm analysis, TMS60444 was sampled every 2 h from 6:00 a.m. to 8:00 p.m. Heat and shade-heat stress experiments were conducted with SC205 and KU50 to explore MeNADP-MEs responses and determine if they have physiologically diverged from the original C₃ wild types. Control samples were collected at 10:00 a.m., followed by heat treatment at 43°C, with samples taken after 24 h and 48 h. Each treatment had three biological replicates, and all samples were immediately stored in liquid nitrogen.

Total RNA from all tissue samples was extracted using the RNAprep Pure Kit (TianGen Biotech., Ltd., Beijing, China; DP441) and treated to remove DNA, following the manufacturer’s instructions. RNA concentration and quality were measured with a NanoDrop, and RNA integrity was confirmed by 1% agarose gel electrophoresis. Only samples with intact RNA bands and a high quality (A260/280 ratio between 2.0 and 2.2) were used for reverse transcription. High-quality RNA was combined with the 5X Evo M-MLV RT Reaction Mix Ver.2 kit (Accurate Biology., Ltd., Wuhan, China; AG11728) to synthesize single-stranded cDNA following the manufacturer’s protocol. For RNA-seq, RNA integrity was further verified using the RNA Nano 6000 kit (Agilent Technologies, Santa Clara, CA, USA; 5067–1511) on the Bioanalyzer 2100 system.

First, mRNA was enriched using mRNA Capture Beads, followed by magnetic bead purification. The mRNA was then fragmented under high-temperature conditions. The fragmented mRNA served as the template for first-strand cDNA synthesis in the reverse transcriptase reaction system. End repair and addition of an A-tail were completed during the synthesis of the second cDNA strand. Adaptors were ligated, and target fragments were selected using Hieff NGS^® DNA Selection Beads (Yeasen Biotechnology Co., Ltd., Shanghai, China; 12601ES). The library was constructed by PCR amplification and sequenced on the Illumina Novaseq X Plus platform. For data processing, raw sequencing data were quality-controlled using fastp with default parameters (Chen et al., 2018). Bowtie 2 was used to align the clean reads to the ribosomal RNA database, removing reads mapped to ribosomal RNA. Unmapped reads were kept for transcriptome analysis (Langmead and Salzberg, 2012). Then, HISAT2 was used to align the paired-end reads to the cassava reference genome (AM560-2-JGI-v8.1) to generate read counts for each gene (Kim et al., 2019). Using the reference genome annotation file, the exon lengths of cassava genes were extracted, and FPKM values were calculated for each gene (Trapnell et al., 2010). The raw data for each sample were uploaded to the National Genomics Data Center, with quality control metrics and mapping rates provided in ( Supplementary Table S1 ).

For quantitative real-time PCR analysis, cDNA samples were diluted 20-fold. The expression levels of target genes were measured on a Bio-Rad CFX96 instrument using the 2X SYBR Green Pro Taq HS Premix II kit (Accurate Biology, Ltd., Wuhan, China; AG11702), with each sample analyzed in four technical replicates. The PCR cycling program began with a 30-s denaturation at 95°C, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Melt curve analysis and cycle threshold (Ct) values were determined using CFX software (Bio-Rad). The cassava Tubling gene was used as the reference gene for relative expression analysis, as described in previous studies (Yao et al., 2014). Primers for MeNADP-ME1, MeNADP-ME2, MeNADP-ME3, and MeNADP-ME4 were designed across exons from the CDS sequence and verified for specificity using Primer-BLAST (Ye et al., 2012). Primer sequences are provided in ( Supplementary Table S2 ). Gene expression levels were calculated using the 2^−ΔΔCt method (Livak and Schmittgen, 2001).

We downloaded the NADP-ME amino acid sequences from Arabidopsis thaliana (TAIR10) and cassava (AM560-2-JGI-v8.1) and analyzed conserved domains using the MEME suite (https://meme-suite.org/meme/), setting the motif number to 15. To construct the phylogenetic tree, BLASTP was used to identify homologous genes (E-value ≤ 1e−5) in Arabidopsis thaliana (TAIR10), rice (MSU-v7.0), maize (MaizeGDB-Zm-B73-REFERENCE-NAM-5.0), potato (MSU-v6.1), and sorghum (JGI-v2), and their protein sequences were retrieved. Homologous NADP-ME genes and sequences from Flaveria robusta (C₃), Flaveria floridana (C₃–C₄), and Flaveria bidentis (C₄) were based on a previous study (Taniguchi et al., 2021). The phylogenetic tree was constructed using the maximum likelihood (ML) method. Multiple sequence alignment was done using the Muscle tool in MEGA v. 11 (Tamura et al., 2021), followed by model selection with default parameters and 1,000 bootstrap replicates. The phylogenetic tree was visualized using iTOL (Letunic and Bork, 2024).

To construct the pG1300-MeNADP-MEs-GFP vector, we used the pCAMBIA1300-GFP vector. Primers were designed to amplify the CDS of each MeNADP-ME family member from the cDNA of mature leaves of the cassava cultivar SC205. CDS cloning was performed using 2X ApexHF CL PCR Master Mix (Accurate Biology, Ltd., Wuhan, China; AG12209). The pG1300-GFP vector was digested with KpnI-HF^® and SalI-HF^® (New England Biolabs). Both the cloned CDS and the digested vector were purified using 1% agarose gel electrophoresis. The ClonExpress Ultra One Step Cloning Kit V2 (Vazyme, Nanjing, China; C116) was used for homologous recombination, ligating the CDS fragments into the vector. The recombinant vectors were transformed into Agrobacterium tumefaciens strain GV3101. After reaching an OD600 = 0.9, the bacterial cells were resuspended in MES buffer (10 mM MgCl₂, 10 mM MES, 100 µM AS) to OD600 = 9.0 and injected into the abaxial side of tobacco leaves. After 48 h in the dark, 1-cm² sections of the transformed leaves were sampled for imaging on a FV3000 confocal microscope (Olympus, Japan), capturing bright field, chloroplast fluorescence (680 nm), and GFP fluorescence (488 nm) images, which were merged using FV31S-DT software (Olympus, Japan).

Genes with the top 7,000 average FPKM values were selected and normalized using log2(FPKM+1). Sample clustering and soft threshold analysis were performed ( Supplementary Figure S1 ). A soft threshold of 18 was chosen to construct the weighted gene co-expression network analysis (WGCNA) co-expression network using the R package version of WGCNA (Langfelder and Horvath, 2008). Photosynthesis-related gene sets were created by extracting cassava genes involved in light-harvesting antennae (KEGG pathway ID: map00196), photosynthetic electron transport (KEGG pathway ID: map00195), and carbon fixation (including the Calvin–Benson cycle and C₄ genes; KEGG pathway ID: map00710) from the KEGG pathway database. The number of photosynthesis-related genes in each module was calculated, and the module with the highest number of photosynthetic genes was selected for further analysis. KEGG enrichment analysis was performed on all genes in this module ( Supplementary Figure S2 ). Using the Plant TF Motifs Shift tool in TBtools, protein sequences in the module (AM560-2 JGI-v8.1) were analyzed with Arabidopsis thaliana as a reference to identify transcription factor binding motifs. The 3,000-bp upstream promoter regions of photosynthesis-related genes were extracted and scanned using the Fimo: Binding Motif Scan tool in TBtools to generate a predicted transcription factor–promoter interaction network. The network was visualized using Cytoscape 3.7.2 software.

Each qPCR sample was analyzed in four technical replicates. Tukey HSD significance analysis of the qPCR results was conducted using SPSS software (IBM Corp., Armonk, NY, USA; version 8.0). Differences were indicated by letter labeling, where different letters represent significant differences (p < 0.05). Figures were created using GraphPad Prism software (GraphPad Software, San Diego, CA, USA; version 8.0).

A BLASTP search using Arabidopsis NADP-ME (TAIR10) as a query against the cultivated cassava AM560-2 dataset (Phytozome 13, E-value ≤ 1e−5) identified four distinct members of the MeNADP-ME family: MeNADP-ME1, MeNADP-ME2, MeNADP-ME3, and MeNADP-ME4, located on chromosomes 16, 16, 11, and 4, respectively ( Figure 1A ). Physicochemical analysis ( Table 1 ) revealed that MeNADP-ME1 has the fewest amino acids, lowest molecular weight, and smallest isoelectric point, while MeNADP-ME2 exhibited the highest values for these properties. MeNADP-ME3 and MeNADP-ME4 had intermediate characteristics. Subcellular localization predictions indicated chloroplast localization for all four isoforms. Promoter analysis ( Figure 1B , Supplementary Table S3 ) identified numerous light-responsive cis-regulatory elements, particularly in MeNADP-ME1 and MeNADP-ME3, suggesting their involvement in photosynthetic regulation. Additionally, these genes contain hormone- and stress-responsive elements, although light-responsive elements dominate their regulatory mechanisms. These findings suggest that the MeNADP-ME family plays a role in the photosynthetic physiology of cultivated cassava.

We conducted a MEME analysis to investigate conserved motifs in the NADP-ME proteins of cassava and Arabidopsis ( Figure 2A ). The results showed that the primary structural differences between the two species’ NADP-ME proteins occur mainly in the N-terminal region. Specifically, MeNADP-ME2, MeNADP-ME1, and MeNADP-ME4 shared high homology with similar Arabidopsis subtypes, while MeNADP-ME3 exhibited a unique N-terminal domain, lacking Motif15 compared to MeNADP-ME4 and AtNADP-ME4. We also constructed a phylogenetic tree with NADP-ME genes from species such as maize (C₄), Arabidopsis (C₃), Flaveria species, potato (C₃), rice (C₃), sorghum (C₄), and cassava (C₃–C₄) ( Figure 2B ). The analysis revealed four distinct subgroups: monocot, both monocot and dicot, cytosolic dicotyledonous, and plastidic dicotyledonous types. The phylogenetic distances correlated with both evolutionary relationships and subcellular localization. All MeNADP-ME members were classified into the plastidic dicotyledonous subgroup, consistent with localization predictions. Notably, cassava NADP-ME genes were most closely related to AtNADP-ME4, providing insights into the functional adaptation of these genes across species.

To investigate the subcellular localization of the MeNADP-ME gene family in cultivated cassava, we cloned the CDS regions of the four identified members from the AM560-2 reference genome, using TMS60444 as a template, and ligated them into a GFP vector driven by the CaMV 35S promoter. Transient expression in tobacco leaves ( Figure 3A ) revealed chloroplast localization for all MeNADP-ME members, consistent with phylogenetic and predictive analyses. Expression patterns in various tissues of cassava cultivar SC205 and wild species A4047 ( Figure 3B ) showed that in SC205, MeNADP-ME1 was most expressed in mature leaves, MeNADP-ME2 in young leaves, MeNADP-ME3 in mature leaves, and MeNADP-ME4 in young leaves and stems. In A4047, MeNADP-ME1 was highly expressed in mature leaves and tuberous roots, while MeNADP-ME2 and MeNADP-ME3 were highest in mature leaves. MeNADP-ME4 exhibited significantly higher expression in tuberous roots. These results highlight significant differences in the expression patterns of MeNADP-ME family members between cultivated and wild cassava, suggesting that the family has been selected to enhance leaf photosynthetic function in cultivated cassava, reflecting C₃–C₄ intermediate photosynthesis evolution.

To investigate the diurnal response of the MeNADP-ME gene family in cultivated cassava and identify key members involved in photosynthetic activity, we analyzed their expression dynamics at different time points in the cultivar TMS60444 ( Figure 4A ). The expression patterns of MeNADP-ME1 and MeNADP-ME2 were similar, showing low levels in the morning, increasing after noon, and peaking at night. In contrast, MeNADP-ME3 and MeNADP-ME4 exhibited a strong positive correlation with light intensity, particularly between 6:00 and 8:00 a.m., when light intensity sharply increased, leading to a peak in expression. After 10:00 a.m., expression gradually decreased, with MeNADP-ME4 showing a second peak at 20:00. Given the involvement of NADP-ME genes in C₃ species’ response to abiotic stress, we further examined the family’s response to heat and combined heat and shading stress in cassava cultivars KU50 and SC205 ( Figure 4B ). In KU50, gene expression was generally upregulated at PG24h, with MeNADP-ME2 significantly upregulated at PG48h, while other genes were downregulated. Under shade-heat stress, all genes, except MeNADP-ME4, were downregulated. In SC205, MeNADP-ME1, MeNADP-ME2, and MeNADP-ME3 peaked at PG24h, significantly decreasing by PG48h. Under shade-heat stress, MeNADP-ME2 was upregulated at SPG24h. In conclusion, MeNADP-ME3 and MeNADP-ME4 were identified as key contributors to photosynthesis in response to increased light intensity. Heat stress induced upregulation of MeNADP-ME genes, peaking at 24 h, while no significant changes occurred at 48 h or under combined heat and shading stress. These findings suggest that the MeNADP-ME gene family in cultivated cassava reflects an evolutionary trend toward C₄ photosynthesis while retaining characteristics of its C₃ ancestral species, underscoring the unique physiological features of its intermediate C₃–C₄ MeNADP-ME genes.

We conducted co-expression network analysis of the transcriptome under shading treatment in the cultivated cassava varieties SC16 and SC205, identifying 18 modules with distinct expression patterns. Notably, the Turquoise module, which included all members of the MeNADP-ME gene family, was enriched with photosynthesis-related and C₄ photosynthesis genes ( Figures 5A, B ). KEGG analysis revealed significant enrichment in the “carbon fixation in photosynthetic organisms” pathway ( Supplementary Figure S2 ). Under shading conditions, cassava C₄ genes, photosynthetic electron transport chain genes, and photosynthetic carbon assimilation genes, including the MeNADP-ME family, were regulated by transcription factors such as GATA1 (Manes.03G154500), YABBY1 (Manes.02G035700), Myb-like (Manes.15G163100), HY5 (Manes.12G040300), and MADS (Manes.10G099000) ( Figure 5C ). Heatmap analysis further showed that MeNADP-ME3 exhibited an expression pattern consistent with most photosynthesis-related genes and their transcriptional regulators ( Figure 6D ), while MeNADP-ME1, MeNADP-ME2, and MeNADP-ME4 exhibited opposite trends in SC16. This highlights MeNADP-ME3 as a key gene in the intermediate C₃–C₄ photosynthetic physiology of cultivated cassava. Transcription factor-promoter prediction suggested that MeYABBY1 regulates MeNADP-ME3 and binds to its systemic promoter Indel site, potentially playing a critical role in recruiting MeNADP-ME3 into the C₃–C₄ photosynthetic regulatory network ( Figures 6A–C ).