Identification of Genomic Regions Associated with Flowering Traits in Cassava (Manihot esculenta Crantz)

1. Introduction

Cassava (Manihot esculenta Crantz, 2n = 36) is a valuable food crop in fighting hunger and poverty in developing countries [1,2]. Additionally, it’s a source of animal feed as well as raw materials for a diverse industries [3]. Because of this economic importance of the crop, much of the ongoing breeding efforts are directed towards its genetic improvement so as to develop cassava varieties targeting different industry and food products [4,5]. This, however, is being hampered by a number of constraints including poor flowering and limited knowledge of the genetics and inheritance of flowering traits [5,6]. To meet the growing demand for cassava as well as overcoming the challenges to its genetic improvement due to poor flowering, it is imperative to expand knowledge on the genetics of flowering so as to customize solutions to this problem.

Cassava is preferentially propagated by vegetative means using stem cuttings [7,8] due to the difficulty in obtaining botanical seed, low germination and seedling establishment and survival rates which often result in small and weak plants with comparatively low root yields [8]. Despite these seed-related challenges, seed remains an indispensable and exceedingly important resource in cassava breeding pipeline, to generate new breeding lines with genetic variabilities which may be beneficial in breeding programs [7,8,9]. In conventional cassava breeding, genetic variabilities are derived through hybridization that involve crosses under controlled or open-pollinations, resulting into full or half sib progenies [8]. Introgression of economically important traits requires selection of progenitors with desirable traits and with timely and synchronized flowering capabilities. However, many of the progenitors with desirable farmers’ preferences are often poor at flowering while others never flower at all. Besides, there is great disparity in time and amount of flowering among many cassava genotypes, with some flowering early and profusely while in others the flowering is late and scarce [5,6,10]. The early flowering or branching types tend to have shorter main stem heights while the late branching ones have taller stem heights at first branching. All this indicates that cassava flowering is complex, exhibiting varied patterns, that make synchronization of crosses difficult [11].

Flowering in cassava begins with morphological changes in plant architecture that involve branching of the main stem [10]. Flower bud formation occurs at branch points in the stem [8,10] as axillary buds below the inflorescence allow upward growth of the plant. This implies, every branching event results in flower formation [12]. The successive branching may be di-, tri-, and tetra- or penta- chotomous, resulting in several branching events [8]. The highly branching genotypes flower more prolifically than those with sparse branching, while the non-branching (erect) ones are non-flowering. Most farmer-preferred cassava genotypes are either erect or non-branching or late flowering and thus rarely produce seeds [5,12]. The correlation between branching behaviour and flower production observed in cassava suggests that these traits might share a similar genetic basis. However, knowledge of genes controlling this kind of flowering correlation and/or behaviour in the crop is scanty.

Flowering process in angiosperms is evoked by an interaction of environmental and endogenous signals [13,14,15,16]. This is known to be regulated by a sophisticated network that monitors changes in the environment to ensure that flowering occurs under suitable conditions for maximization of reproductive success and seed production [17,18]. Like in other plants, studies have indicated that flowering response in cassava is under the influence of environmental factors, such as photoperiod and temperature [19,20,21]. The understanding of the control mechanisms of flowering in many plants is mainly based on the model plant Arabidopsis thaliana, in which various genes including FLOWERING LOCUS T (FT) and proteins have been implicated [17]. Although FT orthologs such as MeGI, Cassava GI and CO-like genes (MeCOL1, MeCO, and MeCOL2) have been associated with photoperiodic induction of flowering in cassava [21,22], the precise mechanism underlying the flowering behaviour exhibited in the crop, particularly branching, branch type, branching levels, number of nodes and stem height at first branching level, is not clearly known. Thus, deciphering the genetic basis of this flowering behaviour and how it is regulated will provide important insights that could be useful in overcoming hybridization challenges, especially flowering synchronization.

Molecular marker-assisted selection (MAS) has been a primary approach in plant breeding for understanding the genetic basis as well as discovery of functional genomic regions (genes) associated with target phenotypic traits in many plant species [23]. It has been used to probe QTL (quantitative trait loci) for developing genetic linkage maps in a number of crops [24]. Several QTLs for flowering time have been detected in several crop species, for example, rapeseed (Brassica napus L.) [25]. MAS has been used to detect QTLs or genomic regions for various traits in cassava breeding [23], for example, root yield and starch content [26], plant and first branch height [27], starch biosynthesis [28], biotic stress resistance [29]. However, knowledge of genes controlling flowering and related traits in cassava is still limited.

Genome-wide association analysis (GWAS) is an alternative genomic approach that has been designed to overcome drawbacks of QTL mapping [30]. It is used to mine natural genetic variations in organisms for revelation of variabilities in genetic architectures among traits in crop plants [31,32]. It is robust and thus has been used to identify many true genotype-phenotype associations [33] in several plants, for example, Arabidopsis thaliana [32,34,35], Oryza sativa [36], Zea mays [37,38], Brassica napus [39] and Setaria italica (foxtail millet) [40]. GWAS has also been successfully utilized to unravel genetic architecture of several phenotypic traits in cassava, for example, starch content [41], dry weight content [42], beta carotene content [42,43,44], disease-resistant traits such as cassava mosaic disease (CMD) [45] and cassava brown streak disease [46,47], and drought tolerance [48].

Relatedly, GWAS has been used to elucidate the genetic basis of flowering traits in several crops. For example; flowering time variabilities in pearl millet (Pennisetum glaucum L.) [49], maize (Zea mays) [50], common bean (Phaseolus vulgaris L.) [51,52], Sesame (Sesamum indicum L.) [53], rapeseed (B. napus L.) [54] flowering and plant height under different levels of nitrogen in Indian mustard (Brassica juncea) [55] and Soybean (Glycine max) [56]. In cassava, reports on utilization of GWAS to determine the genetic basis of flowering and related traits are scanty, except first branch height [57].

Thus, in this study, we aimed at identifying genomic regions associated with flowering traits, namely; first branching; type of branching; number of branching levels; number of nodes and stem height at first branching.

2. Results

2.1. Phenotypic Variation, Distribution and Heritability Estimates of Flowering Traits

Phenotyping was done for flowers (staminate and pistillate) and five flowering traits on a total of 293 cycle two (C2) cassava accessions. This was done for two growing seasons (2019 and 2020) and two locations. Branching ranged from 0 (did not branch) to 1 (branched), branch type or habit from zero to four, branching levels were maximum at 12, both number of nodes and stem height at first branching ranged from zero to 148 and 231cm respectively (Table 1). Coefficient of variation (CV) values varied from 26.87% for branching to 54.66% for nodes.

Broad-sense (H₂) and SNP (genomic) (h₂) heritabilities were estimated for each trait across season and location (Table 1). Highest broad-sense heritability (H₂) was exhibited by pistillates (0.83) followed by stem height at 1st branch (0.60), while staminates had the lowest heritability (0.00). SNP-based heritabilities ranged from lowest (0.00) in branch type to highest (0.35) in stem height (at 1st branch). Negatively skewed distribution trends were observed for branching and branching levels, while branch-type, nodes and stem height at first branching had positively skewed distributions (Figure 1). Results of ANOVA demonstrated that the differences among the cassava accessions (P ≤ 0.001) used in this study were highly significant for all traits, and magnitude was even higher for the seasons and locations, except branch type (Table 2).

2.2. Correlation between the Flowering Traits

Spearman’s correlation coefficient was applied to detect any correlations among the flower and flowering traits considered in this study; 1) first branching of the main stem, 2) branching type, 3) number of branching levels, 4) number of nodes, 5) stem height at the first branching, 6) flowers (pistillates and staminates). All the traits showed significant positive correlations ranging between 0.66 to 1.00 (Figure 2). Generally, relatively lower correlation values were noted between flowering traits and pistillates compared to other traits.

2.3. Genotyping and SNP Identification

Genotyping enabled the identification of 24,040 imputed SNP markers, all of which were mapped onto the 18 cassava chromosomes. Upon filtering, markers with minor allele frequency (MAF) at 1%, 22,103 SNPs were retained, distributed on all the 18 chromosomes (Figure 3). The SNP coverage per chromosome ranged from 806 SNPs on chromosome 16 to 2,065 SNPs on chromosome 1. The average values of the polymorphism information content (PIC) and MAF for the 22,103 SNP markers were 0.25 (ranging from 0.02 to 0.38) and 0.23 (ranging from 0.01 to 0.50) respectively.

2.4. Population Structure and Kinship

The population structure of the 293 accessions was estimated by using principal components analysis (PCA). The first principal component (PC) accounted for 37.3% of the observed variation (Figure 4A) while the second up to the fifth PC contributed 11.9%. Cumulatively, the first five PCs explained 49.2% of the genetic variation. Kinship among accessions estimated using a marker inferred Kinship matrix revealed familiar relationships between the accessions (Figure 4B). Both PCA and kinship analyses revealed clustering of germplasm into subgroups, three major subgroups according to the kinship heatmap (Figure 4B).

2.5. Linkage Disequilibrium Estimation

Linkage disequilibrium (LD) decay rate was measured as the chromosomal distance at which the average pairwise correlation coefficient (r²) dropped to half its maximum value. The LD and LD decay distance in all the 18 chromosomes is summarized in Figure 5. The average r² value for 18 chromosomes was found as 0.32 at a distance of ~2850kb (kilo base pairs) which slowly dropped to its half (0.16) at distance ~4828kb, with maximum decay distance of ~32012kb on chromosome 7 and minimum decay distance of 0.009kb on chromosome 10. The lowest r² value was 0.12 with maximum decay distance of ~39193kb on chromosome 8 and minimum decay distance of 0.006kb on chromosome 15. Meanwhile, the highest r² value was 0.95 with maximum decay distance of ~796kb on chromosome 3 and minimum decay distance of ~33kb on chromosome 10.

2.6. Marker-Trait Association Mapping

Marker-trait association (MTA) using 293 cassava accessions and 22,103 SNP markers was undertaken using the Mixed Linear Model (MLM), which took genetic relatedness into account, greatly reduced false positives, as shown in the Manhattan and quantile–quantile (Q-Q) plots (Figure 6). The Manhattan plots of five flowering traits were based on the combination of >1000 trait de-regressed BLUPs. A total of 53 significant MTA (involving 27 SNP) was identified as being associated with flowering traits, and these were distributed mainly on chromosome 18 which had 20 of the 53 SNPs (Table 3). Branching (Branch) had 15 marker associations on chromosomes 5, 8, 11, 15, 16 and 18; branch type at first branching (Branch1_No) was associated with 08 SNPs on chromosomes 5, 8, 14, 15 and 18. Similarly, number of branching levels (Branch_Levels) had 14 marker associations on chromosomes 5, 8, 9, 11, 14, 15, and 18., while number of nodes at first branching (Branch1_Nodes) had the highest associated markers (16) on chromosomes 3, 5, 8, 9, 14, 15 and 18. Meanwhile, stem height at first branching (Branch1_Ht) had no MTA based on the Bonferroni threshold (Figure 6). The Q-Q plot showed the distribution of observed and expected p-values (association test statistics). The Q-Q-plot analysis for all traits, except stem height, showed a fairly mid-way deviation of the observed p values from those expected, towards the top left.

2.7. Identification of Candidate Genes

Gene annotation information according to Phytozome v13 (Manihot esculenta v7.1) (https://phytozome-next.jgi.doe.gov/) and National Centre for Biotechnology Information (NCBI) of cassava, Manihot esculentum genome (Manes) (https://blast.ncbi.nlm.nih.gov/) were used to determine the putative functionality of genes around the associated marker loci. Out of the 27 SNPs involved in MTA, 05 SNPs with Phenotype_Variance_Explained (PVE) ≥ 5% were selected to predict genes associated with the flowering traits. Consequently, a total of 37 candidate genes located within 50kb upstream or downstream of the five markers in the cassava genome were found to be associated with four of the five flowering traits; branching, branch type, branching levels and nodes at first branching (Table 4). SNPs S5_29309724 (on chromosome 5) and S15_11747301 (on chromosome 15) were associated with all the four traits and had 05 and 12 putative genes respectively. SNP S18_1832353 (on chromosome 18) was associated with three traits involving 10 genes, while the rest of the SNPs were associated with one trait each, but with varying numbers of putative genes. Detailed information about all significant associations is summarized in Table 5.

Marker R² values ranged from 0.017 to 0.177, with P-values varying from 1.15×10⁻⁹ to 1.45×10⁻⁷, while PVE values considered for gene identification ranged from 4.85 to 26.85% (Table 5). Averagely, 26.86% was the highest phenotypic variance, attributed to SNP S5_29309724 (associated with four flowering traits), while 4.85% was the lowest PVE, attributed to SNP S5_22566689 (associated with Nodes at 1st branch). Whereas 10 genes (06 on SNP S3_21330906 and 04 on S5_22566689) were unique to Nodes at 1st branch, the remaining 27 genes were associated with all the four flowering traits (branching, branch type, branching levels and Nodes at 1st branch). A total of 18 gene loci were each associated with four traits, for example, Manes.15G140500.1 (eukaryotic translation initiation factor 3 subunit I) and Manes.15G140900.1 (myb family transcription factor MOF1). Similarly, 09 candidate genes such as Manes.18G016700.1 (aldehyde oxidase GLOX) and Manes.05G186700.1 (protein DETOXIFICATION 48) were associated with three traits.

On the basis of the 53 associated markers identified using GWAS, candidate genes with known functions and orthologous with Arabidopsis or other plant species were identified using five SNPs which had PVE ≥ 5% (Table 6). Among the genes found to be associated with these SNPs, we focussed only on those which were reported to be functionally controlling flowering or flowering responses in various plant species. As such, 08 genes were identified as candidates and all were associated with branching, branch type and branching level traits. For example; Manes.18G016700.1 (aldehyde oxidase GLOX) as well as Manes.18G016725.1 and Manes.18G016800.2 (lysine-specific histone demethylase 1 homolog 3), all associated with SNP S18_1832353 on chromosome 18, are involved in anther and pollen development [58] and promotion of floral transition [59] respectively in Arabidopsis thaliana. Manes.05G186700.1 (protein DETOXIFICATION 48) associated with SNP S5_29309724 on chromosome 5 is involved in specifying the lateral organ initiation rate in Arabidopsis thaliana [60]. Associated with SNP S15_11747301 on chromosome 15 are Manes.15G140500.1 (eukaryotic translation initiation factor 3 subunit I) which was reported to regulate negatively translation during flower development [61], Manes.15G140900.1 (myb family transcription factor MOF1) which plays a role in the regulation of organ identity and spikelet meristem determinacy in rice [62] and Manes.15G140300.1 (non-specific lipid-transfer protein 4.1) involved in seed and ovule maturation and development in Hordeum vulgare (Barley) and Arabidopsis thaliana [63]. Manes.15G140600.1 (short-chain dehydrogenase reductase ATA1), also located on chromosome 15 was reported to play a role in tapetum development in Arabidopsis thaliana [64]. The new genes identified in the present study are promising candidates for follow-up studies on the validation of genes controlling flowering and related traits in cassava.

3. Discussion

In this study, GWAS was used to detect genomic regions associated with flowering traits in a diverse population of 293 cassava accessions. This enabled generation of knowledge to enhance understanding of the genetics of cassava flowering traits, i.e. branching, branch type, branching levels and number of nodes of main stem at first branching.

3.1. Variation in Phenotypic Traits Related to Flowering

Ability to flower as well as early and synchronized flowering are important target traits for selection of progenitors in cassava hybridisation programmes. Branching, branch type, branching levels, number of nodes and stem height at first branch are major indicators of flowering ability in cassava. Commonly, it is through flowering and pollination that crucial genetic variation is generated, which is fundamental for success of plant breeding programs [67]. Understanding phenological variabilities exhibited in cassava flowering is essential for accurate selection of parents and planning crosses in hybridization programs. However, these different phenologies are difficult to understand due to the complex genetics of the crop.

Number of nodes followed by stem height exhibited highest variability (Table 1). There is greater potential for selection of progenitors with capacity to flower basing on these attributes compared to others. High and exploitable variability for flower and fruit numbers as well as time to flower in cassava was also reported in a previous study [11]. Conversely, the lowest variation recorded for branching indicates low exploitable variability compared to other traits.

Heritability estimates give insights into the extent of genetic control of a particular trait [67,68,69]. High heritability indicates less environmental influence [67,70]. Broad-sense heritability varied across traits, i.e. pistillates had the highest heritability (H₂ = 0.83) followed by stem height at first branch (H₂ = 0.60). There are not many reports on broad sense heritability estimates for flowering traits in cassava, however, heritability estimates of over 0.80 were reported for plant height at first branch, number of branches, branching levels as well as female flowers [11,71,72,73]. SNP-based heritability (h₂) estimates were relatively higher for stem height at first branch (h₂ = 0.35), followed by branching (h₂ = 0.25) and branching levels (h₂ = 0.19). These estimates could be attributed to linkage disequilibrium (LD) unevenness or heterogeneity among marker regions, due to overestimation or underestimation, as postulated by [74].

3.2. Correlation between the Flowering Traits

Correlation coefficient analysis is important for measuring the degree and direction of relationships between various traits. The positive and significant correlation among the flowering traits in this study, such as branch and branch type (r² = 0.95) or flowers (r² = 0.94); branching levels and pistillate flowers (r² = 0.76), branching levels and height of first branch (r² = 0.88), etc. was consistent with the previous findings by [75]. Thus, this reflects that selection for one trait would directly affect the expression of the other trait, hence facilitating the selection of progenitors to include in a breeding program. Relatedly, a correlation between branching levels and height of first branch was reported in a previous study on induction of cassava flowering [76]. Overall, all these traits are dependable indicators of flowering response, upon which selection of progenitors can be based.

3.3. Population Structure and Kinship

Population stratification and the degree of kinship of the samples are critical factors in genetic studies including GWAS [77,78]. Whereas population structures show presence of distinct genetic subgroups within a population [79], kinships provide valuable insights into the degree of genetic relatedness between the individuals within a population [77]. PCA and kinship analysis revealed clustering of germplasm into subgroups, three major subgroups according to the kinship heatmap (Figure 4). The clusters (within the PCA and kinship plots) show genetic relatedness within members of the subgroup, suggesting a common ancestry within the subgroup. Clusters along the diagonal line in the kinship plot show unrelatedness among groups, while pairs of individuals along the diagonal line show individuals' self-kinship. Shorter connecting lines between individuals or clusters indicate closer relationships while longer connecting lines indicate more distant relationships.

The Q-Q-plot analysis further revealed subgroups in the study population. Plots for all traits, except stem height, showed a fairly mid-way deviation of the observed p values from those expected, towards the top left (Figure 4). This indicated that, firstly, the SNPs associated with these traits were significant, secondly, the existence of genetic variations between the subpopulations. Thus, suggesting that genetic variants have different allele frequencies among subgroups of the population investigated. For stem height at first branching, the little or non-segregation between the observed and expected p values suggests that the sub populations are genetically homogeneous, indicating a similarity in the distribution of allele frequencies between the subpopulations.

3.4. Linkage Disequilibrium

Linkage disequilibrium (LD), a non-random association or correlation between alleles at different genetic loci within a population plays an integral role in genetic studies, such as genome-wide association studies (GWAS) and genomic selection [80,81]. When LD exists between a marker (e.g. a SNP) and a phenotypic trait, it implies that the marker and the causative variant of that trait are physically close to each other on the same chromosome [81].

There is limited knowledge about genetic architecture and potential markers specifically associated with flowering traits in cassava due to little research that has been conducted on the LD patterns and LD decay for these traits. In this study, the average LD among the SNP marker pairs of 293 cassava accessions was 0.32 at a distance of ~2850kb which slowly declined to 0.16 at distance ~4828kb where it became more or less stable (Figure 5). A study by [82] on evaluation of genetic diversity among 1,580 cassava accessions presented a lower LD value average of 0.014, though the initial LD decline was more rapid before it stabilized at 15 to 20 kb. Also, analysis of a panel of 876 cassava accessions belonging to the Cassava Germplasm Bank of Embrapa, a LD of r² < 0.1 which declined to between 0.3 and 2.0 Mb was reported [83]. While a whole-genome LD decay peaked at r² of 0.349 and dropped to an r² of 0.212 at a distance of 10 kb in a study on genetic architecture of defensive, agro-morphological and quality-related traits involving a panel of 5130 cassava clones [84]. These observations show that the LD estimates and decay rates in the cassava genome are low and inconsistent.

The low and slow decline in LD observed in this study and the observed inconsistencies among the reported LD values of the different studies could be attributed to the genetic complexity (high heterozygosity) of the crop as well as the vegetative mode by which the crop is predominantly propagated. Additionally, the artificial selection pressures imposed by preferential selections based on specific traits, like disease resistance and high yield, tends to favour specific alleles linked to the desired traits, leading to extended LD around these beneficial variants. In potato (Solanum tuberosum L.), a clonally propagated crop, but outcrossing and highly heterozygous like cassava, LD decay was reported to be relatively fast in the short range (r² = 0.208 at 1 kb) but slowed afterward (r² = 0.137 at 70 kb) [85]. Perennial or clonally propagated plants including cassava, generally have long breeding cycles and, thus, show a limited number of recombination cycles. Hence, their LD decays are relatively slow [86] in spite of the outcrossing nature of these crops.

3.5. Marker–Trait Association Mapping

Association or linkage disequilibrium (LD) mapping, also known as association analysis or genome-wide association study (GWAS), is an important tool in plant breeding for understanding the genetic basis of complex traits. It involves surveying genetic variations in the whole genome to find signals of association between genetic markers, such as single-nucleotide polymorphisms (SNPs), and various phenotypic traits of interest within a diverse population [80]. It is also used in identifying candidate genes and genomic regions which are presumed to be associated with traits of interest [80,87].

Using the MLM model, we conducted a GWAS of five flowering traits for 293 cassava accessions based on 22,103 (DArTseq-derived) filtered SNP markers and detected 53 MTA on nine of the 18 chromosomes, with chromosome 18 possessing up to 20 of the 53 MTA. All traits, except stem height at first branching, had MTA on different chromosomes, suggesting pleiotropic effects. For example, branching was associated with 15 markers on chromosomes 5, 8, 11, 15, 16 and 18; while branching levels had 14 MTA on chromosomes 5, 8, 9, 11, 14, 15 and 18 (Table 3). This suggests that these traits are polygenic in nature. Meanwhile, some chromosomes had more than one MTA for some traits, for example, chromosome 18 was associated with 08, 04, 06 and 03 MTA for branching, branch type, branching levels and number of nodes at first branching respectively. This suggests that some traits are regulated by multiple genes located in close proximity to each other on the chromosome, and that, the genes may be interacting to contribute to the expression of the trait. This can be advantageous because it would enable narrowing down the search for traits of interest. A study by [83] also revealed multiple MTA on chromosome 18 for cassava starch pasting properties, indicating high LD of markers on this chromosome. Multiple MTA have also been detected for flowering traits in Indian mustard (Brassica juncea) [55] and grain yield and related traits in durum wheat (Triticum turgidum) [88]. Other traits had only one MTA on a chromosome, such as S11_32333764 on chromosome 11 (associated with branching); S15_11747301 on chromosome 15 (associated with branching, branch type, branching levels and number of nodes). This suggests that the marker may be associated with a potential gene or genomic region primarily responsible for influencing the given traits on that particular chromosome. This may be helpful particularly in genome editing.

3.6. Putative Candidate Genes Linked to Marker Loci for Flowering Traits

Using MLM in GWAS, different genomic regions (putative candidate genes) associated with flowering traits were predicted. All the predicted genes were associated with branching characteristics. These genes are novel and a new discovery, as they have not been previously reported.

Branching of shoots involves axillary bud initiation in which the shoot apical meristem is converted to axillary meristem, expansion as well as sustained branch growth, which is regulated by an interplay of genetic programs with the environmental signals, phytohormones and sugars [13,89]. Genes modulating reproductive branching have been identified in many plant species and several are broadly conserved in between species, though detailed phenotypic and developmental effects are species specific [89]. For example, LATERAL SUPPRESSOR in tomato, its Arabidopsis ortholog LAS [90], and the rice ortholog MONOCULM1 (MOC1) [91], the TEOSINTE-BRANCHED1 (TB1)/BRANCHED1(BRC1) (TB1/BRC1) gene in Arabidopsis, peas, tomato [89], BRC1 was also predicted in modulating branching in soybean [92]. In cassava, branching of the main stem is associated with transition from vegetative growth to reproductive development. This indicates on-set of flowering with flower-bud formation being preceded by a highly genetically planned apical branching [8].

In this study, three branching traits, i.e. branching, branch type and branching levels had commonality in 27 of the 37 predicted genes (Table 5). Of these, eight were reported to be functionally annotated for processes directly or indirectly leading to flowering responses in different plant species (Table 6). For example, Manes.18G016725.1 and Manes.18G016800.2 (lysine-specific histone demethylase 1 homolog 3), associated with SNP S18_1832353 on chromosome 18 was reported to influence the transition from vegetative growth to the flowering phase in Arabidopsis thaliana [59]. It reduces levels of histone methylation in chromatin of key flowering genes such as the floral repressor FLOWERING LOCUS C(FLC) which must be preferentially expressed in shoot and root apical regions during a plant’s vegetative development. This enables achievement of precise levels of FLC required for onset of a flowering response.

Manes.18G016700.1 (aldehyde oxidase GLOX, Germin-like oxalate oxidases) in close proximity with SNP S18_1832353 on chromosome 18. Though its role in plants is not yet fully understood, it is believed to be involved in anther development during which it plays a role in tapetum and pollen development in Arabidopsis thaliana [58]. It is also involved in catalyzation of the oxidation of aldehydes to the corresponding carboxylate by coupling the reaction to the reduction of dioxygen to hydrogen peroxide [65].

Manes.05G186700.1 (protein DETOXIFICATION 48) associated with SNP S5_29309724 on chromosome 5, a member of multidrug and toxic compound extrusion (MATE) gene family, is involved in specifying the lateral organ initiation rate in Arabidopsis thaliana [60]. It has been implicated in axillary bud formation in the shoot apical meristems resulting in accelerated organ formation including leaf formation and early flowering. It is also involved in detoxification or extrusion of toxins and regulation of plant development including fruit development [93].

Manes.15G140500.1 (eukaryotic translation initiation factor 3 subunit I) associated with SNP S15_11747301 on chromosome 15 is a member of eukaryotic translation initiation factor 3 (eIF3) complex [61]. It was reported to regulate negatively translation during flower development in Arabidopsis thaliana. It is responsible for initiating protein synthesis during development and germination of pollen grains and embryogenesis. Its expression was detected in plant cells and tissues differentiating into organs including inflorescences, implying its involvement in influencing transition from vegetative to reproductive growth [61]. Manes.15G140900.1 (myb family transcription factor MOF1), also associated with SNP S15_11747301, plays a role in the regulation of organ identity and spikelet meristem determinacy in rice [62]. It causes delayed transition from the spikelet to the floral meristem, which may result in defective or sterile spikelets.

Manes.15G140300.1 (non-specific lipid-transfer protein 4.1) associated with S15_11747301 on chromosome 15 plays a role in seed and ovule maturation and development, probably by regulating the fatty acids homeostasis during suberin and sporopollenin biosynthesis or deposition [94]. And Manes.15G140600.1 (short-chain dehydrogenase reductase ATA1) also near SNP S15_11747301 was reported to play a role in tapetum development in Silene latifolia and Arabidopsis thaliana [64].

4. Conclusions

The present study was undertaken to determine the genomic regions associated with flowering traits in cassava under field conditions. The interest in unravelling the genetic basis of cassava flowering is based on the premise that many breeding programmes find it still difficult to synchronize flowering time among different cassava genotypes. To date, this still remains as a significant impediment to cassava breeding. Thus, increasing our understanding of the genetics of the flowering process and related traits in cassava would greatly assist the ability to accelerate and synchronize flowering among different cassava varieties through marker-assisted breeding and/or gene editing. To this end, this study identified 53 significant MTA distributed mainly on chromosome 18, as being associated with flowering traits. Five SNPs, i.e. S3_21330906; S5_22566689; S5_29309724; S15_11747301 and S18_1832353 on chromosomes 3, 5, 15 and 18 respectively, were found to be significantly associated with number of nodes and branching traits. Eight genes, i.e. Manes.18G016700.1 (aldehyde oxidase GLOX); Manes.18G016725.1; Manes.18G016800.2 (lysine-specific histone demethylase 1 homolog 3); Manes.05G186700.1 (protein DETOXIFICATION 48); Manes.15G140500.1 (eukaryotic translation initiation factor 3 subunit I); Manes.15G140900.1 (myb family transcription factor MOF1); Manes.15G140300.1 (non-specific lipid-transfer protein 4.1) and Manes.15G140600.1 (short-chain dehydrogenase reductase ATA1) were predicted to be involved in controlling branching, branch type and branching level traits. However, further studies are required to validate the significant markers and to verify the identified candidate genes using a larger population.

5. Materials and Methods

5.1. Plant Material and Field Trials

A diverse panel, consisting of 446 cassava accessions with unknown and/or diverse flowering behaviour were used in this study. These accessions were derived from a cycle two (C2) population that was developed at NaCRRI (National Crops Resources Research Institute). Briefly, the C2 population resulted from successive cycles of selection and hybridization [95,96]. These accessions were established in trials laid out in augmented design with three checks of Mkumba, TME 204 and NAROCASS 1. The clones Mkumba and NAROCASS 1 are known for early forking and/or profuse flowering, while TME 204 is a late or non-forking. Each plot had one row of 10 plants spaced at 1m between plants and rows. The trials were established at two locations, Namulonge in Central Uganda (0.52166458^oN, 32.608997564^oE) at an altitude elevation of 1150m above sea level (asl) and at National Semi Arid Resources Research Institute (NaSARRI), Serere in Eastern Uganda (1.4994^oN, 33.5490^oE) at 1140m asl. The trial was conduvcted for two seasons, 2019/2020 and 2020/2021. All trials were managed according to the Institute’s standard agronomic practices for cassava.

5.2. Flowering Traits Evaluated

Plant-based counts of flowers (staminate and pistillate) were commenced at about three months after planting (MAP), and this was continued every seven days. However, this collection was done only for one season due to interruption by COVID19 lockdown in 2019 and 2020. Meanwhile, at 12 MAP data on the five traits related to flowering were collected on a per plant basis. The traits included: i) first branching of the main stem (Branch) which was scored as “1” for branched or “0” for non-branched, ii) branching type or habit (Branch1_No) as number of branches at first branching event or tier1, iii) number of branching levels of the main stem (Branch_Levels), iv) number of nodes (Branch1_Nodes), counted with tally counters up to first branching, and v) stem height to first branching (Branch1_Ht) from ground to the first branch of the main stem, measured with a meter rule or tape measure. These traits were prioritized basing on reports which inferred that stem branching was an indication and precondition for flowering in cassava [7,8,10]. Data on these traits were collected for two consecutive seasons from 2019–2021 in two locations. All data were electronically collected on tablets using a Field Book application [97].

5.3. Genomic DNA Extraction and Genotyping

Extraction of genomic DNA and genotyping was carried as described in [46]. Briefly, two young top leaves were collected from the study genotypes, folded, punched using a 5mm hand puncher and placed in 96-well plates before being shipped to Intertek, Australia for genomic DNA extraction and purification. Extraction was done using a DNA extraction protocol routinely used at Intertek-AgriTech (http://www.intertek.com/agriculture/agritech/) and the extracted DNA samples were purified and quality as well as quantity evaluated before being genotyped. Genotyping was done using the Diversity Array Technology Pty Ltd. (DArTseq) genotyping platform (https://www.diversityarrays.com/technology-and-resources/dartreseq/).

5.4. SNP Calling and Annotation

The sequences of the genomic representations were aligned to cassava reference genome v6.1 resulting in the selection of 28,434 raw SNP markers. The genotype data were then subjected to stringent filtration measures to remove genotypes with >10% and SNPs with >5% missing data or with minor allele frequency of less than 5%. Thus, a total of 24,040 high quality SNP markers were used for GWAS. The SNPs were annotated on the basis of the cassava reference genome v6.1. The markers were converted to the dosage format of 1, 0, 2, which represented alternative allele homozygotes, heterozygotes, and reference allele homozygotes, respectively.

5.5. Statistical Analyses

All phenotypic data were subjected to various statistical analyses using R Statistical Software v4.2.2 [98]. Accessions that did not have a complete set of data points across seasons and location were excluded from the statistical analyses, thus leaving 293 out of the 446 initial accessions. BLUPs for each trait were computed and used for statistical and GWAS analyses. The phenotypic datasets; number of branches, branching levels and nodes were analysed, as recommended for count data, using negative binomial distribution, [99,100] while analysis of variance (ANOVA) was performed for the parameter of stem height. Generalized Linear Mixed Models and lme4 packages for R statistical software were used. Broad-sense heritability (H²) per trait was computed using accession and location variance components extracted from the ANOVA models as:

$$H^2=\frac{V_G}{(V_G+{\frac{{\mathrm{V}}_{G}xE}{L}+\frac{{\mathrm{V}}_R}{LxS}}_{ }^{ })}$$

where ${\mathrm{V}}_G$ was the accession variance; ${\mathrm{V}}_G x E$ was the variance attributed to accession by location interaction; ${\mathrm{V}}_R$ was the model residual variance; $L$ was number of locations; and $S$ was number of seasons.

SNP-based heritability (h²) was computed as:

$$h^2=\frac{{\mathsf{\sigma }}_{\mathrm{a}}^2}{({\mathsf{\sigma }}_a^2+{\mathsf{\sigma }}_{\mathrm{e}}^2)},$$

where ${\mathsf{\sigma }}_{\mathrm{a}}^2$ was the additive genetic variance and ${\mathsf{\sigma }}_{\mathrm{e}}^2$ was the residual variance.

Correlation analysis for flowering (1. first branching of the main stem; 2. branching type; 3. number of branching levels; 4. number of nodes and 5. stem height at the first branching as well as number of staminate and pistillate flowers) traits was conducted based on phenotypic BLUPs. The analyses were performed using the cor function, and visualization of the correlation matrices was done using the ‘corrplot’ package in R [98].

5.6. Population Structure and Kinship

Population structure analysis was performed to determine the genetic relatedness between individuals using 24,040 polymorphic SNP markers distributed over the whole genome. Using the Genome Association and Integrated Prediction Tool (GAPIT) 3.0 in R, we computed a genomic relationship matrix (K) (kinship) from the SNP data. Principal component analysis (PCA) on all markers was done on the genetic relationship matrix, using the prcomp function in R. The first three principal components (PC1, PC2 and PC3) values were plotted to visualize population structure.

5.7. Linkage Disequilibrium Estimation

Linkage disequilibrium (LD) (in terms of r²) analysis was performed for the whole genome as well as individually for each of the eighteen cassava chromosomes using package LDcorSV in R. To estimate LD decay, mean r² value of each chromosome was computed and plotted against the chromosome’s physical distance to obtain a non-linear regression curve, and the LD decay distance was estimated as the physical distance between SNPs where average r² reduced to half of the maximum LD value.

5.8. Marker-Trait Association Mapping

Search for marker-trait associations between genotypic and phenotypic data was performed for 22,103 filtered SNPs and the BLUP values obtained for the five flowering traits among the 293 filtered accessions. This was done using the Genomic Association and Prediction Integrated Tool (GAPIT) in R [98] and Mixed Linear Model (MLM) which was best at controlling type-I errors (false positives) for all of the traits evaluated. The first five PCs and marker-estimated kinship matrix in each case, were used. The genome-wide significance thresholds for each association study were assigned using Bonferroni correction set as 0.1. Significance of associations between traits and markers was assessed on the basis of p-values corrected for multiple testing. Manhattan plots as well as quantile-quantile (Q-Q) plots for association mapping were used to examine the validity of the SNPs and were visualized using the qqman package in R [98].

5.9. Identification of Putative Candidate Genes

To identify putative candidate genes residing at the close vicinity of high confidence SNPs, the associated SNPs with PVE ≥ 5% were mapped to the reference cassava genome, Manihot esculenta, v7.1 using the Phytozome tool v13 (https://phytozome-next.jgi.doe.gov/jbrowse/index.html?data=genomes/Mesculenta_v7_1). Transcripts present within 50Kb regions from both sides of associated SNPs were fetched along with their description. Identified genes were characterized for gene ontology including molecular and biological functions using UniProtKB tool (https://www.uniprot.org/). Additionally, gene and protein functions were searched using Alliance of Genomes Resources [101].