The F₂ population (AM1588) used in this study originated from the self-pollination of an F₁ plant (CM8996-199) previously developed for genetic mapping of resistance to whiteflies (Bohorquez-Chaux et al., 2025). The two parents from which the F₁ originated (ECU72 and COL2246) exhibited contrasting flowering phenotypes, with ECU72 male sterile, barely producing any male flowers. The CM8996–199 F₁ however, produced both male and female flowers. The mapping population consisting of 109 individuals was first evaluated in 2020, then in 2025. This population was established at the International Center for Tropical Agriculture (CIAT) in Palmira, Colombia in a single diagonal arrangement with five plants per genotype. While the specific field plots differed between years to adhere to crop rotation practices, they were located within the same station featuring similar soil characteristics and environmental conditions. Similar agronomic management regimes were applied in both years. For validation, progenitors comprising 304 diverse landraces and breeding clones from CIAT breeding program were utilized.

Flowering data were collected for all clones of the mapping population at 4 months after planting (MAP) in 2020 and four consecutive time points in 2025: 6, 7, 8, and 9 MAP. Phenotyping in 2025 was performed blind to genotype identity since scoring was decoded post-evaluation. At each time point, we recorded the independent presence of male and/or female flowers separately. This detailed descriptive data was then synthesized into a categorical trait score for each plant, on a 0–2 scale as follows: 0 = absence of flowering; 1 = onset of flowering (without clarity on whether the resulting flower will be either male or female); and 2 = presence of flowering (Supplementary Figure S1). To obtain the score of each accession, Best Linear Unbiased Predictions (BLUPs) adjusted by the grand mean were calculated from the five individual plant scores. The accessions and their corresponding flowering scores are summarized in Supplementary Table S1. Spearman’s correlation coefficient was applied to estimate correlations among the different time points (MAP). Phenotyping of the validation population was done at 6 MAP.

Using the previously constructed linkage map on this AM1588 F₂ mapping population (Bohorquez-Chaux et al., 2025), QTL mapping was performed using BLUPs for five datasets, corresponding to flowering evaluated at 4 (2020), 6, 7, 8, and 9 MAP (2025). Composite interval mapping (CIM) (Zeng, 1994) was performed in WinQTLCart 2.5 (Wang et al., 2007), with significance thresholds determined by 1,000 permutation tests (α = 0.05) (Churchill and Doerge, 1994). Analyses were run with a 10 cM window, model 6, a 1 cM walk speed, and five marker cofactors selected via forward–backward regression. Candidate genes within 2-LOD intervals of significant QTLs were identified using the Manihot esculenta v6 genome (https://phytozome-next.jgi.doe.gov/info/Mesculenta_v6_1) which the genetic map was based on. From the list of genes in these intervals, candidate genes were prioritized based on functional annotations and literature searches for homologs of known flowering-time regulators, previously characterized in model species.

Markers closest to the peak regions were selected for validation. Since the region on chromosome 13 was highly significant, we identified three other markers around the peak region. The significant markers were validated in the cassava progenitors (N = 304) exhibiting variable flowering. The allele frequencies and quality of these single nucleotide polymorphisms (SNPs) were assessed using a metric that estimates false positive (FPR) and false negative (FNR) rates (Bohorquez-Chaux et al., 2025; Mbanjo et al., 2024; Platten et al., 2019). Favorable alleles in the homozygous state were identified in genotypes with a score of 2 (presence of flowering). In contrast, unfavorable alleles in the homozygous state were identified in genotypes with a score of 0 (absence of flowering). The validation population and corresponding flowering scores and marker data are summarized in Supplementary Table S2.

For the 4 MAP evaluation in 2020 and the 6, 7, 8 and 9 MAP evaluation in 2025, the presence of male and female flowers was perfectly correlated (in all cases, r² > 0.99, p < 0.001); whenever flowering occurred, both male and female flowers were present. Therefore, we analyzed the trait as general flowering. The calculated correlations (r²) among the different months were between 0.54 - 0.95 with the least correlation between 4 MAP and 9 MAP, and the highest observed correlation between 8 MAP and 9 MAP (Supplementary Figure S2). At 9 MAP, when evaluations were completed in the AM1588 F₂ population, 62 genotypes (67.6%) exhibited both female and male flowers (score 2), seven genotypes (7.63%) showed the onset of flowering (score 1), and 40 genotypes (43.6%) displayed no flowering (score 0). The validation population, comprising 304 parental genotypes, was evaluated at 6 MAP using the same 0–2 scale applied to the AM1588 F₂ population. Among the 304 parental genotypes evaluated, 148 (48.7%) exhibited both female and male flowers (score 2), 39 (12.8%) showed the onset of flowering (score 1), and 117 (38.5%) exhibited no evidence of flowering (score 0) (Supplementary Table S2).

For the five traits mapped in the AM1588 F₂ population, QTL were detected on chromosomes 1, 7, 13 and 16, with phenotypic variance explained (R²) values ranging from 5.45% to 42.63% (Figure 1, Table 1). Among these, the QTL on chromosome 13 (Meflwr13) was consistent and overlapped across all five traits, while Meflwr1 and Meflwr7 were present in four and three traits, respectively. Meflwr16 was only identified at 8MAP. Meflwr13 exhibited the highest LOD scores (maximum LOD = 20.82), and phenotypic variance explained (R² = 42.63%). Table 1 summarizes the detected QTL for each trait along with their corresponding LOD scores and R² values.

To validate the QTL identified in the analysis, selection focused on SNPs located at or nearest to the peak regions. Figure 2 illustrates the marker-trait association across the candidate genomic regions. For the minor QTL regions, one representative SNP was selected from the peak area of chromosomes 1, 7, and 16. As shown in Figure 2A, these markers exhibited lower association with the trait. Three SNPs were selected closest to the peak of the major QTL on chromosome 13 (Meflwr13) (Figure 2B), along with a further three markers selected based on their proximity to the peak markers (Figure 2C), allowing for a finer view of the locus. The analysis revealed that the SNPs on Meflwr13 exhibited the strongest association with the trait, displaying low p values consistent with a major effect locus, in contrast to the SNPs on the other QTL. These results confirm the high potential of this locus for marker-assisted selection (MAS).

The six Meflwr13 SNPs were further evaluated for their potential use in MAS by checking their allele frequencies and phenotypic association in the breeding progenitors of diverse backgrounds (Figure 3). The six SNPs exhibited dominant allelic effects as no significant difference was observed between genotypes homozygous for the favorable allele and heterozygotes. In contrast, the homozygous unfavorable genotype (phenotypic score of 0) was significantly distinct in all cases, indicating that the unfavorable allele could be effectively selected against during MAS. Notably, markers C13_634483 and C13_658450 showed the strongest differentiation between homozygous genotypes carrying the favorable and those with the unfavorable alleles (p ≤ 0.0001) in the progenitors. A total of 78.2%, 64.4%, and 76.8% of the evaluated genotypes carried at least one copy of the favorable alleles (C, G, and G) for markers C13_889929, C13_634483, and C13_658450, respectively. In contrast, for C13_643226, C13_743917, and C13_620577, only 45.0%, 33.4%, and 34.5% of the genotypes carried at least one copy of the favorable alleles (G, T, and A, respectively). C13_889929, C13_634483, and C13_658450 showed acceptable FPR and FNR values, while the other three had very high FNR values (Table 2). Based on validation in both the AM1588 F₂ mapping population and the progenitor validation population, these three SNPs were identified as the most promising markers for selection.

Further examination of Meflwr13 was conducted to pinpoint potential candidate genes (Supplementary Table 3). Because the locus mapped to the proximal end (start) of the chromosome, the research focused on the region spanning from 0 Mb to the right-most flanking marker, located at approximately 2.5 Mb (2,495,605 bp). Flowering-related genes were identified within this region (Figure 4), including two major genes of interest: Manes.13G000800 (260 kb) which encodes a FLOWERING LOCUS T (FT) gene and Manes.13G011900 (1.12 Mb) which encodes a TERMINAL FLOWER 1 (TFL1) gene. The SNP explaining the highest phenotypic variation in the mapping population was C13_889929_G. This SNP mapped on an intron of a protein of unknown function (DUF1664) and is located only 6 kb from Manes.13G009100, a WRKY transcription factor (TF). WRKY TFs are implicated in modulating flowering pathways to regulate flowering time (Song et al., 2024). Interestingly, an Apetala-like ethylene-responsive transcription factor {AP2/ERF TF: Manes.13G001800 (376 kb)} and TEOSINTE BRANCHED 1 {TB1: Manes.13G008300 (815 Mb)} were also identified in the region. AP2/ERF TFs are known to repress flowering by direct repression of FT expression (Sgamma et al., 2014; Yant et al., 2010), in addition to their involvement in determining the identity of flower organs (Luo et al., 2024; Wang et al., 2023). TB1 also interacts with FT to delay flowering (Colleoni et al., 2024; Feng et al., 2022; Mimida et al., 2011). On Meflwr1, an APRR5 response regulator of two component system (Manes.01G043200) that controls flowering time was also identified in the peak region (7.05 Mb). Genes within the various identified QTL that have been previously reported as being involved in flowering time, male sterility/fertility, or gametophyte and pollen development in other species are summarized in Table 3.

The development of precise phenotyping methods is foundational to genetic studies. This work utilized a categorical 0–2 scoring scale for the absence, onset and presence of flowering, providing a simple yet effective tool to quantify the vegetative-to-reproductive transition in the AM1588 F₂ population. Flowering scales have been utilized across diverse plant species to quantify this transition. For instance, in model species like Arabidopsis, detailed, multi-stage scales are used to precisely link visible flowering stages to the molecular expression profiles of key floral regulators (FT, FLC, SOC1) (Boyes et al., 2001). In barley developmental flowering stages determined are used in accordance with the Zadoks’ scale (Z55) (Parrado et al., 2024; Zadocks et al., 1974), while soybean defines flowering by the R1–R2 stages (beginning to full bloom) (Fehr et al., 1977; Plumblee and Harrelson, 2022). In the context of cassava, a presence and absence score of 1 and 0, respectively, has been used for flowering phenotyping (Fukuda et al., 2010). Our study employed a categorical 0–2 scoring scale to enhance resolution of the flowering phenotype for quantitative trait locus (QTL) mapping.

The time-series phenotyping provided crucial insights into the trait’s dynamics. The high correlation of 0.95 observed between later stages (8 MAP and 9 MAP) indicates that the plant’s final flowering status is largely established by eight months. Conversely, the lower correlation (0.54) between the earliest (4 MAP) and latest time point suggests that the genetic or environmental control governing early flowering initiation may be partially distinct from the factors regulating the established, later flowering phenotype, or more sensitive to environmental influence. In our analyses, there was perfect correlation between the presence of male and female flowers. The F₁ parent of the mapping population was fertile, with both male and female flowers, despite having a male sterile donor (Bohorquez-Chaux et al., 2025). This may suggest that the male sterility trait is likely controlled by a recessive gene but remains to be investigated.

Efforts to identify genomic regions associated with flowering in cassava have so far focused primarily on the measurement of flowering-related traits, such as branching type and level, number of nodes, height to first branch, and the numbers of pistillate and staminate flowers (Baguma et al., 2024; Boonchanawiwat et al., 2011; Ewa et al., 2021; Zhang et al., 2018a). However, unlike previous efforts that focused on secondary traits or proxies (such as branching architecture), this study represents the first reported QTL mapping analysis directly based on a categorical flowering time scale (the 0–2 scoring system).

This study successfully identified and validated QTL associated with flowering in cassava, providing critical genomic resources for accelerating breeding efforts in this staple crop. Our results confirm the quantitative and polygenic nature of flowering time in cassava, consistent with other crops, while highlighting a major-effect locus on chromosome 13 (Meflwr13) as the primary determinant of flowering variation within the AM1588 F₂ population. The identification of Meflwr13, which consistently exhibited the highest LOD score and explained the largest proportion of phenotypic variance (up to 42.63%) across all five evaluation time points (4, 6, 7, 8, and 9 MAP) was key. The stability and magnitude of this QTL strongly suggest it harbors major genes controlling flowering. Besides Meflwr13, we also identified other QTL for flowering on chromosomes 1, 7 and 16. Meflwr1 was consistent across four time points, indicating a role in general flowering initiation, albeit with a smaller phenotypic variation explained (up to 13.26%). The region harbors an APRR5 response regulator homolog (Manes.01G043200), which is part of the circadian clock. It is known to regulate flowering time by influencing CO/FT expression (Hayama et al., 2017; Sato et al., 2002). Meflwr7 was mapped at 4, 6 and 7 MAP with R² values up to 13.18%, suggesting it modulates early-to-mid stage reproductive development before its influence diminishes in later stages. Meflwr16 was only mapped at one time point (8MAP) suggesting it is likely transient or more stage-specific, and may function in a more environment-specific manner.

Meflwr13 QTL encompasses two genes central to the global floral network: Manes.13G000800, a homolog of FLOWERING LOCUS T (FT) and Manes.13G011900, a homolog of TERMINAL FLOWER 1 (TFL1) (Figure 4). Both FT and TFL1 belong to the phosphatidylethanolamine-binding protein (PEBP) family and are highly homologous (Ahn et al., 2006; Colleoni et al., 2024; Nadal Bigas et al., 2025). Despite this structural similarity, they perform opposite functions in plants (Azevedo et al., 2025; Jin et al., 2021; Kaneko-Suzuki et al., 2018; Liu et al., 2021; Matsoukas et al., 2012). FT acts as a florigen, promoting flowering, while TFL1 is a floral repressor by antagonizing the activity of FT. This antagonistic relationship is critical in regulating the transition from vegetative growth to reproductive phase. The proximity of these genes in the QTL—and the functional relationship between the FT and TFL1 orthologs—indicates a high likelihood of functional interactions between these genes in the region to control the flowering phenotype. This co-localization supports previous functional studies in cassava that demonstrated the central role of FT homologs in controlling flowering and branching (Adeyemo et al., 2011, 2017, 2019).

This regulatory complexity is further supported by the inclusion of two key transcription factors that control the core FT/TFL1 switch: Manes.13G001800, an AP2-like ethylene-responsive transcription factor (AP2/ERF TF), and Manes.13G008300, which encodes a TEOSINTE BRANCHED 1 (TB1) gene. The AP2/ERF TFs are known to act as transcriptional repressors of flowering, often directly inhibiting FT gene transcription in the leaves to reduce the florigen signal (Sgamma et al., 2014; Yant et al., 2010). Similarly, TB1 is a founding member of the TCP (TEOSINTE BRANCHED 1/cycloidea/proliferating cell factors) family of transcription factors that are critical regulators of plant architecture. TB1 and related TCP proteins have been described to interact directly with FT in other crops including apples (Mimida et al., 2011), rye (Zhan et al., 2023), and Brassica juncea (Feng et al., 2022). Their primary mechanism is thought to involve disrupting the florigen activation complex, thereby delaying flowering.

Additionally, the region includes a WRKY75 transcription factor (Manes.13G009100), identified near the most significant SNP from QTL mapping. Various WRKY TFs have been identified to play a crucial role in modulating flowering pathways to regulate flowering time (Song et al., 2024). Specifically, WRKY75 orthologs have been repeatedly implicated in promoting flowering across different species. In Brassica juncea, BjuWRKY75 promotes flowering by activating FT expression (Feng et al., 2022). In the gibberellin pathway, AtWRKY75 binds to FT to promote its expression, resulting in early flowering (Zhang et al., 2018b). Furthermore, overexpression of CpWRKY75 from Chimonanthus praecox is involved in multiple pathways to promote flowering time in transgenic Arabidopsis (Huang et al., 2022). Studies established that both MdFT1 and MdTFL1 interact competitively with MdWRKY6 protein to facilitate and inhibit, respectively, activation of an apple LEAFY-like gene, ultimately regulating apple flower bud formation (Zuo et al., 2024). The presence and proximity of these antagonistic and promoting factors in Meflwr13 indicates a highly complex and integrated molecular switch driving the primary flowering phenotype.

The genetic background of the mapping population, derived from a male-sterile grandmother (ECU72), adds a layer of complexity. Although the F₁ parent was fertile, the subsequent perfect correlation between male and female flower presence in the F₂ is noteworthy. The identification of several genes related to sterility, pollen development, and fertility restoration within the Meflwr13 interval (Table 3)—such as PPR repeats, beta-ketoacyl reductase 1, and pectate lyase (Beaudoin et al., 2009; Bisht et al., 2025; Durand et al., 2021; Gaborieau et al., 2016; Liu et al., 2017; Smirnova et al., 2013; Zheng et al., 2018) —suggests a possible pleiotropic role for this locus. However, it remains to be determined whether this arises from true biological pleiotropy or tight linkage drag between flowering time regulators (e.g., FT/TFL1) and distinct fertility genes. This makes the region a strong candidate for future investigation into the genetic basis of male sterility and fertility restoration in cassava, which could further refine breeding strategies. The region could be further dissected through high-resolution fine-mapping by increasing marker density within the interval and screening a larger F₂ or F₃ population to identify recombinants that potentially uncouple these traits.

Marker selection focused on SNPs within the peak regions to confirm their association with the trait, with the Meflwr13 SNPs exhibiting the strongest statistical association (Figure 2). This high association was successfully confirmed by validation in the independent progenitor population (Figure 3), which is a crucial step toward practical application. The analysis revealed that three of the six key SNPs—Chr13_889929, Chr13_634483, and Chr13_658450—are the most robust candidates. They exhibited a dominant segregation pattern for the favorable (flowering) allele and showed the most acceptable False Negative Rate (FNR) values (Table 2). While these markers exhibited higher FPR compared to the excluded markers, they provide the most effective compromise for breeding, prioritizing the retention of favorable flowering alleles. In contrast, the alternative markers displayed FNR values exceeding 60%, which would result in the inadvertent discarding of the majority of favorable genotypes. Dominant markers are highly prized in breeding because they simplify the selection process. While the validation in the diverse progenitor population demonstrates the utility of these markers, their broader applicability across distinct genetic backgrounds and geographic locations warrants further investigation to ensure stability. Flowering is heavily influenced by environmental cues such as photoperiod, temperature, and altitude. It will be essential to determine if the Meflwr13 QTL effect remains stable or if significant genotype-by-environment interactions influence the efficacy of marker-assisted selection.

Current flower-induction technologies including photoperiod extension, branch pruning, and the application of growth regulators are highly labor-intensive and costly to implement. For instance, expenses may range from approximately USD 1,000–3,000 for a dozen progenitors in a national breeding program. Furthermore, the effectiveness of these manual methods has only been evaluated in a limited number of genotypes, underscoring the urgent need for robust, large-scale genetic tools (Barinas et al., 2023).

The successful validation of these dominant markers within the major Meflwr13 locus has immediate and significant translational implications. Developing molecular markers that enable early selection—at the nursery stage— is the ideal solution. Breeders can use these markers in marker-assisted selection (MAS) to rapidly and cost-effectively screen seedlings before they are transplanted to the field. Early flowering genotypes, confirmed by the presence of the favorable dominant Meflwr13 allele, can be immediately allocated to standard crossing lots, saving time and resources by eliminating the need for expensive flower-induction treatments. These resources can then be strategically redirected toward those genotypes expected to flower late, maximizing the efficiency of the entire breeding pipeline.