For a subset of 22 Cleomaceae species selected for the present study, comprising herbaceous and shrub species, seeds were collected in situ from 20 locations distributed in 11 Brazilian states (Figure 1 and Supplementary Table S1). The species used in this study were sampled from four Brazilian biomes (Amazon, Atlantic Forest, Cerrado, and Caatinga, or transitional areas between these two last) (Figure 1 and Supplementary Table S1). These biomes differ from each other mainly by the predominant climate and vegetation. Thus, the Amazon Forest has a hot and humid climate, with dense vegetation characterized by large trees. The Cerrado (Brazilian Savanna), however, has a predominantly seasonal tropical climate, with periods of rain and drought, and vegetation characterized by gnarled trees, grasses and shrubs. The Caatinga (Brazilian semi-arid region), in turn, has a semi-arid tropical climate and medium-sized shrubby vegetation, with twisted branches and leaves adapted to periods of drought. The Atlantic Forest, on the other hand, has a predominantly tropical-humid climate with high temperatures and rainfall and vegetation marked by the presence of large and medium-sized trees forming a dense and closed forest (IBGE, 2021).

After morphological identification, seven accessions were classified as T. hassleriana (shrub species), collected from different environments [locations: Viçosa, São Miguel do Anta, Canaã, and Piau in the state of Minas Gerais (MG); Domingos Martins in the state of Espírito Santo (ES); and Joinville and Canoinhas in the state of Santa Catarina (SC)]. Six accessions were classified as T. longicarpa (shrub species), collected from different environments [locations: Manaus in the state of Amazonas (AM); Ibimirim, Afrânio and Arcoverde in the state of Pernambuco (PE); Picos in the state of Piauí (PI); and Lavras in the state of Ceará (CE)]. The other species were classified as C. paludosa (shrub species) and T. microcarpa (herb species), which were collected from Belém in the state of Pará (PA); T. aculeata (herb species) and T. diffusa (herb species) collected from Feira de Santana in the state of Bahia (BA); T. rosea (shrub species) collected from Colatina (ES); T. spinosa (shrub species) collected from Teresina (PI); T. parviflora (shrub species) collected from Pombal in the state of Paraíba (PB); T. siliculifera (shrub species) collected from Rio Pardo (MG); and G. gynandra (herb species) collected from Mossoró in the state of Rio Grande do Norte (RN) (Figure 1 and Supplementary Table S1). The exsiccates are deposited in the VIC herbarium of the Universidade Federal de Viçosa (UFV).

The studied genera can be easily distinguished based on their morphological traits. Cleoserrata species are herbs with a glabrous surface, bearing leaves with three leaflets but without stipules, ebracteate inflorescences, and large smooth black seeds. Gynandropsis species are annual herbs with a puberulent-glandular surface, bearing leaves with three to five leaflets, ebracteate inflorescences, and small smooth black seeds. Meanwhile, Tarenaya species are herbs or shrubs with a puberulent-glandular or pubescent surface, bearing three to seven foliolate leaves with epinescent stipules, bracteate inflorescences, and seeds with longitudinal and transverse streaks that are not very well developed (Stevens, 2001; Bercht and Presl, 2021).

Furthermore, the different Tarenaya species in Brazil can also be easily distinguished from one another. T. aculeata bears tree foliolate leaves, small white petals, and slightly moniliform capsules; T. diffusa bears three to five foliolate leaves, small white petals, and fusiform capsules; T. microcarpa bears three to five foliolate leaves, glandular-puberulent indumentum, and small purplish petals; T. hassleriana bears armed petioles, glandular-puberulent indumentum, spines on the veins, and pink flowers; T. parviflora bears glabrescent indumentum, three to seven foliolate leaves, armed petioles, and white to deep pink petals; T. longicarpa bears glandular-puberulent indumentum, completely white petals or petals white at base and purple at the apex (depending on accession), and long capsules, with the fruit being two to three times longer than the gynophore (which differs from T. spinosa); T. spinosa bears glandular-puberulent indumentum and small white petals (compared with T. longicarpa); T. rosea bears glandular-puberulent indumentum and petals white at base and pink at apex; and T. siliculifera bears glandular-puberulent indumentum without stipular spines, small white petals, obovoid capsules (almost all species of the genus have cylindrical fruits). The different accessions of T. hassleriana and T. longicarpa could also be distinguished from one another based on height, hairiness, spine density and disposition, and inflorescence as well as flower size.

The collected seeds were germinated in 5 L plastic containers, filled with a commercial substrate supplemented with 14 g of NPK (4:14:8) per pot [0.56 g of N, 0.86 g of P₂O₅, and 0.93 g of K₂O]. During spring and summer, five plants of each accession were grown in a greenhouse under semi-controlled conditions (maximum photosynthetically active radiation of 1,500 μmol photons m^–2 s^–1, relative air humidity of 60%, and temperature of 30 ± 2°C) and were irrigated daily. Under these conditions, the plants were grown for 2 (herbaceous species) to 5 months (shrubby species) until physiological characterization and harvesting of leaf material for anatomical and biochemical analyses.

Total DNA was extracted from 22 fresh samples using a modified CTAB method (Doyle and Doyle, 1990). For phylogenetic analyses, the ITS regions were sequenced using the universal primers ITS 1 and ITS 4, as previously described (White et al., 1990). The fragments of interest were amplified using polymerase chain reaction (PCR). PCR products were visualized using 1% agarose gel electrophoresis, isolated, and cleaned with the Wizard SV Gel and PCR Clean-Up System (Promega). The purified fragments were then sequenced (Myleus; Belo Horizonte, Brazil). The electropherograms were visually assessed, edited, base-called, assembled, and manually aligned using Sequencer 4.1 (GeneCodes, Ann Arbor, MI, United States).

We, therefore, performed a comprehensive phylogenetic analysis by adding the ITS sequences from Cleomaceae species, recovered from the National Center for Biotechnology Information (NCBI) (Supplementary Table S2), which has shown promise for the family in previous works (see Patchell et al., 2014; Reeves et al., 2018; Soares-Neto et al., 2020). Based on these sequences, we generated a matrix with 135 taxa (648 bp in length). The results were interpretated using Bayesian inference (BI). BI analyses were performed in MrBayes using the best-fit evolutionary model, selected according to the Akaike information criterion (Posada and Buckley, 2004), with MrModeltest 2.2 (Nylander, 2004). The Markov chain Monte Carlo algorithm was executed using four runs with 50 million generations each and sampling every 5,000 generations. The first 25% trees were discarded as burn-in. The remaining trees were used to construct the majority rule consensus tree, and the posterior probability (PP) for each node was calculated.

Leaf fragments (∼2 cm²) from each Cleomaceae accession (sample) and Solanum lycopersicum (internal standard, 2C = 2.00 pg; Praça-Fontes et al., 2011) were co-chopped in 0.5 mL of OTTO-I nuclear extraction buffer (Otto, 1990) supplemented with 2 mM dithiothreitol and 50 μg mL^–1 RNAse (Praça-Fontes et al., 2011). Then, 0.5 mL of the same buffer was added, and the suspensions were filtered through a 30 μm nylon mesh, placed into a microtube, and centrifuged at 100 × g for 5 min. The precipitate was resuspended in 100 μL of OTTO-I buffer and incubated for 10 min. The suspensions were stained with 1.5 mL of OTTO-I:OTTO-II (1:2) buffer supplemented with 2 mM dithiothreitol, 75 μg mL^–1 propidium iodide, and 50 μg mL^–1 RNAse (Praça-Fontes et al., 2011). The suspensions were incubated in the dark for 30 min, filtered through a 20 μm nylon mesh, and analyzed using a flow cytometer (BD Accuri C6, Accuri Cytometers, Belgium) equipped with a 488 nm laser source to promote propidium iodide emission at FL2 (615 – 670 nm) and FL3 (>670 nm). Fluorescence peaks of the G₀/G₁ nuclei of each Cleomaceae accession and S. lycopersicum were analyzed based on histograms using BD Accuri™ C6. The G₀/G₁ nuclear peaks exhibiting a coefficient of variation of ≤5% were considered for genome size measurement, calculated as 2C_D = (C1/C2) × 2C_S, where 2C_D is the 2C value (pg) of each Cleomaceae accession, C1 is the mean G₀/G₁ peak channel of the Cleomaceae accession, C2 is the mean G₀/G₁ peak channel of S. lycopersicum, and 2C_S is 2.00 pg of S. lycopersicum. Based on previous results, the S. lycopersicum internal standard was replaced with Raphanus sativus ‘Saxa’ (2C = 1.13 pg, Praça-Fontes et al., 2011) to measure the 2C value of G. gynandra (GG, Mossoró-RN) and T. siliculifera (TSI, Rio Pardo-MG).

Leaf material was oven-dried at 60°C for 48 h and used to determine δ¹³C and the C/N ratio. The material was analyzed using the Isoprime 100 isotope ratio mass spectrometer coupled to an isotope cube elemental analyzer (Elementar, Hanau, Germany), as described by Gowik and Bra (2011).

Fully expanded leaves from the third node of non-flowering plants were used for gas exchange and chlorophyll a fluorescence analysis, which were carried out with an open-flow infrared gas exchange analyzer system equipped with an integrated fluorescence chamber (IRGA, LI-COR Inc. LI-6400XT; NE). The measurements were performed during the light period between 8 and 12 h (solar time), using a 2 cm² leaf chamber at 25°C, flow rate of 400 μmol CO₂ s^–1 and saturating light intensity of 1,000 μmol photons m^–2 s^–1 at the leaf level. The leaf-to-air vapor pressure deficit was maintained at 1.2–2.0 kPa, and the amount of blue light was set to 10% PPFD to optimize the stomatal aperture. Corrections for CO₂ leakage into and water vapor from the leaf chamber of LI-6400 were applied to all gas exchange data, as previously described (Rodeghiero et al., 2007). Dark respiration (R_d) was measured using the same gas exchange system described above, on the same leaves used for determination of the gas exchange parameters, after at least 1 h of acclimation in the dark (Rodeghiero et al., 2007).

Fully expanded leaves from the third node of non-flowering plants were harvested in the middle of the light period, soon after gas exchange and fluorescence analyses. At this point, the herbaceous plants were 2-month-old and the shrub species were 5-month-old. The samples were snap-frozen and stored at –80°C until analysis. Metabolites were extracted by grinding the material in liquid nitrogen and adding the appropriate extraction buffers. Approximately 25 mg of the freeze-dried leaf matter was subjected to ethanolic extraction at 70°C for 30 min. After centrifugation (at 21,800 × g for 5 min), the total chlorophyll content (a + b) and chlorophyll a/b ratio were determined as previously described (Porra et al., 1989).

The concentration of glucose, fructose, and sucrose was determined from the liquid phase, as described previously (Fernie et al., 2001). Total free amino acid (Cross et al., 2006), malate, and fumarate (Nunes-Nesi et al., 2007) contents were determined as previously described. From the ethanol insoluble fraction, the concentration of starch (Fernie et al., 2001) and total soluble proteins (Bradford, 1976) was determined according to standard protocols.

Fully expanded source leaves from each sampled plant were collected in the middle of the light period, wrapped in aluminum foil, flash frozen in liquid nitrogen, and stored at –80°C until analysis. After lyophilization and grinding, approximately 15 mg dry weight was extracted for gas chromatography–mass spectrometry (GC–MS) in 1.5 mL of extraction solution (water:methanol:chloroform, 1:2.5:1 v/v) as described by Fiehn (2007), using 5 μM ribitol (Sigma-Aldrich) as the internal standard. The extract (15 μL) was dried in a vacuum concentrator and derivatized in two steps using the MPS Dual Head autosampler (Gerstel), as described by Gu et al. (2012). After incubation for 2 h at room temperature, the samples were analyzed as described previously (Shim et al., 2020) using the 5977B GC-MSD system (Agilent Technologies). The peaks were integrated using MassHunter Quantitative (v b08.00; Agilent Technologies). For relative quantification, metabolite peak areas were normalized to the corresponding dry weight and peak area of the internal standard.

Anatomical analysis was conducted to evaluate the VD of fully developed leaves in the upper third. Diaphanization of the plant material was performed as described previously (Zsögön et al., 2015). Subsequently, the samples were assembled on glass slides, and images of the adaxial epidermis were obtained using the Zeiss Axio Scope A1 photomicroscope coupled to a color image capture system (Axiovision ^® 105). Five photographs of each slide were obtained to cover the entire leaf length. In addition, five replicates (five leaves per individual) were collected and analyzed. Vein density was calculated as the length of the second-, third-, and fourth-order veins in a given leaf area (mm mm^–2). Quantification was performed using Image Pro-Plus ^® (v 4.5; Media Cybernetics, Silver Spring, MD, United States). In addition, the median region of the fully expanded leaves from the third node was collected and fixed in formaldehyde:acetic acid:ethanol 50% (FAA50) at a ratio of 5:5:9 (v/v) for 48 h. The material was dehydrated in an ethanolic series and embedded in historesin (Leica, Heidelberg, Germany). Transverse sections (thickness, 5 μm) were obtained using a self-advancing rotary microtome (RM 2155, Leica, Heidelberg, Germany) with glass razors. The sections were stained with 0.05% toluidine blue in acetate buffer (pH 4.7) (O’Brien et al., 1964) for 1 min and assembled between slides and coverslips with synthetic resin. Photomicrographs were obtained using a photomicroscope (Zeiss, MC-80). The micrometric scales were photographed and enlarged under the same optical conditions. Leaf thickness was measured considering the entire length of the epidermis. For each slide analyzed, five measurements were obtained in different regions to cover the entire area. The distance between the vascular bundles was measured between the sides of the veins. Bundle sheath cell width (BSCW) was measured radially in five cells at random on each slide (cells from different vascular bundles were selected when possible). Leaf thickness, palisade parenchyma, spongy parenchyma, intervein distance, and BSCW were measured (Supplementary Figure S1) using Image Pro-Plus ^®.

Data were obtained from five plants per accession, each placed in an individual pot; thus, each pot represented a biological replicate. All plants were completely randomized. The effect of accession was determined using analysis of variance (P < 0.05). Means were analyzed using Tukey’s test.

Phylogenetic generalized least squares (PGLS) regression was applied. We used a phylogenetic tree constructed using BI based on the ITS sequences of 21 novel Brazilian accessions (Supplementary Figure S2) and a matrix of 74 variables. The species G. gynandra (C₄ species) was removed from the analysis, as its genome size values (2–4 times greater compared to the other species in this study) and results obtained from other analyzes could bias/influence the regression result. The script available at https://lukejharmon.github.io/ilhabela/instruction/2015/07/03/PGLS/ was used. PGLS uses the common statistical mechanism of generalized least squares and applies it to phylogenetic comparative data (Felsenstein, 1985a). This method allows for the estimation of the impact of phylogeny on the covariance among residuals, thereby controlling for relatedness. In the case of Brownian motion, residues with variances and covariances follow the structure of the phylogenetic tree (Felsenstein, 1985b). All statistical analyses were performed using Statistica and R.

In flow cytometry histograms with G₀/G₁ peaks of nuclear suspensions, the coefficient of variation for all Cleomaceae species/accession and the internal standards (S. lycopersicum or R. sativus) was < 3.70%. Based on genome size, the 22 accessions were separated into seven groups: (i) T. hassleriana accessions (THC, THCS, THDM, THJ, THP, THS, and THV) with 2C = 0.55 pg ± 0.003; (ii) T. diffusa (TD) with 2C = 0.59 pg ± 0.004; (iii) T. aculeata and T. microcarpa (TA and TM) with 2C = 0.66 pg ± 0.004; (iv) T. siliculifera (TSI) with 2C = 0.77 pg ± 0.003; (v) C. paludosa (CP) with 2C = 1.08 pg ± 0.004; (vi) nine Tarenaya species (TS, TP, TL, TR, TAF, TAM, TARC, TC, and TIB) with 2C = 1.30 pg ± 0.014; and (vii) G. gynandra (GG) with 2C = 2.20 pg ± 0.002 (Figure 2A). This pattern found for genome size was somehow unexpected, mainly because most species belong to the same genus. Accordingly, indicating the high degree of variation, mainly of Tarenaya species/accessions, which are spread into five different groups. Thus, we observed intrageneric and interspecific variation in nuclear genome size among the Brazilian Cleomaceae species. Based on the mean 2C values and phylogeny (Figure 2B), we speculate that increase and decrease in genome size contributed to the diversification of Cleomaceae species, specifically the C₄ species G. gynandra (2C = 2.20 pg, fourfold higher than that of T. hassleriana). Considering the correlation between the 2C value and the 2n chromosome number of Tarenaya species and G. gynandra, we hypothesize that euploidy and disploidy (ascendant and/or descendant) resulted in interspecific variations in the 2C value and C₄ photosynthetic mechanism of G. gynandra. Alternatively, genome size differences may be attributable to whole genome duplications (Mabry et al., 2020) and/or changes in repetitive DNA content (Beric et al., 2020, 2021).

Our Bayesian phylogenetic reconstruction based on the ITS sequences of Cleomaceae species indicated that most of our analyzed accessions were positioned in the Tarenaya genus. The exceptions were the sequences of C. paludosa and G. gynandra, which were placed into the Cleoserrata and Gynandropsis clusters, respectively (Figure 2B). Moreover, Tarenaya is sister to two American genera, namely Cleoserrata and Melidiscus, whereas Gynandropsis is sister to the African cluster comprising Coalisina and Iltisella as well as to Cleomella (Figure 2B). These results are somehow contrary to the findings reported by Bayat et al. (2018), who used five molecular markers (ITS, matK, ndhF, ycF, and rps) and another set of species/sequences. Interestingly, Gynandropsis, Areocleome, and Coalisina, which harbor species with C₄ photosynthetic metabolism, formed different clusters (Figure 2B), reinforcing that this photosynthetic mechanism has evolved independently at least three times within the family, as was previously reported (Patchell et al., 2014; Bayat et al., 2018).

The analyzed genera, namely Gynandropsis, Cleoserrata, and Tarenaya are monophyletic. Within the Tarenaya group, the series Rosea, Parviflorae, and Aculeatae are monophyletic; the Spinosa series, which includes the type species of the genus (T. spinosa), is polyphyletic; and the Aculeatae series, which includes herbaceous species bearing small flowers, is placed more externally to the other series of the genus (Figure 2B). Within the Spinosa series, T. hassleriana accessions form different sub-clusters, indicating high genetic diversity among the species (Figure 2B). Overall, the 20 accessions of Tarenaya studied here cover a great genetic diversity into this genus, mainly considering that species from all series were sampled and were represented in all analysis.

The carbon isotope composition of plants has been widely used to identify the photosynthetic mechanisms, and it can, therefore, assist phylogenetic studies exploring the evolution of C₄ photosynthesis (Von Caemmerer et al., 2014; Yang et al., 2017). Hence, to help identifying the photosynthetic mechanism of the accessions studied here and to merge these data with molecular phylogeny, we performed carbon isotope composition (δ¹³C) analysis using fully expanded leaves of plants growing under the same conditions. In the Cleomaceae accessions tested, the δ¹³C value ranged from –15‰ to –32‰, representing at least two types of photosynthetic mechanisms, namely C₃ and C₄ (Figure 2C). As expected, the Brazilian accession of G. gynandra exhibited the highest δ¹³C value (–15‰), and the other species showed values around –32‰. Meanwhile, T. siliculifera and T. longicarpa (Lavras, CE) exhibited δ¹³C values close to –29‰, being significantly different from the other accessions (Figure 2C), highlighting a certain degree of variation among the C₃-like species found into the genus Tarenaya.

In order to determine the correlations of genome size with the other study variables, independently, in the phylogenetic topology of the Brazilian Cleomaceae accessions, we decided to perform PGLS regression analysis (Table 1 and Supplementary Figures S3, S4). As the size of the G. gynandra genome is 2–4 times greater compared to the other species in this study, we removed this species from the analysis to avoid bias/influence the regression result.

It is noteworthy that a phylogeny, based on the studied species (except G. gynandra), was built exclusively for PGLS. Thus, with a smaller sample of species and a single molecular marker (ITS), it is possible to see that the Tarenaya series were no longer monophyletic (Supplementary Figure S2). PGLS analysis revealed that 32 out of 74 variables analyzed here were significantly correlated to genome size (P < 0.05) (Table 1). Based on the significant regressions, the models with gas exchange parameters (A_N, A_gross and R_D), leaf anatomical traits (size of palisade parenchyma and VD), organic acids (malate, isocitrate, and gluconate), amino acids (aspartate, α-alanine, methionine, phenylalanine, tyrosine, glutamine, serine, glycine, and leucine), other nitrogen-metabolism related compounds (γ-aminobutyric acid, GABA), putrescine and ornithine), total clorophyll, glucose, and glycerol were positively correlated with genome size. Interestingly, the δ¹³C value, intrinsic water use efficiency (WUEi), malonate, myo-inositol, proline, starch content, proteins content, seeds weight and specific leaf area were negatively correlated with genome size (Table 1 and Supplementary Figures S3, S4).

To completely characterize and understand the physiological basis for the differences among the studied Cleomaceae accessions, we performed detailed physiological analyses of diffusional, photochemical, and biochemical constraints on photosynthesis (Figures 3A–D and Supplementary Figure S5). Under the tested environmental conditions (as described in the section “Materials and Methods”), Cleomaceae accessions showed natural variation, even within the same species, in terms of physiological parameters related to photosynthesis. In general, the studied Cleomaceae accessions exhibited a high A_N, with values greater than 20 μmol CO₂ m^–2⋅s^–1 (Figure 3A and Table 2). However, the herbaceous accessions of T. aculeata, T. diffusa, and T. microcarpa exhibited the A_N of ∼10 μmol CO₂ m^–2⋅s^–1. While the shrubby accession of T. siliculifera showed the A_N of <10 μmol CO₂ m^–2⋅s^–1. T. siliculifera exhibited the lowest stomatal conductance (g_s), followed by G. gynandra, herbaceous species (T. aculeata, T. diffusa, and T. microcarpa), and Tarenaya species (Spinosa I and II, Parviflorae, and Cleorosea) (Supplementary Figure S5). For the internal CO₂ concentration (C_i), which is related to g_s, the same pattern as that for g_s was observed (Supplementary Figure S5). Regarding transpiration (E) and R_d, the studied Cleomaceae accessions exhibited substantial differences with no clear pattern between the groups (Figure 3 and Supplementary Figure S5). Similarly, WUE_i was remarkably higher in G. gynandra, albeit without a clear trend in the rest of the species (Figure 3).

To unveil the natural variations among the different Cleomaceae accessions studied, we performed detailed metabolite analyses using leaf samples. First, we examined the levels of the major carbon- and nitrogen-containing metabolites in leaves harvested at the middle of the light period. Overall, despite clear variations in physiological traits (Figures 3E,F and Supplementary Figure S6), the steady-state levels of soluble sugars (glucose, fructose, and sucrose), starch, amino acids, proteins, and chlorophylls did not differ among the plants of different species (Figures 3E,F and Supplementary Figure S6). The only exception to this trend was T. siliculifera, which exhibited higher steady-state levels of sugars and starch than the other species (Figure 3E and Supplementary Figure S6). Meanwhile, G. gynandra exhibited higher levels of amino acids (Figure 3E) and proteins than T. diffusa (Figure 3F). There were no significant differences in the chlorophyll a/b ratio among the studied species (Supplementary Figure S6).

Furthermore, we performed metabolite profiling of leaf samples from all studied accessions using GC–MS. This analysis allowed for the identification and quantification of 33 metabolites (Supplementary Table S3 and Supplementary Figure S7). For 19 metabolites analyzed, we did not observe significant differences among the studied accessions. Among these metabolites, there were no differences in the levels of the amino acids aspartate, glycine, and serine or the organic acids citrate and isocitrate. In contrast, T. siliculifera exhibited significantly higher levels of myo-inositol, succinate, malonate, malate, threonate, glycerate, glycerol, glucose, fructose, and sucrose. In addition, T. hassleriana exhibited higher levels of methionine, glutamate, and alanine.

In leaves, the formation of venation, plasmodesmata, and other barriers to gas diffusion is important for photosynthetic efficiency. Thus, the variability in leaf anatomical traits, which are related to the carbon concentrating mechanism, was determined in the Cleomaceae accessions studied. Leaf anatomical parameters were evaluated by diaphanization and transverse sectioning to determine VD and BSCW. T. hassleriana (THC, THCS, THP, and THS), T. longicarpa (TP, TC, and TIB), T. rosea, and T. spinosa showed higher VD than the other analyzed accessions (Figure 3G). Meanwhile, T. hassleriana (THDM and THV), T. longicarpa (TAF, TAM, and TARC), and T. parviflora showed intermediate VD, whereas G. gynandra, T. hassleriana (THJ), T. diffusa, T. aculeata, T. microcarpa, T. siliculifera, and C. paludosa exhibited lower VD (Figure 3G).

Bundle sheath cell width was higher in G. gynandra, T. hassleriana (THC and THCS), T. siliculifera, and T. longicarpa (TAM) than in the other accessions studied (Supplementary Figure S8). Meanwhile, T. diffusa, T. aculeata, T. microcarpa, T. parviflora, and T. longicarpa (TARC and TAF) exhibited lower BSCW. The remaining 11 accessions exhibited intermediate BSCW (Supplementary Figure S8).

The interveinal distance (ID) was higher in T. microcarpa but lower in G. gynandra and T. siliculifera. The other accessions showed intermediate ID (Supplementary Figure S8). T. hassleriana (THP) showed the highest leaf thickness, whereas T. diffusa and T. microcarpa showed the lowest leaf thickness. The other accessions exhibited intermediate values and did not differ in terms of thickness (Supplementary Figure S8). Regarding the thickness of the palisade (Figure 3H) and spongy parenchyma (Supplementary Figure S8), T. hassleriana (THP) exhibited higher values and G. gynandra showed lower values than the remaining accessions, which showed intermediate values (Supplementary Figure S8).

Surprisingly, we could not separate the accessions into groups based on these morphoanatomical characteristics, which would be expected considering that the C₄ species G. gynandra could be clearly distinguished from the other C₃ species in terms of most of the analyzed parameters (e.g., VD and BSCW) (Supplementary Figures S9, S10). In this context, our results suggest that the analyzed accessions diversified at the species level. Since T. hassleriana and T. longicarpa differed in terms of the analyzed parameters, anatomy would provide additional evidence of this diversification. Although these traits indicate the possible diversification in photosynthetic mechanisms, physiological and biochemical data provide additional evidence.

Cleomaceae is a highly diverse family, both morphologically (e.g., monoecious, dioecious, polygamous species; e.g., Omondi et al., 2017; Zohoungbogbo et al., 2018; Riaz et al., 2019) and physiologically (species with different photosynthetic metabolisms: C₃, C₃–C₄, and C₄). This diversity was represented by species studied here, which were morphologically (Tarenaya, Cleoserrata, and Gynandropsis; Stevens, 2001; Bercht and Presl, 2021) and physiologically (e.g., Tarenaya and Gynandropsis; Marshall et al., 2007; Bayat et al., 2018) distinct. The three genera analyzed here, namely Tarenaya, Cleoserrata, and Gynandropsis, were monophyletic (Figure 2). Of note, however, Gynandropsis is a monotypic genus (G. gynandra). Tarenaya, one of the most diverse genera in terms of the number of species within Cleomaceae, also includes many morphologically distinct species, which allows for the separation of the genus into series (small groups of species) (e.g., Iltis, 1957, 1959; Hall et al., 2002; Sanchez-Acebo, 2005; Inda et al., 2008; Feodorova et al., 2010; Iltis and Cochrane, 2014). These series have previously been shown to be monophyletic (see Soares-Neto et al., 2020), and our results are largely congruent with this report, with the exception of the Spinosa series, which we found to be polyphyletic. Of note, our study is limited by the number of samples of this rather diverse series (see Soares-Neto et al., 2020). Thus, studies with increased taxon sampling or markers are warranted to resolve the monophyly of the Spinosa series.

The morphological and physiological diversity can be linked to variations in genome size, among other factors (Paterson et al., 2010; Roddy et al., 2019). Thus, an interspecific variation at the diploid level, for example, as found for the species T. hassleriana, T. diffusa, T. aculeata, and T. microcarpa, may suggest that the DNA content may be a parameter that can be used to differentiate the species (e.g., Miscanthus 2C = 3.9 to 6.9 pg – Sheng et al., 2016). However, the observed species and genome size diversity was not consistent with the molecular phylogenetic results (Figure 2; Patchell et al., 2014; Bayat et al., 2018). This pattern is different, for example, from that observed for Lupinus (Fabaceae). In this genus, the genome size (0.97–2.44 pg) data supports the generally accepted taxonomic classification of the Old World lupins (Naganowska et al., 2003). In this sense, considering members of Cleomaceae, we observed that different series have species with the same genome size (Parviflora, Rosea, and Spinosa). Furthermore, within the same series it was also possible to observe representatives with different genome sizes (e.g., T. hassleriana and T. spinosa – Spinosa series, and T. rosea and T. siliculifera – Rosea series) (Figure 2). It is noteworthy that the studied Cleomaceae species have an average genome size of 0.5 – 2.2 pg, which is larger than that recorded for Arabidopsis. However, in line to most species that have their genome quantified (Bai et al., 2012), such as rice, tomatoes and fruit–cherry, mango, papaya, orange, and peach (Arumuganathan and Earle, 1991).

Based on the variations in genome size among Cleomaceae genera and even among species of the same genus (e.g., Tarenaya; Figure 2A), even our small sample of the family reflects the relatively large interspecific variability. The intraspecific variations in the ITS sequence as well as the interspecific variations in this sequence and the 2C value may be related to different selective environmental pressures. Polyploidy, disploidy (ascendant and/or descendant), and genome size are associated with life-history traits, including vegetative form, flowering time, and adaptation to particular ecological niches, generating a basis for evolutionary novelty (Gregory, 2002; Chris Pires et al., 2006; Defoort et al., 2019). It is noteworthy that the species studied here have different forms of life (herbaceous vs. shrubby) and occur in different environments/biomes and consequently occupy different niches (Figure 1 and Supplementary Table S1). Nonetheless, only a few Cleomaceae species have been characterized in terms of their chromosome number (Ammal, 1933; Subramania and Susheela, 1988; Inda et al., 2008) and/or genome size (Omondi et al., 2017), and further efforts are required to understand the karyotype evolution in these taxa.

The events of gene/genome duplication observed for members of Cleomaceae contributed to several aspects of the evolution of C₄ photosynthesis in G. gynandra (Huang et al., 2021). Accordingly, this species has retained the duplicates of alanine aminotransferase and glutamine oxoglutarate aminotransferase. Besides almost all known vein-development-related paralogous genes derived from the genome duplication, this event also facilitated the evolution of C₄ enzyme genes and their recruitment into the C₄ pathway. In addition, several genes encoding photosystem I proteins were derived from the genome duplication (for more details see Huang et al., 2021). Furthermore, there are some examples of the duplication of gene copy number and its outcomes in the evolution of C₄ metabolism, such as the phosphoenolpyruvate carboxylase kinase gene families and the carbonic anhydrase genes (Wang et al., 2009). However, the number of enzymes required for C₄ metabolism is limited. In addition, the duplicated genes and other genomic regions often produce different effects (Lynch and Conery, 2003) and do not always result in functional innovation (Nielsen et al., 2005). In this context, some studies have proposed the involvement of parallel evolution in changes of gene expression and amino acid sequences (Christin et al., 2012a,b; Williams et al., 2012). Therefore, the cis-elements that direct cell specificity in C₄ leaves are present in C₃ orthologous genes recruited to the C₄ pathway, probably facilitating parallel evolution (Williams et al., 2012).

The increase and/or decrease in DNA content can promote phenotypic changes in the biochemical, anatomical, and physiological aspects of photosynthesis (Warner and Edwards, 1993; Vyas et al., 2007; Oa et al., 2011; Roddy et al., 2020). The importance of species relationships in trait-based studies is well appreciated in the field of comparative biology (Guy et al., 2019). Accordingly, G. gynandra, the species with C₄ photosynthesis, possesses a larger genome (Figure 2A). Nonetheless, a critical component in the study of adaptations is the identification of evolutionary correlations between phenotypic characteristics in phylogenetically close species (Adams, 2014). In this context, phylogenetic regression analysis provides a flexible analytical tool for assessing the degree of evolutionary association between variables while accounting for phylogeny (Adams, 2014).

The species with intermediate genome size (2C = 1.30 pg) (Figure 2A) had shorter ID and BSCW values, similar to the C₄ species G. gynandra (Figure 3 and Supplementary Figure S8). The lowest ID observed in G. gynandra can be attributed to the decline in the mesophyll (M) cell number, as verified, for instance, by Marshall et al. (2007) and Reeves et al. (2018), and size, which are considered essential for optimal C₄ function, since fewer M cells reduce the mean diffusion distance for C₄ metabolites (Hattersley, 1984; Khoshravesh et al., 2020). In addition, the enlargement of BSCs (equivalent to that in C₄ plants) is considered essential for C₄ function; as such, larger BSCs allow for a greater organelle volume and create larger vacuoles, serving as a barrier to facilitate CO₂ trapping in the sheath tissues (Von Caemmerer and Furbank, 2003; Khoshravesh et al., 2020). Consistent with previous reports, our results (Figure 2A) indicated that increase in the genome size affected leaf morphological and anatomical traits, in addition to physiological function (Roddy et al., 2020). Although the mechanism underlying the effect of genome size on cellular morphology remains unknown (Roddy et al., 2020), its effects on the vascular transport network, vessel size, and density have already been demonstrated (Maherali et al., 2009; Hao et al., 2013; De Baerdemaeker et al., 2018; Simonin and Roddy, 2018). In agreement, it has been suggested that genome size is a better predictor of guard cell size and stomatal density (Simonin and Roddy, 2018). Therefore, genome size can act as a first-order restriction on carbon gain, which directly affects the upper limit of allocation for growth, reproduction, and defense (Roddy et al., 2020).

Based on this premise, genome size may affect the tolerance of environmental changes related to latitude, altitude, temperature, precipitation, salinity, and desiccation (Grime and Mowforth, 1982; Otto and Whitton, 2000; Díez et al., 2013; Wang et al., 2018; Silva Artur et al., 2019; Zhang et al., 2019; Meyerson et al., 2020). Therefore, in habitats that can support high productivity and primary metabolic rates (e.g., Atlantic Forest and Amazon rainforest), species with smaller genomes (e.g., T. hassleriana, C. paludosa, and T. microcarpa; Supplementary Table S1) predominate, as they can maintain higher metabolic rates and rapidly adjust their physiology to match the environmental conditions (Simonin and Roddy, 2018; Roddy et al., 2020). It is thought to occur, since the genome size restricts the minimum cell size and the maximum cell compaction densities (Simova and Herben, 2012; Roddy et al., 2020). The increased cell size represents a direct physical constraint on the number of cells that can occupy a given space and, as a result, on the distance between cell types and tissues. Furthermore, reductions in cell size, necessary to pack more veins and stomata into leaves, effectively bring actual primary productivity closer to its maximum potential (Simonin and Roddy, 2018). Thus, leaves with many small stomata and a high density of veins can maintain higher rates of gas exchange than leaves with fewer, larger stomata and larger, less numerous, veins (De Boer et al., 2012). Accordingly, variation in cell size can drive large changes in potential carbon gain (Simonin and Roddy, 2018). Regarding the effects of stress, environmental and physiological factors influence the final sizes of somatic eukaryotic cells. In this sense, the minimum size of meristematic cells and the rate of their production are strongly constrained by nuclear volume, more commonly measured as genome size (Simonin and Roddy, 2018). Meanwhile, arid habitats (e.g., Cerrado, Brazilian savanna, and Caatinga, Brazilian semi-arid region; Figure 1) are characterized by low productivity and can support species with large genomes (Roddy et al., 2020), such as G. gynandra, T. spinosa, and T. longicarpa, since a larger genome may be associated with greater genetic diversity (heterozygosity).

Although we hypothesize that genome size can predict the metabolic rate, its effects may be probably more nuanced, since numerous studies that have correlated the genome size with ecological factors have often produced conflicting results, and this correlation may change according to the group analyzed. Notably, studies addressing the diversity of genome size among closely related species and the association of genome size with phenotypes and ecological factors are scarce. In this light, the present study is the first to address variations in genome size within Cleomaceae, and further research is warranted to draw conclusions.

C₄ photosynthesis arose from C₃ photosynthesis through a series of events/phases (Sage, 2004, 2021; Gowik and Westhoff, 2011). Combined physiological, molecular, and morphological changes played crucial roles in the evolution of the C₄ photosynthetic mechanism. The Cleomaceae accessions studied here were collected from warm and humid regions (Figure 1). For instance, C. paludosa, T. microcarpa, and most T. hassleriana accessions were collected from riverbanks in abandoned fields, except T. hassleriana (Joinville, SC), which was on an access road. The other species studied were collected from abandoned fields (T. longicarpa or G. gynandra). Similar to the species collected in the present study, some species with intermediate metabolism (e.g., Flaveria linearis; Holaday et al., 1984) grow in comparable environments and at disturbed sites, such as abandoned fields and roads. Thus, taking the steps essential for the development of the C₃–C₄ and C₄ photosynthetic mechanisms into account (see Sage, 2004; Gowik and Westhoff, 2011) and given that the studied Cleomaceae species occur in favorable environments (Figure 1) for the development of the C₃–C₄ photosynthetic mechanism, our results suggest that Tarenaya species as well as Cleoserrata may be pre-conditioned to evolve the characteristics associated with carbon concentrating mechanisms. Similar trends have been observed in a previous study on Cleome foliosa and Cleome africana (Marshall et al., 2007), which exhibited increased VD and enlarged BSCs. These alterations incorporate some of the most important changes required for C₄ photosynthesis.

Some Cleomaceae species exhibit phenotypic plasticity in leaf development and cell biology (Table 2; Marshall et al., 2007; Reeves et al., 2018), and some species have already undergone gene duplication events (Schranz and Mitchell-olds, 2006; Bergh et al., 2014). Furthermore, accessions of the same species (T. hassleriana and T. longicarpa) used in this study have a great genetic diversity, which can be observed by the different physiological, biochemical and anatomical data. This intraspecific diversity has already been observed for nine accessions of G. gynandra (Reeves et al., 2018) and can also be observed when comparing our accessions of T. hassleriana, T. spinosa, and G. gynandra with those from other continents (Table 2). However, modern phylogenetic methods allow us to analyze much more than just phylogeny as a statistical control. Specifically, they allow us to appreciate the evolutionary history of species, to better understand the patterns of biological diversity, to trace character traits over the evolutionary time, and to draw inferences regarding the evolution of these traits (Adams, 2014). As such, feature-based studies offer an important tool to better understand the ecological drivers of biological diversity. First, as observed in the present study, genome size was significantly correlated to variables essential for the C₄ photosynthetic mechanism, such as the δ¹³C value, WUEi, A_gross, A_N, malate, aspartate, starch, proteins, VD and size of palisade parenchyma (Table 1). As demonstrated for Grasses (Wang et al., 2009) and for G. gynandra (Huang et al., 2021) and reported in review papers (e.g., Monson, 2003; Sage, 2004, 2021; Gowik and Westhoff, 2011), the gene duplication or whole-genome duplication is thought to facilitate the evolution of C₄ photosynthesis from C₃ photosynthesis (e.g., Wang et al., 2009; Bianconi et al., 2018; Tronconi et al., 2020; Huang et al., 2021). The gene duplication or whole-genome duplication creates multiple copies of the gene, allowing for modification of the copies without the original function of the transcribed protein, so to reduce selective constraint and to acquire beneficial morphological or biochemical modifications (neofunctionalization) (Monson, 2003; Huang et al., 2021). It should be noted that gene duplications occur more frequently, since it is through sexual recombination, and thus are more likely to accumulate in short-lived annuals and perennials (Sage, 2004).

Furthermore, the second step in the development of the C₄ mechanism is characterized by increased VD (Sage, 2004; Gowik and Westhoff, 2011). The vascular system integrates with the photosynthetic tissues, and provides a flexible mechanical framework (Nelson and Dengler, 1997; Mckown and Dengler, 2010). In this way, in C₄ plants, the vascular system takes on an additional functional role, particularly C₄ biochemical cycling (Nelson and Dengler, 1997), due the ability to rapidly assimilate and reduce CO₂, and survive conditions that promote the loss of carbon and energy to photorespiratory processes (Mckown and Dengler, 2010). In the present study, 15 accessions exhibited higher VD than the C₄ species G. gynandra; four accessions exhibited values similar to G. gynandra (as verified for two Flaveria species C₃ and C₄, which presented similar VD values – Table 2); and only two accessions exhibited values lower than G. gynandra (Figure 3G). Of note, the VD in G. gynandra was equivalent to that in other accessions of the species (6–10 mm.mm^–2) (Marshall et al., 2007; Reeves et al., 2018). However, in the remaining C₃ accessions studied here, VD was higher than the previously reported values in the C₃ species C. violaceae, C. isomeris, C. hirta, and C. africana as well as the C₃-C₄ intermediate species Coalisina paradoxa (Marshall et al., 2007; Table 1). These values are too higher than those observed for Brassicaceae, Arabidopsis (C₃) and Moricandia (C₃ and C₃–C₄), and Asteraceae, Flaveria (C₃ and C₄) (Table 2).

Further, differences in VD observed in the present study may also be related to increase in the mechanical integrity of leaves, which may be beneficial in windy habitats or may improve water supply to leaves in dry and hot biotopes (Sage, 2004). The minor vein density, for example, is a key determinant of leaf hydraulic capacity and photosynthetic rates, and there would be strong benefit, independently of leaf size, in an ability to vary across a wide range of environments (Sack et al., 2012). In this regard, the characteristics that may initially be an adaptation to the environment, such as the combination of shorter ID, higher VD, and larger BSC, may predispose a species to developing the C₄ photosynthetic mechanism, as evidenced in grasses (Christin et al., 2013). Likewise, in Brassicaceae plants, VD has never been reported at the level of the C₄ species, indicating that this anatomical attribute may be one of the major constraints to the evolution of C₄ traits in this family (Schlüter et al., 2017). Therefore, increase in VD and enlargement of BSCs are the key alterations of foliar anatomy occurring in the C₃ context, preceding the emergence of the C₄ syndrome (Christin et al., 2013).

The third step in the development of the C₄ mechanism is marked by increase in the number of organelles in BSCs, resulting in cell enlargement. Although the activation of BSCs may be a secondary effect of increase in VD (Gowik and Westhoff, 2011), some accessions such as T. hassleriana (Canaã, MG, and Canoinhas, SC), T. siliculifera, and T. longicarpa (Manaus, AM) (Supplementary Figure S8) showed BSC size comparable to G. gynandra. The fourth step in the evolution of the C₄ mechanism is the appearance of glycine shuttling, a process in which the photorespiratory intermediate glycine is shuttled from the mesophyll tissue to the BSCs, where it is metabolized to CO₂ and serine by the action of glycine decarboxylase (Sage, 2004). Plants exhibiting this trait are considered C₃-C₄ intermediates because they present intermediate characteristics between the C₃ and C₄ mechanisms.

According to a previous carbon isotope discrimination analysis, T. siliculifera may present a modified photosynthetic mechanism (Voznesenskaya et al., 2007). However, in the present study, the T. siliculifera accession exhibited a δ¹³C value equivalent to the C₃ species, although anatomical and metabolic data indicated that it may be in transition from the C₃ to C₃–C₄ metabolism. Together, these results indicate that the studied T. siliculifera accession warrants further investigation, at both physiological and molecular levels, to define its photosynthetic metabolism.