We used 242 plant genomes, covering plant evolution from green algae to eudicots (Supplementary Table 1), to perform an extensive genomic and synteny analysis of the CesA superfamily in plants, of which the first step was the study of its genomic architecture. The 222 plant genomes with a BUSCO representation ≥75% contained 7,997 CesA superfamily genes, with an average of 36 members per species (CV = 48.6%; Supplementary Table 5; CV: Coefficient of Variation). The copy number of CesA/Csl genes correlates with both the ploidy level and the number of genome duplications of each species (ρ = 0.57 and ρ = 0.55, respectively; p < 0.001 for both; Supplementary Figure 1). Moreover, the CesA superfamily size increased considerably during plant evolution, with an acceleration at the rise of angiosperms. In fact, while charophytes, bryophytes, lycophytes, ferns, and gymnosperms contain a relatively similar number of CesA/Csl genes (∼18 per species, CV = 36.9%), this figure doubles in angiosperms despite similar inter-species variability (∼40 CesA/Csl genes per species, CV = 35.9%, p < .001; Supplementary Figure 2).

Cellulose synthase A is by far the most abundant CesA/Csl family, with 12 genes per species on average (CV = 43.1%; Supplementary Figure 3). This is roughly two times the size of the CslD, CslA, and CslC families (p < 0.001), which all display ∼6 genes per species. These three families are in turn significantly larger than the CslB, CslE, CslG, and CslM families (1-3 genes per species each; p < 0.001). Finally, the monocot-specific CslF and CslH families typically display 7 and 1 genes per species on average, respectively. Overall, the copy number of all the CesA/Csl families varies extensively across plant families, with standard deviations often equal to the mean size of gene families across species.

The relative positions of the CesA/Csl genes within each genome of the study were also assessed, highlighting two main patterns (Supplementary Tables 6–8). On the one hand, the CesA, CslA, CslC, and CslD families display a scattered genomic distribution among multiple singleton loci. On the other hand, the CslB, CslE, CslG, and CslM clades are usually organized in tandem arrays of 2–4 genes. These arrays were detected in both monocots and dicots, and synteny analysis showed their extensive conservation across species (see section “Gene Synteny”). A large conserved gene array was also found for the CslF family in all the Poaceae species evaluated [as previously reported by Schwerdt et al. (2015) for rice, sorghum, Brachypodium and barley only]. However, we also found 20–30% of the CslF genes of each Poaceae species organized as conserved singletons or tandems at separated loci. Interestingly, this second group of CslF genes was syntenic with CslD members of several eudicot species (see section “Gene synteny” and “Synteny Reveals the Genomic Trace of CslF Evolution”). Finally, the CslH and CslJ families displayed both configurations – tandem arrays and singleton loci. However, these families are present at too low frequencies and in too few species to allow generalizations.

The syntenic conservation of the CesA superfamily was also assessed in detail by using network synteny analysis. Results indicated that CesA superfamily genes are highly syntenic across diverse plants, more than what observed for other plant gene families, including regulatory genes. In fact, of the 7,005 CesA superfamily genes from the 193 genomes with a BUSCO representation ≥75% and at least five genes per scaffold on average, 6,262 (89%) were included into the CesA/Csl synteny network (Supplementary Table 9). This is a remarkably high proportion, 9 and 24% higher than what is found in comparable studies on the MADS-box and the APETALA2 gene families, respectively (Zhao et al., 2017; Kerstens et al., 2020). Furthermore, the synteny of the CesA/Csl genes is dense and extensive across diverse plant species and families. In fact, each CesA/Csl gene is syntenic to another 69 homologs from 46 different plant species on average. Moreover, the degree of gene synteny does not significantly differ between intra- and inter-family syntenic comparisons (p = 0.054, Figure 1A). Remarkably, for angiosperms only, syntenic conservation is even higher between than within plant families (p < 0.001, Figure 1B) and crosses the boundaries of monocots and dicots. This trend is clearly divergent from the patterns commonly observed for plant genes (Zhao and Schranz, 2019).

While synteny is dense and extensive across different plants, the same does not hold across different CesA/Csl families, which are organized in separate conserved genomic contexts. Accordingly, 99% of syntenic connections within the CesA/Csl synteny network takes place between genes from the same CesA/Csl family (Figure 1C). Moreover, each of the 48 syntenic communities identified by decomposing the CesA/Csl synteny network in groups of highly syntenic genes typically contains genes from only one CesA/Csl family (Figure 2). Furthermore, multiple communities were detected for most of the CesA/Csl families, revealing distinct conserved genomic contexts even at the level of gene subclades within CesA/Csl families (Figure 2). In this respect, the marked differentiation in the genomic organization of the CslA genes between monocots (especially Poaceae) and dicots is noteworthy (group B and parts of Groups A and C, Figure 2). Moreover, syntenic communities for specific plant families were detected for some CesA and CslA genes within Salicaceae, Fabaceae, and Brassicaceae (Group A of Figure 2). Next to these examples of intra-family sub-organization of gene architecture, we also found communities that cross the boundaries of single CesA/Csl gene families (Group C of Figure 2). Specifically, two communities, including 580 syntenic connections between CslF and CslD genes, appeared particularly relevant as they grouped ∼40% of all the syntenic connections involving CslF genes (Figure 2).

To further characterize the genomic patterns described above and to study their relevance for the evolution and sub-functionalization of the CesA/Csl genes, we investigated how such patterns relate to phylogenetic and functional gene data. For the CesA family, our results revealed striking correspondence between the genomic organization of genes, the syntenic conservation of gene architecture, the phylogenetic relationships between genes, and the diversification of genes into distinct isoforms with different functional properties. Specifically, the CesA phylogeny was divided into six main clades supported by moderate-to-high bootstraps, corresponding to the six main groups of CesA isoforms known in Arabidopsis thaliana and Oryza sativa (Figures 3A, 5). These six clades are grouped into two separated super-groups of three clades each, corresponding to the subdivision into primary and secondary CW CesA genes (Figures 3A, 5 and Supplementary Table 3). Interestingly, each of the six phylogenetic clades is differentially organized and differentially conserved at the genomic level by spanning only one of the six largest CesA syntenic communities found in the CesA/Csl synteny network (Figures 3B,D). Overall, these six clades and communities span 84% of the syntenic CesA genes and cover 93% of the angiosperm families included in the synteny analysis. The correspondence between CesA phylogeny and synteny is therefore striking, with only one major exception found for clade 1 of Figures 3A,B, which is subdivided into a main syntenic community largely dicot-specific (community 19 of Figure 3B) and other two smaller communities, of which one is specific to monocots (community 42 of Figure 3B) and one is a group of diverse CesA members (community 21 of Figure 3B). Besides this clade, only few other minor groups of genes specific to certain plant families deviate from the common phylogenetic and syntenic structure described above by being organized in distinct genomic contexts. Examples are a Fabaceae-specific community from clade 6 of Figure 3B (community 22 of Figure 3B), a Brassicaceae-specific community from clade 3 of Figure 3B (community 13 of Figure 3B), and a community from clade 4 of Figure 3B specific to Brassicaceae and Malvaceae (community 7 of Figure 3B).

Overall, the six different CesA phylogenomic clades are genomically independent from each other, with seldom inter-clade syntenic connections (Figures 3B,D). An exception is represented by the 236 syntenic links between clades 2 (primary CW CesAs) and 4 (secondary CW CesAs) of Figures 3B,C. A minor part of these connections is spread across different angiosperm families and can be regarded as background noise of the synteny analysis. However, 194 connections specifically involve the Malvaceae family and three cotton species therein: G. hirsutum, G. arboreum, and G. raimondii. These connections mainly involve three primary and five secondary cotton CesA genes from branches 2 and 4, respectively, which are syntenic to a total of 104 CesA genes from 78 different angiosperm genomes placed in opposite phylogenetic groups (branch 4 for the cotton CesAs from branch 2, and branch 2 for the cotton CesAs on branch 4) (Figure 3C). This suggests that the cotton secondary CW CesA genes involved in these connections changed their genomic position relative to the other angiosperm genes from the same phylogenetic group, ending up in a genomic context typical of primary CW CesAs (clade 2 of Figures 3B,C).

Remarkably, the CesA genes of G. hirsutum involved in the genomic inversion turned out to be associated with distinct functional properties at the basis of the massive fibre production of cotton. These are three different G. hirsutum CesA members: an already known CesA7, an already known CesA8, and a previously uncharacterized homolog of GhCesA8 (99.6% sequence identity) (Supplementary Tables 10, 11). Interestingly, these three isoforms have been proved to assemble together into functional rosettes for cellulose deposition in G. hirsutum fibre cells, with CesA8 specifically acting as enhancer for massive fibre deposition (Li et al., 2016, 2013). In addition, the biological mechanism of fibre synthesis involving G. hirsutum CesA7 and CesA8 appears conserved across different cotton species, including the ones included in our analysis and showing the same genomic inversion (Li et al., 2013). Accordingly, all the secondary CW CesA isoforms from G. arboreum and G. raimondii that are involved in the genomic context shift share 95.6% sequence identity with the G. hirsutum CesA7 and CesA8 for which functional data were available in scientific literature (Supplementary Table 11). Overall, these results suggest that the genomic positioning and organization of CesA genes might be critical for determining gene function, representing an important factor at the basis of cotton fibre deposition.

Striking correspondence between genomic organization, syntenic conservation, phylogenetic relationships, and functional diversification was also found within the CslD family. This family is structured into four main phylogenetic clades (Figure 4A). In addition, previous studies demonstrated the involvement of different CslD members into three main plant processes – pollen development, determination of vegetative organ size, and root hair formation – across several species (Wang et al., 2001; Kim et al., 2007; Bernal et al., 2008; Li et al., 2009; Luan et al., 2011; Hunter et al., 2012; Hu et al., 2018; Li et al., 2018). Strikingly, the functional subdivision of the CslD members overlapped with the main phylogenetic CslD clades (Figure 4). In fact, known CslD genes involved in determination of vegetative organ size and root hair formation were only found within clades 1 and 4, respectively, while CslD members known to be involved in pollen development were only observed in clades 2 and 3 (Figure 4). In turn, the different phylogenetic and functional CslD groups corresponded to distinct and highly conserved CslD syntenic communities, highlighting their independent organization and conservation across plants (Figure 4B). However, while one syntenic community was found for each of the phylogenetic clades grouping CslD genes involved in organ size determination and root hair formation, respectively, multiple communities specific to different taxonomic groups were detected within the CslD branch associated to pollen development (Figure 4). Of these, one is broad and spans both monocot and dicot species (community 39), one is large but restricted to dicots (community 41), and four are minor and restricted to specific taxonomic clades.

While the subdivision of CslD genes described above is shared by all angiosperms, the same does not hold for earlier land plants. In fact, CslD members from bryophytes and lycophytes are found in only a single phylogenetic group placed between the CslD branches associated to organ size determination and pollen development. Moreover, fern CslDs are divided into two groups that are closest to angiosperm CslD members involved in organ size determination and root hair formation, respectively. Finally, CslD genes from bryophytes, lycophytes, ferns, and gymnosperms do not display any syntenic connection with the angiosperm ones.

Commonalities between phylogenetic, syntenic, and functional patterns were also found in other Csl families. However, for the evolutionarily close CslB/H/G/J/M/E families, such commonalities are observed for whole families rather than gene subclades within families (Supplementary Figure 4), with deviations observed only in a CslM syntenic community specific to Fabaceae (community 23 of Supplementary Figure 4) and a small CslE community with only Capsicum genes (community 17 of Supplementary Figure 4). Finally, CslA and CslC genes also display syntenic suborganization of their phylogenies (Supplementary Figure 5). However, the functional information available for these families are not abundant and do not reveal any clear correspondence with the phylogenetic and/or syntenic structure observed.

Phylogenetic and synteny data were used to also study the evolution of the main CesA/Csl families. For the CesA genes, the taxonomic profiling of the six main CesA clades (see section “Distinct Phylogenomic Features for Distinct CesA Isoforms”) revealed that CesA sequences of bryophytes and lycophytes are phylogenetically closest to the primary CW CesA branches, and specifically to the primary CW CesA genes homolog to the redundant Arabidopsis CesA2/5/6/9 and the O. sativa CesA3/5/6 (clade 1 of Figure 5). This observation suggests that primary CW CesAs are the oldest CesA sequences of land plants, and that early land plants likely only had primary CW CesA (-like) genes. The positioning of fern CesAs suggests that the diversification of CesA genes toward the six main angiosperm clades described in section “Distinct Phylogenomic Features for Distinct CesA Isoforms” started sometime during the late lycophyte or early fern evolution. In fact, while a group of fern CesAs is also positioned close to the bryophyte and lycophyte sequences of clade 1 of Figure 5, another group of fern genes precede the split of the other two clades of angiosperm primary CW CesAs (clades 2 and 3 of Figure 5). Moreover, a third group of fern CesAs is closest to the three groups of angiosperms secondary CW CesAs (clades 4, 5, and 6 of Figure 5), just preceding their diversification. Therefore, CesAs diversification toward the different phylogenomic groups observed in higher plants significantly progressed during fern evolution. However, gymnosperms are the first group of plants whose CesA sequences are found in all the six CesA clades genomically conserved across all the angiosperms. To conclude, data on gene copy-number within each of the six CesA clades also support the evolutionary model discussed. In fact, the early-diverging clade 1 contains by far the most CesA sequences in both Figures 3A, 5 (771 sequences in the full CesA tree of Figure 3A), followed by clades 2 and 3 of primary CW CesAs (394 genes on average, CV = 4.3%), and finally by the three secondary CW CesA clades (273 genes each on average, CV = 6.3%). These data, together with the presence of redundant Arabidopsis CesA copies within two primary CW CesA clades (1 and 2 of Figure 5), agree with the view that oldest CesA families evolved first into multiple isoforms, of which some underwent full sub-functionalization within the cellulose synthesis machinery [a model known as constructive neutral evolution, recently advanced by Haigler and Roberts (2019) for the CesA genes – see discussion].

In the CslF family, our phylogenomic analysis detected 580 syntenic connections between the CslF genes organized as singletons or tandems in the Poaceae genomes (see section “Gene Distribution Along Genomes”) and a subset of CslD members, all included within the same syntenic community (community 35 in Figure 6 and green dots in Figure 2). Interestingly, 327 of these connections involve CslD genes from 53 different eudicot species (27 different plant families) that display simultaneous synteny with a grass CslD gene and a grass CslF member in 17 of the 21 Poaceae genomes of our study (Figure 6B). The frequency of these connections is negatively correlated with the evolutionary distance of the species harbouring “seed” CslD sequences from Poaceae (r = –0.57). Moreover, the “seed” CslD genes all belong to the CslD phylogenomic clade grouping CslD genes involved in root hair formation (see section “Gene Distribution Along Genomes”). All together, these observations reveal the genomic trace of the CslF origin in grasses. On the one hand, they confirm that the CslF genes originated through the duplication of some CslD members (Yin et al., 2009; Schwerdt et al., 2015; Little et al., 2018). On the other hand, they show that the CslF family is nested within the CslD one (and specifically within the CslD clade involved in root hair formation). This in turn proves that the CslF genes took origin when the CslD family was already formed, thanks to a relaxed evolutionary pressure on duplicated CslD members in grasses. To conclude, it is noteworthy that no synteny was detected between the CslF genes syntenic to CslDs (red arrows in Figure 6) and the other CslF members included in community 43 (Figure 6), which form the conserved CslF genomic array (see section 2.1.2). The fact that only some CslF genes – corresponding to one of the two architectural configurations of this family – display synteny with CslD members questions whether the CslF family originated through a single CslD duplication, or if the family underwent different rounds of expansion in Poaceae. In this respect, the phylogenomic positioning of the barley CslF members recently shown to synthesize (1,4)-β-linked glucoxylans was checked in the tree of Figure 6A (blue arrow). Interestingly, both HvCslF3 and HvCslF10 are positioned within the phylogenetic clade spanned by community 43, which does not display synteny with CslD members.