Introduction

Coffee is one of the most economically important beverage crops worldwide and ranks among the most highly traded agricultural commodities (Murthy and Madhava Naidu 2012). Global coffee production, however, is increasingly threatened by climate change. Elevated temperatures, irregular precipitation patterns, and extreme climatic events have reduced suitable cultivation areas and intensified vulnerability to pests and pathogens (Bunn et al. 2015; Velásquez et al. 2018). A recent report from the Food and Agriculture Organization indicates a significant price increase in 2024, largely driven by climate-associated supply instability (FAO 2025). These challenges underscore the urgent need to enhance genetic resilience to both biotic and abiotic forms of stress in coffee.

The genus Coffea (Rubiaceae), belonging to the Euasterid I clade, comprises approximately 90–124 species (Davis 2011). Among these, Coffea arabica and C. canephora account for nearly all global production of coffee (Clarindo et al. 2013; ICO 2018). Most Coffea species are diploid (2n = 2x = 22) and self-incompatible, whereas C. arabica is a natural allotetraploid (2n = 4x = 44) derived from hybridization between C. canephora and C. eugenioides approximately 350–610 thousand years ago (Lashermes et al. 1999; Salojärvi et al. 2024). Due to its recent origin and predominant self-fertilizing nature, C. arabica exhibits limited genetic diversity relative to its diploid progenitors, possibly contributing to its heightened susceptibility to environmental stresses and pathogens (Spinoso-Castillo et al. 2021).

Despite its economic importance, genetic improvements through molecular breeding and genome-editing applications in coffee remain underdeveloped relative to those in other major economic crops. The recent availability of chromosome-level genome assemblies for C. canephora (Denoeud et al. 2014), C. eugenioides (Salojärvi et al. 2024), and two C. arabica cultivars (Geisha and Typica) (Lyu et al. 2025; Medrano et al. 2025) provides an unprecedented opportunity for genome-wide gene family analyses. Several important gene families related to hormone signaling and stress responses have been characterized (Acevedo-Garcia et al. 2014; Zanin et al. 2022; Naeem et al. 2024; Park et al. 2025). However, disease-resistance and susceptibility-associated gene families remain poorly explored in coffee.

Powdery mildew (PM) caused by fungi of the order Erysiphales (Ascomycota) (Glawe 2008; Hacquard 2014) affects a wide range of angiosperms. While PM is not currently as devastating in coffee as coffee leaf rust (Hemileia vastatrix) or coffee berry disease (Colletotrichum kahawae), infections can reduce photosynthetic efficiency and plant vigor, ultimately impacting yields and quality levels negatively, as reported in many other important crops (Scholes and Rolfe 1996; Glawe 2008; Fondevilla and Rubiales 2012). Natural PM resistance was first reported in barley (Hordeum vulgare), where loss-of-function mutations in the Mildew Resistance Locus O (MLO) gene confer broad-spectrum resistance to PM (Büschges et al. 1997; Kusch and Panstruga 2017). Later studies demonstrated that mlo-based resistance is recessively inherited and non-race-specific, providing long-lasting protection against diverse pathogen isolates (Pavan et al. 2010; Kusch and Panstruga 2017), unlike race-specific, dominant R-gene-mediated resistance (Moffett 2009). Further research showed that mlo-based resistance is conserved across distantly related dicot species, including tomato (Solanum lycopersicum) (Bai et al. 2008), pea (Pisum sativum) (Pavan et al. 2011), pepper (Capsicum annuum) (Zheng et al. 2013) and Arabidopsis (Consonni et al. 2006; Acevedo-Garcia et al. 2014). In Arabidopsis thaliana, three paralogs (AtMLO02, AtMLO06, and AtMLO12) act redundantly as susceptibility genes (Consonni et al. 2006), and the triple mutant mlo2mlo6mlo12 confers near-complete immunity against PM (Acevedo-Garcia et al. 2017). More recently, targeted mutagenesis of TaMLO genes in wheat and SlMLO01 in tomato using genome-editing technologies successfully generated broad-spectrum resistance (Wang et al. 2014; Nekrasov et al. 2017). Given the demonstrated utility of MLO loss-of-function alleles for durable PM resistance, systematic genome-wide identification and evolutionary analyses of the MLO gene family in Coffea species are highly warranted.

In this study, we undertook genome-wide identification and conducted a comparative analysis of the MLO gene family in two commercially important C. arabica cultivars, Geisha and Typica, and their parental progenitors, C. canephora and C. eugenioides, utilizing recently available chromosome-level genome assemblies (Supplementary Table S1). In total, we identified 90 MLO genes across the four Coffea genomes. By integrating comparative evolutionary and structural analyses with model species, we characterized the diversification and regulatory features of the coffee MLO gene family. Our findings provide evolutionary insights into the coffee MLO gene family and putative susceptibility MLO candidates that can be targeted for resistance breeding against PM in coffee.

Materials and Methods

Genome-wide identification of the MLO gene family

Genome assemblies, predicted protein sequences, and corresponding annotation files (GFF3) for C. arabica cv. Geisha were obtained from Zenodo (https://zenodo.org/records/10059814), and those for C. arabica cv. Typica were retrieved from Figshare (https://figshare.com/articles/dataset/A_chromosome-level_genome_assembly_of_Coffea_arabica_L_var_Kona_Typica/28425329). The genome data for C. canephora (GCA_900059795.1) and C. eugenioides (GCF_003713205.1) were downloaded from the NCBI Genome database (Supplementary Table S1). To ensure dataset consistency and to evaluate annotation quality across all species, assembly completeness was independently assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO). For each genome, assembly completeness was independently calculated using the eudicots_odb10 database for C. canephora and the brassicales_odb10 database for A. thaliana.

Coffee MLO genes were identified using TBLASTN (version 2.17.0+) (Altschul et al. 1997) searches with the 15 previously characterized A. thaliana MLO protein sequences (NP_192169.1, NP_172598.1, NP_566879.1, NP_001320916.1, NP_180923.1, NP_176350.1, NP_179335.3, NP_565416.1, NP_174980.3, NP_201398.1, NP_200187.1, NP_565902.1, NP_567697.1, NP_001319086.1, NP_973686.1) as queries. Searches were conducted against the genomes of C. arabica cv. Geisha, C. arabica cv. Typica, C. canephora, and C. eugenioides using an E-value threshold of 1e-5.

In parallel, MLO-related sequences were independently retrieved using InterProScan (version 5.59–91.0) (Jones et al. 2014) based on the presence of conserved domains annotated as InterPro entry IPR004326 and/or Pfam entry PF03094 (Paysan-Lafosse et al. 2025). The two datasets were subsequently merged, and redundant or partial hits were manually curated to ensure high-confidence gene models.

The amino acid sequences of the retrieved candidate genes were analyzed using the ExPASy-ProtParam (Gasteiger et al. 2003) tool to determine their physicochemical properties, including the Isoelectric point (pI) and Molecular weight (Mw). Subsequently, only gene models containing at least ten exons, possessing at least five transmembrane domains predicted using TMHMM online software (version 2.0) (Krogh et al. 2001), and predicted to localize to the plasma membrane according to both Plant-mPLoc (version 2.0) (Chou and Shen 2010) and DeepLoc (version 2.1) (Thumuluri et al. 2022) were retained as coffee MLO candidates. When multiple isoforms were identified for a single locus, the isoform with the longest transcript was selected for downstream analyses. The final set of genes was renamed according to the chromosomal position within each genome (Supplementary Table S2).

Phylogenetic tree construction

To classify MLO genes in coffee, MLO protein sequences from A. thaliana, O. sativa, and S. lycopersicum were aligned with those identified in C. arabica cv. Geisha, C. arabica cv. Typica, C. canephora, and C. eugenioides using MAFFT (version 7.520) with the –globalpair and –maxiterate 1000 parameters (Katoh and Standley 2013). Alignment sites containing more than 10% gaps were removed using trimAl (version 1.4.rev15) (Capella-Gutiérrez et al. 2009). A maximum-likelihood phylogenetic tree was reconstructed using IQ-TREE (version 2.3.6) with the -alrt 1000 and -B 1000 parameters (Minh et al. 2020), with the best-fit substitution model (Q.plant+I+R6) determined by ModelFinder (Kalyaanamoorthy et al. 2017). The MLO clades were assigned according to previously described phylogenetic classification and taxonomic distributions (Acevedo-Garcia et al. 2014). These clades were strongly supported by ultrafast bootstrap values exceeding 90%. The phylogenetic tree was visualized using iTOL (v7.3) (Letunic and Bork 2024).

To construct a species-specific phylogeny for the gene structure analysis, coffee MLO protein sequences were aligned using MAFFT (version 7.520) with the –globalpair and maxiterate 1000 parameters (Katoh and Standley 2013), and alignment positions with more than 10% gaps were removed with trimAl (version 1.4.rev15) (Capella-Gutiérrez et al. 2009). Maximum-likelihood phylogenetic trees were reconstructed using IQ-TREE (version 2.3.6) with the -alrt 1000 and -B 1000 parameters (Minh et al. 2020). The best-fit substitution models identified by ModelFinder (Kalyaanamoorthy et al. 2017) were Q.pfam+I+G4 for C. arabica cv. Geisha, C. canephora, and C. eugenioides, and Q.plant+I+G4 for C. arabica cv. Typica.

Structural comparisons of coffee MLO proteins in Clade V-ii

Protein structures of coffee MLO proteins were predicted using AlphaFold3 (Abramson et al. 2024) with default modeling parameters. For each protein, the highest-ranking model was selected based on predicted local distance difference test (pLDDT) confidence scores. Structural comparisons between the predicted coffee MLO models and the Arabidopsis AtMLO02 structures were performed using the TM-align algorithm implemented in the RCSB Protein Data Bank (Berman et al. 2000) pairwise structure alignment tool. Structural similarity was quantitatively assessed using the root-mean-square deviation (RMSD) and template modeling score (TM-score) metrics.

Gene structure and motif analysis of coffee MLO genes

Exon–intron structures were visualized based on GFF3 annotation files using custom Python scripts with pandas (McKinney 2010) and matplotlib (Hunter 2007). The conserved motif analysis was performed using MEME (v5.5.7) (Bailey et al. 2015), specifying up to 15 motifs with a maximum width of 50 amino acids. Predicted transmembrane domains were also identified during the motif analysis. An additional discriminative motif analysis was undertaken to identify motifs enriched in Clade V-ii proteins using MEME (v5.5.7) (Bailey et al. 2015) with the parameters -objfun de, -mod zoops, -nmotifs 40, -minw 2, and -maxw 9.

To investigate functionally relevant motifs further, sequences corresponding to the calmodulin-binding domain (CaMBD) and conserved peptide domains I and II were examined. Motif occurrences were scanned using FIMO with a significance threshold of 1 × 10^-4. The MEME output file was used as input, and the search was conducted against the combined protein sequence dataset.

Ka/Ks analysis of coffee MLO genes

To investigate the evolutionary dynamics of the MLO gene family in coffee, a Ka/Ks analysis was performed separately within each of the four Coffea genomes using all pairwise combinations of MLO paralogous gene pairs identified in each species. Ka/Ks ratios were calculated using KaKs_Calculator v3.0 (Zhang 2022) with the following parameters: genetic_code = 1, method = MYN, iteration = 5, and sleep_time = 10. Gene pairs with Ks values greater than 2.0 or P-values greater than or equal to 0.05 were excluded to improve the reliability of substitution rate estimates and to focus on relatively recent duplication events.

Chromosomal mapping and physical localization of MLO genes

Genomic coordinates and chromosome lengths were extracted from the GFF annotation file using custom Python scripts with pandas v2.3.3 (McKinney 2010), matplotlib v3.10.7 (Hunter 2007), numpy (Harris et al. 2020) and adjustText (version 1.3.0) (https://github.com/Phlya/adjustText). The physical positions of the identified MLO genes were mapped onto their respective chromosomes. Global gene density was calculated based on the number of genes within 1 Mb non-overlapping windows and visualized as a heatmap to provide genomic context. To ensure structural alignment and accurate comparison, the chromosomes of C. arabica cv. Typica were inverted to match the orientation of the C. arabica cv. Geisha assembly based on their respective genome annotations as references.

Synteny analysis of coffee MLO genes

To identify orthologous collinearity among MLO genes in coffee genomes, a synteny analysis was conducted using MCScan (Python version) implemented in the JCVI tool kit (Tang et al. 2024) with a cutoff parameter of CSCORE = 0.99.

In C. arabica cv. Geisha, two gene pairs (CagMLO17-CcMLO09 and CagMLO21-CcMLO12) lacking clear collinearity or exhibiting overlapping regions were additionally examined by means of a chromosomal position analysis and through macro-synteny comparisons. The genomic context surrounding these regions was inspected in detail, and the results were visualized using the pyCirclize package (v.1.10.1) (Shimoyama 2022).

Cis-elements analysis of coffee MLO genes

Promoter regions extending 1,500 bp upstream of the annotated gene start site for all coffee MLO genes were extracted and analyzed for cis-regulatory elements (CREs) using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al. 2002).

Results

Genome-wide identification and physiochemical properties of the MLO gene family in Coffea species

To identify MLO gene family members in coffee genomes, Arabidopsis AtMLO protein sequences were used as TBLASTN queries (Altschul et al. 1997) against the genomes of C. arabica cv. Geisha, C. arabica cv. Typica, C. canephora, and C. eugenioides (Supplementary Table S1). A total of 40 and 36 putative MLO-like genes were initially identified in the Geisha and Typica genomes, respectively. To ensure accurate identification, candidates with nine or fewer exons were excluded because plant MLO genes typically possess complex exon-intron structures with at least ten exons (Liu and Zhu 2008; Zheng et al. 2016; Davis et al. 2017; Hussain et al. 2025). The remaining 28 and 34 MLO candidates in Geisha and Typica, respectively, were filtered further based on transmembrane domains and predicted subcellular localization. After filtering, 26 and 32 membrane-localized MLO proteins with at least five predicted transmembrane domains were retained in Geisha and Typica, respectively. Applying the same criteria to the diploid progenitors identified 16 MLO genes in each of C. canephora and C. eugenioides (Supplementary Table S2). The higher gene number in C. arabica relative to its diploid progenitors is consistent with its allotetraploid origin and suggests substantial gene retention after hybridization and genome duplication. In total, 90 putative MLO genes were identified across the four Coffea genomes.

To validate the identification strategies, the same procedure was applied to A. thaliana and O. sativa recovering all previously reported MLO genes (Chen et al. 2006; Liu and Zhu 2008), respectively, supporting the robustness of the pipeline. A total of 14 candidate MLO genes were identified in S. lycopersicum cv. Micro-Tom (GCF_036512215.1) (Table 1). Likewise, the updated genome assembly for S. lycopersicum cv. Heinz 1706 (SL5.0) yielded 14 candidates, showing a slight discrepancy with the 16 candidates previously reported in the SGN release 2.4 database (Zheng et al. 2016). This discrepancy is likely due to genome annotation differences between genome assemblies or database versions.

The average MLO protein length and molecular weight were 534 aa/61.02 kDa in Typica and 540 aa/61.41 kDa in Geisha. In the diploid progenitors, averages were 525 aa/59.82 kDa in C. canephora and 540 aa/61.58 kDa in C. eugenioides (Supplementary Table S2). Predicted isoelectric point (pI) values varied considerably (Supplementary Table S2), although most proteins exhibited basic pI values (8.5–9.0), a common feature of integral membrane proteins (Wallin and von Heijne 1998). This alkalinity likely facilitates electrostatic interactions with negatively charged membrane phospholipids (Kim et al. 2002a). In Typica, 25 MLO proteins had basic pI (8.5–9.38) and seven were near-neutral (6.87–7.68). Similarly, Geisha was found to contain 19 basic (8.5–9.39) and six near-neutral (6.87–7.70) MLO proteins. In C. canephora, 13 proteins were basic, two near-neutral, and one acidic, whereas C. eugenioides contained 14 basic and two neutral proteins. Previously characterized MLO proteins exhibited divergent pI values, e.g., acidic AtMLO12 (pI 5.5) and basic HvMLO (pI 9.2), while maintaining conserved core functions (Supplementary Table S2) (Zhu et al. 2005). Accordingly, pI variation among coffee MLO proteins may reflect differences in subcellular localization or regulatory roles in plant immunity and development.

Phylogenetic classification and conservation of MLO gene architectures

MLO proteins in land plants are typically classified into seven phylogenetic clades (Clades I - VII) based on full-length amino acid sequences (Kusch et al. 2016). However, recent comparative studies have revealed lineage-specific expansions, contractions, and the emergence of additional clades in certain taxa, highlighting the dynamic evolutionary history of this gene family (Acevedo-Garcia et al. 2014; Kusch et al. 2016; Kusch and Panstruga 2017). To investigate the evolutionary relationships of coffee MLO proteins, a comprehensive phylogenetic analysis was conducted using representative sequences from dicots and monocot model species (Supplementary Table S2). MLO proteins from C. arabica cv. Geisha and Typica, C. canephora, and C. eugenioides were aligned with 15 A. thaliana MLO (AtMLO), 12 O. sativa MLO (OsMLO), 14 S. lycopersicum cv. Micro-Tom MLO (SlMLO), and 14 S. lycopersicum cv. Heinz 1706 MLO (SlHMLO) proteins (Supplementary Table S2). Physcomitrium patens MLO-like protein 1 was used as an outgroup. The maximum-likelihood tree resolved seven major clades (I-VII), consistent with previous classifications (Kusch et al. 2016). Coffee MLO genes from both allotetraploid cultivars and their diploid progenitors were distributed across all seven clades (Fig. 1), indicating that diversification of major MLO lineages predates speciation within the Coffea genus.

Clade V, which in dicots contain MLO proteins associated with PM susceptibility, was of particular interest. Six CatMLO proteins and their homologs from Geisha, C. canephora, and C. eugenioides clustered within this clade (Fig. 1), which was further divided into two subclades (V-i and V-ii). Clade V-i comprised ten coffee MLO proteins and four tomato homologs, including two SlMLO proteins (SlMLO07 and SlMLO14) and two SlHMLO proteins (SlHMLO07 and SlHMLO14). Notably, Clade V-ii included experimentally validated PM susceptibility genes from Arabidopsis and tomato, such as AtMLO02, AtMLO06, AtMLO12, SlMLO05, SlMLO06, and SlHMLO05 (Consonni et al. 2006; Acevedo-Garcia et al. 2014; Zheng et al. 2016). This suggests that the six coffee MLOs in this subclade are strong susceptibility gene candidates.

Fig. 1.

Phylogenetic analysis of the MLO gene family across coffee and representative plant species. A maximum-likelihood (ML) phylogenetic tree was constructed using IQ-TREE based on full-length protein sequences. Outer circular strips classify MLO proteins into distinct clades, while inner strips further resolve two subgroups within Clade V, with Clade V-ii highlighted in pink. Species and cultivar identities are indicated by prefixes and symbols: Arabidopsis thaliana (At), Solanum lycopersicum cv. Micro-Tom (Sl), Solanum lycopersicum cv. Heinz 1706 (SlH), Oryza sativa (Os), Coffea canephora (Cc), C. eugenioides (Ce), C. arabica cv. Geisha (Cag), and C. arabica cv. Typica (Cat). Nodes with ultrafast bootstrap (UFBoot) values 90 are marked with circles. The tree is rooted using Physcomitrium patens MLO-like protein 1 as an outgroup. The scale bar represents the number of amino acid substitutions per site.

Exon-intron organization plays a key role in gene evolution by enabling alternative splicing, exon shuffling, and regulatory diversification (de Souza et al. 1996; Keren et al. 2010). To assess structural variation, exon numbers and gene architectures were analyzed across the four Coffea genomes (Supplementary Fig. S1). In C. arabica var. Typica, CatMLO genes contain 14–16 exons, with CatMLO24 showing the highest exon number (16), while most genes contain 14 or 15 exons. In the diploid progenitors, exon numbers ranged from 13 to 16 in C. canephora and 14 to 16 in C. eugenioides, indicating strong conservation of exon-intron structures between Typica and its parental genomes. This suggests limited post-polyploid structural reorganization of MLO loci in Typica.

In contrast, greater structural variability was observed in C. arabica cv. Geisha, where exon numbers ranged from 11 to 16 (Supplementary Table S2). Six CagMLO genes exhibited exon-intron structures distinct from their Typica orthologs. Notably, CagMLO04 and CagMLO23 contained 11 exons, whereas their Typica orthologs (CatMLO31 and CatMLO26) contained 14 exons. These Typica genes shared identical exon structures with their parental orthologs (CeMLO02 and CcMLO14), indicating that structural divergence occurred specifically in the Geisha lineage.

A phylogenetic analysis further showed that exon-intron organization was not strictly conserved within clades. Although genes within the same clade shared high sequence similarity, their exon numbers and intron arrangements varied (Supplementary Fig. S1). This finding suggests that structural diversification occurred independently of clade-level relationships. Overall, MLO gene architectures are largely conserved in coffee; however, cultivar-specific structural variations, particularly in Geisha, may contribute to regulatory divergence and functional specialization.

Chromosomal distribution, synteny, and evolutionary constraints of the MLO gene family

To investigate the genome organization and evolutionary conservation of MLO genes in C. arabica, the chromosomal distribution and gene order were analyzed first in its diploid progenitors. In both progenitors’ genomes, 16 MLO genes were distributed across seven chromosomes, with no loci detected on chromosomes 4, 8, 9, or 11 (Supplementary Fig. S2). Specifically, three genes were located on chromosome 1, five on chromosome 2, one each on chromosome 3 and 5, and two each on chromosome 6, 7, and 10 (Supplementary Fig. S2). Chromosomal positions, gene order, and orientations were highly conserved between the two species, indicating strong structural collinearity and limited rearrangement.

To confirm conservation in allotetraploid C. arabica, physical mapping was performed in the Geisha and Typica cultivars and compared with the diploid progenitors. In Typica, 32 MLO genes were mapped to chromosomes corresponding to those of the progenitor species, largely maintaining a conserved gene order and orientation (Supplementary Fig. S2). In contrast, only 26 MLO genes in Geisha were mapped to homologous regions. Six loci were absent in Geisha relative to Typica and the progenitor genomes: one from C. canephora chromosome 2C (Chr2C), two from C. eugenioides chromosome 2 (Chr2E), two from chromosome 6E (Chr6E), and one from chromosome 7E (Chr7E).

The previous phylogenetic analysis classified coffee MLO genes into seven clades (Fig. 1). To resolve orthologous relationships and positional conservation further, a genome-wide synteny analysis was conducted across the four genomes (Fig. 2). In Typica, all 32 CatMLO genes formed 16 syntenic pairs with clear one-to-one correspondence to both CcMLO and CeMLO loci (Supplementary Tables S3 and Supplementary Tables S4), confirming extensive retention of orthologous gene pairs in Typica. In contrast, a synteny analysis of the Geisha genome identified only 27 syntenic pairs between CagMLO and the diploid progenitors (Fig. 2). Although these retain one-to-one correspondence with CcMLO and CeMLO loci, orthologs corresponding to CcMLO08, CeMLO07, CeMLO08, CeMLO11, CeMLO12, and CeMLO14 were not detected (Supplementary Tables S3 and Supplementary Tables S4). These missing orthologous relationships indicate that six MLO genes (CatMLO07, CatMLO13, CatMLO16, CatMLO23, CatMLO24, and CatMLO27) were not detected in the current Geisha genome assembly. This pattern may reflect gene loss or pseudogenization during C. arabica evolution following hybridization and genome stabilization, although genome assembly or annotation differences cannot be ruled out.

Fig. 2.

Synteny analysis of MLO genes in Coffea species. Genome-wide synteny relationships of MLO genes among C. canephora, C. eugenioides, and C. arabica are shown. The outer tracks represent chromosomes of C. arabica cv. Geisha (blue), C. arabica cv. Typica (green), C. canephora (purple), and C. eugenioides (orange). Colored lines in the center indicate collinear MLO gene pairs, with yellow lines representing those derived from the C. canephora subgenome and sky-blue lines representing those derived from the C. eugenioides subgenome.

The synteny analysis further supported phylogenetic classification, as genes within the same clade exhibited conserved orthologous relationships. For example, Clade V-ii members showed clear correspondence with parental loci: CagMLO11 and CatMLO11 aligned with CcMLO04, whereas CagMLO14 and CatMLO12 corresponded to CeMLO03. Overall, the MLO gene family was largely conserved between diploid progenitors and allotetraploid C. arabica. The observed differences between Geisha and Typica may reflect cultivar-specific structural variations, or differences in genome assembly or annotation.

To evaluate the evolutionary pressure on duplicated MLO genes, nonsynonymous (Ka) and synonymous (Ks) substitution rates were calculated and Ka/Ks ratios were determined. In Typica, Ka/Ks values ranged from 0.124 to 0.601 (mean 0.258), while in Geisha, they ranged from 0.106 to 0.448 (mean 0.223). In the diploid progenitors, Ka/Ks values ranged from 0.130 to 0.681 in C. eugenioides and from 0.124 to 0.57 in C. canephora, with mean values of 0.303 and 0.261, respectively (Supplementary Table S5). All predicted gene pairs exhibited Ka/Ks ratios below 1, confirming that purifying selection is the predominant evolutionary force acting on this gene family. The lowest Ka/Ks ratios were observed for CatMLO27/CatMLO30 (0.124) and CagMLO07/CagMLO10 (0.106), reflecting strong evolutionary constraints and high sequence conservation. Collectively, these results indicate that the coffee MLO gene family is largely conserved under strong purifying selection.

Potential susceptibility candidates via functional motifs and 3D structural homology

MLO proteins are integral membrane proteins involved in diverse processes, including PM susceptibility, root development, pollen tube reception, and light signaling (Chen et al. 2009; Kessler et al. 2010; Nguyen et al. 2016; Jones et al. 2017). Their functional diversity is thought to be associated with clade-specific conserved motifs (Kusch et al. 2016). A motif analysis across all 90 coffee MLO proteins identified 15 conserved core motifs (Supplementary Table S6). Although the motif number and arrangement varied globally, they were largely conserved within specific phylogenetic clades (Supplementary Fig. S1). Generally, while TM numbers vary among MLO proteins (Supplementary Table S2), they exhibit transmembrane architectures and a conserved C-terminal CaMBD, which mediates Ca²⁺-dependent interaction with calmodulin (CaM) (Kim et al. 2002b).

All Clade V-ii members retained the fifteen core motifs, primarily localized to cytosolic loops or adjacent regions of transmembrane domains, consistent with roles in regulatory complex assembly and protein-protein interaction (Kusch and Panstruga 2017). Furthermore, Clade V-ii members consistently retained the canonical 7-transmembrane structure (Supplementary Fig. S3). The IC2 domain, located between transmembrane 2 and transmembrane 3, was also highly conserved, consistent with its roles in structural stability and biological activity (Reinstädler et al. 2010). Notably, two conserved peptide motifs, designated peptide domains I and II, were identified downstream of the CaMBD in these proteins but were absent in other clades (Fig. 3A; Supplementary Tables S7, Supplementary Tables S8, and Supplementary Tables S9). Peptide domain I contained conserved serine and threonine residues implicated in PM regulation (Panstruga 2005a), while peptide domain II included the conserved DFSF motif, characteristic of susceptibility-associated MLO proteins (Panstruga 2005b; Kusch and Panstruga 2017).

Fig. 3.

Conserved motifs and C-terminal divergence of Clade V-ii MLO proteins. (A) Schematic representation of the protein domain architecture, showing the positions of the calmodulin-binding domain (CaMBD) and conserved peptide domains I and II, along with their sequence logos. (B) Multiple sequence alignment of the C-terminal region (residues > 470). Clade V-ii sequences from Coffea spp. and A. thaliana are compared with representative A. thaliana MLO proteins from Clades I, II, III, and VI. A red horizontal line separates Clade V-ii from other clades. Black boxes indicate conserved regions corresponding to the CaMBD and peptide domains I and II, consistent with the domain representations shown in panel (A).

To identify motifs linked to susceptibility specifically, multiple sequence alignments were performed on Clade V-ii MLO proteins. Focused analysis of these Clade V-ii proteins identified three highly conserved motifs: AKKXXKED within the CaMBD, RPXTPLHG within peptide domain I, and DFSF within peptide domain II (Fig. 3B). Importantly, these motifs were not only consistently present in coffee Clade V-ii proteins but were also highly conserved in known susceptibility-associated MLO homologs from other plants, such as Arabidopsis and tomato, highlighting their shared evolutionary role in disease susceptibility.

Given the strong sequence conservation at the primary sequence level, structural modeling was conducted to assess three-dimensional structural conservation among susceptibility-associated MLO proteins. Predicted 3D structures of the six Clade V-ii MLO proteins were generated using AlphaFold3 (Abramson et al. 2024). All six proteins were modeled across their full-length sequences, with pTM scores ranging from 0.65 to 0.67. Structural superposition with AtMLO02 showed structural alignment coverage ranging from 76% to 80% with RMSD values ranging from 2.37 to 3.22 and TM-scores of 0.75 and 0.77 (Supplementary Fig. S4), indicating high structural similarity. In contrast, a comparison with AtMLO11 (Clade I) resulted in higher RMSD (3.68–4.02) and lower TM-scores (0.69–0.72) (Supplementary Fig. S5), reflecting reduced structural similarity. Collectively, these results suggested that Clade V-ii MLO proteins contain conserved sequence motifs and structural features with known PM susceptibility proteins, supporting their role as candidate susceptibility factors in coffee.

Analysis of MLOcis-regulatory elements

Cis-regulatory elements (CREs) in promoter regions play key roles in modulating gene expression during developmental and stress responses (Hernandez-Garcia and Finer 2014; Marand et al. 2023). To investigate the transcriptional regulation of coffee MLO genes, CREs within 1.5 kb upstream promoter regions were analyzed using the PlantCARE database (Lescot et al. 2002). A total of 44 CREs associated with stress responses, phytohormone signaling, and developmental processes were identified (data not shown) and classified into four categories: light-responsive, hormone-responsive, development-related, and stress-related. Overall, the CRE distribution and frequency were largely comparable across the 90 MLO promoters, with minor variations among the four coffee genomes (Supplementary Fig. S6). Among the identified CREs, Box 4 (ATTAAT) and the abscisic-acid responsive element (ABRE; CACGTG) were the most abundant across all four coffee genomes (Supplementary Fig. S6). Box 4 is associated with light-mediated development and stress responses (Wuthipong and Huehne 2015), while ABRE mediates ABA-dependent transcription via bZIP transcription factors (Hobo et al. 1999), linking MLO regulation to abiotic stress responses. In addition, the methyl jasmonate (MeJA)-responsive TGACG-motif, the salicylic acid (SA)-responsive TCA-element, and the anaerobic-related response (ARE) were also frequently detected in coffee MLO promoters (Supplementary Fig. S6). Stress-related CREs, such as ARE (anaerobic induction), DRE1 (Dehydration-responsive element; GCCGAC), LTR (Low-temperature response), MBS (MYB Binding Site), STRE (Stress response element; CCCCT or AG-rich), WER3 (WER-type MYB protein 3), and TC-rich repeats (defense and stress response), were widely distributed across coffee MLO promoters. In contrast, development-related motifs, including CAT-box, GCN4_motif, O2-site, and circadian elements, were less abundant but consistently present, suggesting roles in tissue-specific and temporal regulation.

The analysis of Clade V-ii MLO promoters revealed enrichment of ABRE, ARE (anaerobic responsive element), and TC-rich repeat (ATTCTCTAAC) (Fig. 4). While ABRE and ARE were common across all MLO promoters, TC-rich repeats were relatively rare overall (Supplementary Fig. S6) but specifically enriched in Clade V-ii promoters in Arabidopsis, tomato, and coffee (Fig. 4), indicating a specialized regulatory role in susceptibility-associated MLO genes. TC-rich repeats are known to be induced under defense and stress conditions and contribute to gene regulation during pathogen responses (Wen et al. 2015; Andolfo et al. 2019).

Fig. 4.

Cis-acting regulatory elements of Clade V-ii MLO gene promoters. Promoter regions were analyzed for cis-acting regulatory elements. The heatmap illustrates the abundance of each element, ranging from white (low) to red (high). Functional categories are indicated at the top, specific element names at the bottom, and gene names on the left.

Collectively, these results indicate that coffee MLO genes are regulated by complex, multilayered networks involving stress-and hormone-responsive CREs. The enrichment of ABRE, TGACG-motifs, and TCA-element suggests the integration of ABA-, MeJA-, and SA-mediated signaling pathways (Uno et al. 2000; Garretón et al. 2002), enabling the dynamic regulation of MLO expression under environmental and pathogenic stimuli.

Discussion

We undertook genome-wide identification and conducted a comparative analysis of the MLO gene family across four Coffea genomes, in this case two allotetraploid C. arabica cultivars and their diploid progenitors. The near doubling of the MLO gene number in C. arabica reflects its recent allotetraploid origin and retention of duplicated loci. Despite this expansion, the coffee MLO family remains highly conserved, as supported by phylogenetic relationships, conserved motifs, structural features, and strong purifying selection, indicating functional constraints and evolutionary stability. Nevertheless, variations in sequences, motif architectures, physicochemical properties, and regulatory regions, particularly between Geisha and Typica, suggest lineage-specific diversification following polyploidization. This potentially contributes to cultivar-specific responses to environmental and pathogenic stimuli.

Functionally, MLO proteins show clade-associated specialization, with PM susceptibility genes predominantly clustered in Clade V in dicots (Acevedo-Garcia et al. 2014; Kusch et al. 2016). Coffee MLO genes in the Clade V-ii subgroup exhibit strong similarity to known susceptibility genes in dicots, supporting their role as primary candidates mediating PM susceptibility. Together, these results indicate that Clade V-ii represents the most conserved susceptibility-associated subgroup of the coffee MLO gene family.

Given the limited availability of characterized PM resistance genes and natural sources of PM resistance in Coffea species, targeting MLO susceptibility genes may provide a promising strategy for developing PM-resistant coffee cultivars. In particular, the conserved evolutionary and regulatory features identified in Clade V-ii provide a basis for prioritizing candidate susceptibility genes for functional validation and genome-editing studies. Overall, this study provides an evolutionary and functional framework for understanding the diversification of the coffee MLO gene family and supports the development of MLO-based strategies for improving PM resistance in economically important Coffea species.

Supplementary Material

Supplementary materials are available at Horticultural Science and Technology website (https://www.hst-j.org).