INTRODUCTION

Evolutionary convergence refers to phenotypic similarities that have arisen in independently evolved species or populations. Convergent evolution is thought to be an adaptive process, where organisms respond similarly to the same environmental challenges. Because it comprises naturally occurring evolutionary replicates, convergent evolution serves as a valuable framework for understanding the flexibility or constraints governing the process of adaptation (Futuyma, 2010;Kitano et al., 2022;Losos, 2011;Orr, 2005;Sackton & Clark, 2019). An open question in biology is whether observed phenotypic convergence is often associated with convergence at the genetic level, and, if so, under what circumstances and to what extent these associations exist (Sackton & Clark, 2019).

In the context of agriculture, abscission, the process by which plants shed organs at a given developmental stage, is a trait that has been consistently altered in response to various environmental settings. The purpose of abscission is to get rid of unnecessary organs, like senescent leaves, or to disperse seeds or fruits. Natural seed shattering in wild plants is an abscission process of great ecological importance for the dispersal of mature seeds and the propagation of future generations (Dong & Wang, 2015). However, this trait is not favored in many cultivated plants, particularly cereals, and has been selected against to facilitate harvesting and to avoid yield loss (Fuller, 2007;Fuller & Allaby, 2009). In contrast, weeds in agricultural fields often leverage the shattering trait to improve their survival and maximize their dispersal (Maity et al., 2021;Vigueira et al., 2013). Thus, abscission has constantly been a target of selection in different groups within the agricultural setting: a reduced shattering during domestication and an elevated shattering level during weed evolution.

Rice (Oryza sativa) is an excellent system to examine the convergence of abscission, as seed shattering has been repeatedly modified throughout its evolutionary history. Stemming from two separate but entangled domestication events (Choi et al., 2017;Civ a n et al., 2015;Jing et al., 2023), five major rice cultivar subgroups, aus (also referred to as circum-Aus), indica, temperate japonica, tropical japonica, and aromatic (also referred to as circum-Basmati) are recognized today (Garris et al., 2005;Wang et al., 2018), and they are all characterized by reduced seed shattering levels compared with the wild ancestor of rice (Oryza rufipogon/nivara) (Thurber et al., 2010). Weedy rice (Oryza spp.) is a related weed of cultivated rice and aggressively outcompetes cultivated rice in crop fields. Weedy rice has evolved multiple times independently, most often from crop ancestors through a process of de-domestication (Cao et al., 2006;Hoyos et al., 2020;Imaizumi et al., 2021;Li et al., 2017Li et al., , 2022;;Qiu et al., 2020;Reagon et al., 2010), each time gaining a competitive advantage through re-acquisition of high seed shattering (Huang et al., 2017;Li et al., 2022;Thurber et al., 2010). The convergence of reduced shattering among different rice cultivars and of regained shattering in various weedy rice populations can provide a robust framework for examining convergent evolution of the abscission process.

We have previously investigated the convergence of seed shattering in rice at the morphological level through examination of the abscission zone (AZ), a specialized cell layer responsible for seed detachment (Li et al., 2024). In rice, the AZ is positioned in the rachilla, above rudimental glumes and below sterile lemma. The AZ-type cells are typically small, cytoplasmically dense, and lack lignin compared with their surrounding cells (Li et al., 2024). We have observed morphological convergence among independently evolved weedy rice groups, with all weedy rice displaying a complete and obvious AZ (Li et al., 2024). In contrast, distinct AZ region morphologies with different levels of disruption are prevalent in different cultivated rice populations (Li et al., 2024). In particular, the number and distribution of AZ-type cells varies greatly and significantly impacts seed shattering level (Li et al., 2024). However, the extent to which there is convergence in the genetic mechanisms governing AZ-type cell identity, cell number, cell distribution, and, ultimately, shattering levels in different rice lineages is unknown.

Here, we conducted a comparative tissue-specific RNA-seq analysis at three developmental stages for two cultivated and two weedy rice genotypes. We selected one cultivated rice genotype from each of two groups, aus and temperate japonica, that are known to stem from different domestication events and differ drastically in AZ morphology (Li et al., 2024); we also selected a BlackHull Awned (BHA) weedy rice genotype, which is known to descend from aus cultivated ancestors, and a Spanish Weedy Rice (SWR) genotype, which is known to descend from temperate japonica cultivated ancestors (Li et al., 2022;Reagon et al., 2010), that converge in both high seed shattering and in AZ morphology. Each crop-weed pair thus represents two independent evolutionary lineages, the aus lineage and the temperate japonica lineage, in which seed shattering and the AZ have been modified. We also examined AZ morphologies across our sampled developmental stages and their correspondence to patterns of gene expression. We found that the AZ-type cell characteristics of lack of lignin, small size, and thin cell walls were consistently throughout development. Transcription factors related to the repression of lignin and secondary cell wall deposition, along with certain known shattering genes involved in AZ formation, likely played a role in maintaining AZ region identity in both lineages. However, the majority of genes correlating with the morphological differences in the AZ between weedy and cultivated rice were unique to each lineage. More genes exhibited shared AZ-region exclusive expression pattern in domesticated and weedy genotypes in the aus lineage than in the temperate japonica lineage; and the modules exhibiting an AZ-region exclusive expression pattern in the weed but no expression in the crop were only found in the temperate japonica lineage. These findings align with the higher level of AZ morphological similarity between BHA weeds and aus cultivars compared to SWR weeds and temperate japonica cultivars. Our study reveals the degree of genetic convergence and divergence shaping AZ morphology, enhancing our understanding of the evolutionary malleability of the abscission process under selection.

RESULTS

Alterations in seed shattering levels and AZ morphology across development in weedy and cultivated rice

We have previously established that seed shattering levels in rice typically stabilize at 30 days after the panicle has emerged halfway from the boot (heading) (Thurber et al., 2010(Thurber et al., , 2011)). Using Breaking tensile strength (BTS), a widely employed measure for assessing the force required to separate florets or seeds from plants (Turra et al., 2023), we examined if each of our weedy and cultivated genotypes of the aus and temperate japonica lineages followed this typical shattering trajectory throughout development.

We measured BTS at heading/flowering time (FT), 1 week after heading (LFT), and 30 days after heading (LFT30). We attempted to measure BTS 1 week prior to heading (BFT) but failed due to the fragility of floral tissue. Large variation was evident in all four genotypes during the two early developmental stages (Table S1) with no significant increase or decrease of BTS values from FT to LFT in the temperate japonica lineage while relatively clear decrease from FT to LFT in the aus lineage (Figure S1). At LFT30, BTS decreased for all genotypes with values consistent with their expected seed shattering abilities. At this point, the weedy BHA and SWR genotypes exhibited high shattering tendencies with median BTS values of 0 gf (Table S1; Figure S1); the aus cultivar genotype, which is known to belong to a relatively easy-threshing cultivated group, presented a median BTS value of 5.50 gf; and the temperate japonica cultivar genotype, which belongs to a generally low shattering cultivated variety, had an LFT30 median BTS value of 87.8 gf, distinctive from the other genotypes (Table S1; Figure S1).

To understand the relationship between shattering levels and AZ morphology throughout development, we also examined the histology of AZs in these four genotypes at three developmental stages: BFT, FT and LFT. Longitudinal sections on seeds collected at LFT30 cannot be conducted due to increased starch content and hardness of the seeds, and the seeds being prone to separation from the pedicel in high shattering rice. In the two high shattering weedy rice genotypes and in the relatively easy shattering aus cultivar, the AZ cell layer was visibly present throughout all three developmental stages (Figure 1d-l). The AZ-type cells distinguished themselves from neighboring cells with their distinctive dark purple color, indicating a lack of lignification, thin walls, and smaller sizes. A small gap existed between the AZ cell layer and middle vascular bundle in the aus cultivar (Figure 1g-i), serving as the major difference in AZ morphology between BHA and aus cultivar and corresponding to their shattering level disparities. In contrast, in the AZ region of the temperate japonica cultivar, no distinct AZ cell layer formed at any developmental stage and only a few scattered AZ-type cells were present (Figure 1a-c). Non-AZ-type cells within the AZ region exhibited a turquoise color, suggesting lignification. Additionally, these non-AZ-type cells featured thick rather than thin cell walls compared to AZ-type cells, underscoring their cellular distinctions.

Effects of developmental stage, lineage, and tissue type on gene expression patterns

To identify the genes contributing to AZ development in these four rice genotypes, we conducted tissue-specific RNA-seq analyses. This involved the manual dissection of the AZ region, of the adjacent tissue toward floret side, and of the adjacent tissue toward pedicel side (Figure S2), at the developmental stages corresponding to our examination of AZ histology: BFT, FT and LFT. For each genotype, we collected three biological replicates of each tissue at each developmental stage. The RNA-seq reads were mapped to an indica reference genome, R498 (Du et al., 2017) with an average overall mapping rate greater than 85%. Principal component analyses (PCA) indicated that developmental stage played the most significant role in distinguishing the expression profiles among samples (Figure 2a; Table S2). Notably, the early developmental stage was clearly separated from the two later developmental stages along principal component 1 (PC1), explaining 62.3% of the variation. The genetic background specific to the evolutionary lineage was the second factor shaping the clustering of the expression patterns (Figure 2b). Along PC2, which accounts for 30.7% of the variation, a clear separation between the aus and temperate japonica lineages was evident. However, the distinction between weedy and cultivated genotypes within the same lineage was not as pronounced as the separation by lineages, despite their shattering differences. Tissue specific expression patterns were observed, with replicates from the same tissues clustering together, particularly in the latter developmental stages (Figure 2c). These patterns also exhibited variation between different lineages. During the early developmental stage, the separation among tissues was more apparent in the aus lineage compared to the temperate japonica lineage (Figure 2c).

AZ-region associated modules and their relation to AZ morphology in weedy and cultivated rice from each lineage

Modules represent groups of genes that share similar expression profiles, often indicative of functional relatedness and co-regulation (Saelens et al., 2018). We employed weighted gene co-expression network analysis (WGCNA) (Langfelder & Horvath, 2008) to identify modules exhibiting exclusive expression in the AZ region, reasoning that genes within these modules are more likely to affect seed shattering abilities. We performed module detection by jointly analyzing the expression profiles of the cultivar and its corresponding weedy counterpart. In the WGCNA analyses for the temperate japonica cultivar and SWR, we used gene counts obtained from alignment to a japonica reference genome (Kawahara et al., 2013) to mitigate potential issues arising from genetic disparities between japonica and indica, which could lead to missing sequence or misalignment (Li et al., 2023).

We first clustered samples from each lineage using their expression profiles. The results mirrored the patterns observed in PCA, clearly distinguishing samples collected in the early developmental stage from those in the two latter stages (Figures S3 and S4). During both the BFT and LFT stages, AZ samples from the aus cultivar and BHA genotypes consistently grouped into a single sub-cluster (Figure S3), suggesting the presence of a similar expression pattern, aligning with their similar AZ histology (Figure 1g-l). In contrast, the AZ samples in the temperate japonica lineage usually clustered with the pedicel or floret samples of the same genotype at each given developmental stage, indicating distinct expression profiles in the AZ regions of the temperate japonica cultivar and SWR (Figure S4). The notable differences in the AZ histology between these two cultivated and weedy genotypes could account for this disparity (Figure 1a-f).

Modules explaining AZ similarities and differences between weedy and cultivated rice

We identified modules displaying an AZ-region exclusive expression pattern in both weedy and cultivated genotypes from the same lineage. These modules are more likely to contribute to the AZ region identity or to AZ morphological features shared between genotypes within a lineage. We correlated the eigengene expression, a representative summary measure of the gene expression profiles within a module, with tissue classification by assigning 1 to samples from the AZ region in both cultivated and weedy genotypes and 0 to all other tissues. Modules with a correlation co-efficiency greater than 0.7 and a P-value less than 0.05 were defined as AZ-region associated. A greater number of AZ-region associated modules and notably more genes were observed within these modules in the aus lineage (six modules; Figure 3a,b,e,f,i,j) compared to the temperate japonica lineage (four modules; Figure 3c,g, h,k) (Tables S3-S6; Figure 3). These findings are consistent with the greater similarity in AZ morphology observed between BHA and aus cultivar than between SWR and temperate japonica cultivar (Figure 1).

We also identified modules exhibiting an AZ-region exclusive expression pattern specific to weedy rice but not to cultivated rice, through a separate correlation analysis assigning 1 to samples from the AZ region in weedy genotypes and 0 to all the other tissues in cultivated and weedy genotypes (Table S3). Genes in these modules may contribute to the development, maintenance, or competence of a complete AZ in weedy rice, and to weedy rice's greater shattering ability; the altered expression patterns of these genes in cultivated rice may also be responsible for the Ó 2025 Society for Experimental Biology and John Wiley & Sons Ltd., The Plant Journal, (2025), 121, e70083 morphological variation observed in the AZ region of cultivars. We identified weed-specific AZ-region associated modules with two types of expression patterns (Figure 3). One type comprised genes showing a distinctive AZ-region exclusive expression pattern solely in the weed, with no expression in the crop (Table S7; Figure 3d,l). The other type of module displayed the same AZ-region exclusive expression pattern in both the weed and crop, but with different expression levels between the weed and crop in the AZ-region (Figure 3i,h); these modules were thus also detected in the correlation analyses involving crop and weedy genotypes concurrently. In the temperate japonica lineage, three modules, one for each developmental stage, exhibited significant correlations with the weedy AZ region, with two being the first type module (Figure 3d,l) and one being the second type (Figure 3h). In the aus lineage, only one module detected in the later developmental stage, was significantly correlated with the weedy AZ region, exhibiting the expression pattern of the second type (Figure 3i). These results suggest that the same genes or genetic networks are involved in AZ development in BHA and aus cultivars, but it is their varying expression levels that lead to the slight differences in AZ morphology and the differences in seed shattering (Figure 1; Figure S1).

Conversely, the restoration of a complete AZ in SWR seems to entail activation of genetic networks at multiple developmental stages that are largely not active in cultivated temperate japonica.

The occurrence of genes related to lignin and cell wall modification within AZ-region associated modules We examined each AZ-region associated module for enriched gene ontology (GO) terms, including modules associated in both cultivated and weedy genotypes and those associated with the AZ region only in weedy genotypes. Many modules did not have any significant GOterm enrichment and the number of enriched terms, when present, was usually modest, and primarily associated with various metabolic processes which did not provide meaningful insights into AZ-specific biological processes (Tables S8 and S9; Figures S5 and S6). Given the observed characteristics of AZ-type cells, we then identified genes in each AZ-region associated module with GO terms related to lignin biosynthesis and cell wall modification. In total we found 8 and 24 genes associated with lignin metabolic process (GO:0009808) and cell wall organization and biogenesis (GO:0071554, GO:0042546 or GO:0009827), respectively (Table S10). In the aus lineage, genes related to lignin modification were present in four out of six AZ-region associated modules for both genotypes, including the module more highly expressed in the BHA weedy genotype AZ region (Table S10). Surprisingly, in the temperate japonica lineage, only two AZ-region associated modules for both genotypes identified in the BFT and FT stages contained lignin modification related genes, and none of the weedspecific AZ-region associated modules, despite the clear difference in lignified cells in the AZ region between cultivated and weedy rice (Table S10; Figure 1). Notably, OsWRKY102 a transcriptional repressor of secondary cell wall formation (Miyamoto et al., 2020), was identified in the AZ-region associated modules in BHA and aus cultivar genotypes during both the FT and LFT stages, and in SWR and temperate japonica cultivar genotypes during the FT stage. This gene was also present in the module associated with the AZ region only in BHA genotype (Table S10). A loss-of-function mutation in this gene results in increased lignin levels in rice (Miyamoto et al., 2020), thus the exclusive expression of OsWRKY102 in the AZ region suggests its necessity in repressing lignin deposition. Its expression in the AZ region of the temperate japonica cultivar genotype, which does not form a distinct AZ cell layer but contains only a few AZ-type cells, suggests this gene is important for maintaining overall AZ region identity.

Cell wall modifications were regulated in the AZ-region throughout all three developmental stages in both lineages. Generally, more genes involved in cell wall modification were found in AZ-region associated modules during the BFT stage compared to the later developmental stages, particularly in the BHA and aus cultivar genotypes (Table S10). During the BFT stage, the modification of pectin, regulating intercellular adhesion, appeared to be important for AZ development. Pectin Methylesterase Inhibitor (PMEI ) genes, OsPMEI 17 and OsPME136, appeared as hub genes, that is, genes with the highest top 10% connectivity within AZ-region associated modules, detected for BHA and aus cultivar genotypes and for SWR and temperate japonica cultivar genotypes during this stage (Tables S6 and S10). Several genes related to cell wall modification were also found within the weed-exclusive AZ-region associated modules (Table S10). For instance, LOC_Os02g32110 in the SWR AZ-region associated module is a homolog of AtIRX10 in Arabidopsis thaliana, which regulates the synthesis of xylan, a hemicellulose in the cell wall, and exhibits a conserved function over a substantial evolutionary distance (Jensen et al., 2014). Another example was a member of the xyloglucan endotransglucosylase/hydrolase (XTH ) gene family, XTH3, in the BHA AZ-region associated module. Genes in the XTH family are involved in the biosynthesis of xyloglucan, another type of hemicellulose in cell wall (Yokoyama et al., 2004).

Genes overlapping in the AZ-region associated modules among different rice lineages

There was limited overlap of genes within AZ-region associated modules between the aus and temperate japonica lineages (Tables S4, S5 and S11). The highest number of overlapping genes, 14, was observed at the BFT stage (Figure 4; Table S11). Notably, an APETALA2-like transcription factor (LOC_Os03g60430) exhibiting a high degree of similarity to known seed shattering genes, SHAT1 and SNB, which are involved in AZ formation (Wang et al., 2019), was detected in AZ-region associated modules in both lineages (Table S11). This gene also served as a hub gene in the AZ-region associated module for BHA and aus cultivar genotypes during this BFT stage (Table S6). The pectin regulation gene OsPMEI36 (LOC_Os10g36500) also overlapped between modules of the two lineages during the BFT stage (Table S11). These results suggest that at this stage, determining AZ region identity and modifying cell wall composition are likely shared processes between these two lineages.

At the FT stage, three genes stood out among those overlapping between the AZ-region associated modules in both lineages: OsWRKY102 (LOC_Os01g08710), GA2-oxidase7 (GA2ox7) (LOC_Os01g11150) and a NAC transcription factor (OsXND1-like) (LOC_Os04g35660) (Table S11). OsWRKY102 functions as lignin biosynthesis repressor (Miyamoto et al., 2020). GA2ox7 codes for an enzyme involved in regulating gibberellin (GA), and modulating GA content can impact lignin deposition and ultimately alter seed shattering (Wu et al., 2023). OsXND1-like is homologous to A. thaliana Xylem NAC domain 1 (XND1), which negatively regulates secondary cell wall synthesis and promotes programmed cell death in xylem (Zhao et al., 2008). Ectopic expression of XND1 in fiber cells failed secondary cell wall deposition and dramatically reduced the expression of secondary wall biosynthesis genes (Zhao et al., 2008). This gene was hub gene in the AZ-associated modules detected in both the BHA and aus cultivar genotypes as well as SWR and temperate japonica cultivar genotypes at this developmental stage (Table S6).

Importantly, this gene ranked as the second most interconnected gene within the AZ-associated module detected in the BHA and aus cultivar genotypes (Table S6). These findings indicate that the primary processes shared among these four genotypes during the FT developmental stage included the repression of lignin deposition and the inhibition of secondary cell wall formation in the AZ region.

During the LFT stage, only two genes, specifically a protein kinase and a membrane protein, exhibited overlap (Table S11). Compared to the other developmental stages, there was a reduced number of shared genes among these genotypes at this stage.

Genes overlapping across development in AZ-region associated modules

There was no overlap observed among AZ-region associated modules across all three developmental stages within each lineage, regardless of whether they were exclusively associated with the weedy AZ region or simultaneously associated with both weedy and cultivated AZ regions (Figure 5). More genes were shared between AZ-region associated modules from two different developmental stages in the aus lineage compared to the temperate japonica lineage (Table S12).

The most substantial overlap (16 genes) was observed between the FT (orange) and LFT (magenta) AZ-region associated modules in the aus cultivar and BHA genotypes (Figure 5), in line with the results of the PCA where samples from these two stages clustered together (Figure 2a). Interestingly, the genes that overlapped between these two AZ-region associated modules included four out of six genes that also overlapped between AZ-region associated modules from different lineages during the FT developmental stage (Table 12). This included one coding for a MATE efflux family protein and three previously mentioned genes, OsWRY102, GA2ox7 and OsXND1-like, suggesting a consistent requirement for repressing lignin deposition and secondary cell wall deposition and the potential involvement of the same genetic networks in regulating these processes in the AZ in both developmental stages in BHA and the aus cultivar. However, only one gene, ubiquitin domain-containing protein 2, overlapped between FT and LFT stages in the AZ-region associated modules identified in the temperate japonica lineage (Table 12).

Previously known seed shattering genes exhibiting tissuespecific expression patterns

Several genes associated with seed shattering levels in rice have been previously identified via crosses or mutagenesis (Turra et al., 2023). We examined modules containing these previously identified genes to determine how they related to the shattering differences observed among crops and weeds in each of our lineages. We particularly focused on the genes that have been functionally verified using transformation and are associated with AZ region development. Eight such genes, which include SH5 (Yoon et al., 2014), SNB (Jiang et al., 2019), SHAT1 (Zhou et al., 2012), OsCPL1 (Ji et al., 2010), OSH15 (Yoon et al., 2017), SLR1 (Wu et al., 2023), OsGRF4 (Sun et al., 2016) and NPC1 (Cao et al., 2016), were mainly identified through mutagenesis, and their roles in loss of shattering during domestication remaining unproven. Three such genes -SH4 (Li et al., 2006), qSH1 (Konishi et al., 2006), and OsSh1 (Li et al., 2020)were primarily identified through crosses, and natural variation within these genes is associated with reduced seed shattering during domestication for some cultivated rice lineages. Cultivated alleles of SH4 occur both in aus and temperate japonica cultivars, but no role for SH4 has been yet demonstrated for seed shattering in weedy rice, which often also carries crop alleles (Thurber et al., 2010;Turra et al., 2023). For qSH1 and OsSh1, there is some degree of evidence of involvement during domestication or improvement in the temperate japonica lineage, however, their role in reduced seed shattering in aus cultivars remains undemonstrated (Li et al., 2020). We additionally verified that weedy rice samples in this study consistently carried the same alleles at the inferred domestication-related mutations in SH4, qSH1, and OsSh1 as the cultivars within the same lineage (Table S13). We found that most of these 11 known shattering candidate genes occurred in modules that were expressed in all three developmental stages, with the exception of SH4 (Figure S7). The eigengenes of the majority of these modules containing candidate genes exhibited a tissue-specific expression pattern, with non-uniform expression across all three tissues (Figure S7), suggesting that the specific location of gene expression was crucial for fulfilling their functions. The interactions among these genes were also reflected in their memberships within the same module. For example, OsSh1, SNB and qSH1 were included in the same module detected in BHA and aus cultivar during BFT (Figure S7). While OsSh1 and qSH1 were in the same module detected in SWR and temperate japonica cultivar during BFT; OsSh1 and SNB were in the same module detected in SWR and temperate japonica cultivar during LFT (Figure S7). These results suggested that although these genes work together in both lineages, the nature of interactions might differ. Because most of these genes have been previously implicated in AZ formation in early floral development, we focused on the modules that contain these known shattering candidate genes at the BFT stage. qSH1 (Konishi et al., 2006) and OsSh1 (Ishikawa et al., 2022;Li et al., 2020) were detected in modules demonstrating expression in the AZ and floret tissues in the BHA, aus cultivar and SWR genotypes, while no expression in the AZ was observed in the temperate japonica cultivar (Figure 6a). A similar situation where expression in the AZ and pedicel tissues was evident in the rice varieties with a well-defined AZ morphology, BHA, aus cultivar, and SWR, while no AZ expression was discernible in the temperate japonica cultivar, was also observed for the modules that contained SH5 (Yoon et al., 2014) (Figure 6b) and NPC1 (Cao et al., 2016) (Figure 6c). The qSH1 and OsSh1 containing modules (Figure 6a) overlapped by 258 genes between the two lineages (2260 genes in the aus lineage while 735 genes in the temperate japonica lineage) and the NPC1 containing modules (Figure 6c) overlapped by 158 genes (1610 genes in the aus lineage while 619 genes in the temperate japonica lineage). These findings suggest that the genes within these modules might form certain genetic networks and be crucial for an AZ layer presence and shattering in both lineages.

We also found that certain AZ genes were part of modules with a consistent expression pattern between the crop and weed belonging to the same evolutionary lineage but exhibited shifted expression patterns in rice from different lineages. For example, SHAT1 (Zhou et al., 2012) and OSH15 (Yoon et al., 2017) were both members of a module exclusively expressed in the AZ in the BHA and aus cultivar genotypes, but were detected in a module expressed in both AZ and pedicel tissues in the SWR and temperate japonica cultivar genotypes (Figure 6d). Thirtyfive genes overlapped between these SHAT1 and OSH15 containing modules (Figure 6d) between the two lineages (476 genes in the aus lineage while 898 genes in the temperate japonica lineage). Conversely, SH4 (Li et al., 2006) occurred in a module displaying exclusive expression in the AZ in the SWR and temperate japonica cultivar genotypes, while the SH4 containing module in the BHA and aus cultivar genotypes did not exhibit a clear pattern, possibly showing higher expression in both AZ and floret tissues (Figure 6e). There were also situations where a gene belonging to a module demonstrating expression in the AZ and floret tissues in one lineage, displayed a relatively random expression pattern in the other lineage; SNB (Jiang et al., 2019) was one such example (Figure 6f). Few genes of these modules (Figure 6e,f) overlapped between the two lineages.

Genes exhibiting differential expression in the AZ between weedy and cultivated rice in the aus and temperate japonica lineages

To identify genes potentially responsible for variation in seed shattering between cultivated and weedy rice within each lineage, we performed differential expression analyses for each tissue type at every developmental stage. Hundreds of genes were differentially expressed (DE) between crops and weeds in every tissue and at every stage. We focused on DE genes within the AZ and found larger numbers between BHA and the aus cultivar -545, 712 and 642 at the BFT, FT and LFT developmental stages, respectively (Table S14)than between SWR and the temperate japonica cultivar -420, 312 and 368, respectively (Table S15). None of known shattering genes exhibited significant expression difference within the AZ between the crop and weed in any lineage. This difference in DE gene number, however, was reduced once we examined genes exclusively DE in the AZ but not in the other tissues; 132, 234, and 133 AZ-exclusive DE genes were identified between BHA and the aus cultivar, while 162, 70 and 152 were found between SWR and the temperate japonica cultivar at the BFT, FT and LFT developmental stages, respectively (Table S16). Among these AZ exclusive DE genes, 28, 30, and 26 genes overlapped at each stage between the weed-crop comparisons detected in two different lineages, although the directions of regulation of some of these genes were not consistent between comparisons (Table S17).

GO term enrichment for weedy-cultivated rice DE genes

Several related GO terms were consistently enriched for AZ-exclusive DE genes between BHA and the aus cultivar across all three developmental stages, such as "cell death", "programmed cell death", "apoptosis" and "response to stress" (Figure 7; Table S18). Numerous GO terms related to metabolic and cellular regulation were exclusively enriched during the FT stage (Figure 7; Table S18). In contrast, no enriched GO terms for the AZexclusive DE genes between SWR and the temperate japonica cultivar were shared across the three developmental stages. During the FT stage, only a single GO term, "response to biotic stimulus", displayed enrichment. During the BFT stage, multiple GO terms related to "photosynthesis" and one GO term related to "cell-wall" were enriched. During the LFT stage, the enriched GO terms encompassed "response to stress/stimulus", "DNA packaging" and "DNA conformation change" (Figure 7; Table S18). Cell death related GO terms were not enriched for AZ-exclusive DE genes between SWR and temperate japonica cultivar. The only shared enriched GO terms with the BHA-aus cultivar comparison were "response to stress" and "response to stimulus." This suggests that different regulatory processes are responsible for the differences in seed shattering and AZ morphology between crops and weeds in each rice lineage. Particularly, the enriched GO terms in the SWR and the temperate japonica cultivar comparison appeared to be more focused on responses while those in the BHA and the aus cultivar comparison were more related with cell identity (Figure 7).

Genes related to lignin modification differentially regulated between weedy and cultivated rice

Similar to the GO term analysis in AZ-region associated modules, we examined genes annotated with GO terms of lignin metabolic process (GO:0009808) and cell wall organization or biogenesis (GO:0071554) among the DE genes exclusively in the AZ region in each weedy-cultivated comparison (Table S18).

During the BFT stage, three laccase genes (OsLACs) exhibited downregulation in the AZ region of BHA compared to the aus cultivar (Table S19). Laccases, a type of oxidative enzyme, play a crucial role in determining the spatial and temporal polymerization of lignin (Emonet & Hay, 2022). The diminished activity of these genes may lead to less lignin being deposited in the AZ of BHA, thereby contributing to its increased shattering compared to aus cultivar. Notably, OsLAC15_2 consistently demonstrated reduced expression in the AZ of BHA compared to the aus cultivar across all three developmental stages, suggesting a persistent reduction in lignin deposition in the AZ region of BHA throughout development. In line with this lignin repression, OsWRKY102, a repressor of lignin deposition, was found to be highly expressed in the AZ of BHA compared to the aus cultivar during the FT stage (Table S19). In contrast, there were fewer DE genes in the AZ region between SWR and the temperate japonica cultivar associated with lignin (Table S19). One of these genes was OsWRKY36, another repressor for lignin and secondary cell wall deposition (Table S19). OsWRKY36 activates the transcription of SLR1 by directly binding to its promoter (Lan et al., 2020); SLR1, in turn, functions as a key repressor of GA signaling directly involved in seed shattering in rice through impacting lignin content in the AZ region (Wu et al., 2023). Interestingly, OsWRKY36, displayed lower expression in the AZ of SWR compared to the temperate japonica cultivar during the LFT developmental stage (Table S19), despite lower levels of lignin expected in SWR.

Pectin methylesterase (PME ), XTH genes, cellulose Synthaselike (CSL) genes, involved in biosynthesis of pectin, hemicellulose, and cellulose respectively, exhibited differential gene regulation between weedy and cultivated rice in both lineages (Table S19). The modification of these components can impact the deposition of secondary cell wall (Coomey et al., 2020). However, examining the precise impact of these genes on cell wall modifications can be challenging, given that the function and interactions of these genes have not be well-characterized in rice.

DISCUSSION

Seed shattering is one of the most consistently differentiated traits that distinguishes weedy from cultivated rice, with high seed shattering significantly increasing the competitiveness of weeds in rice fields (Vigueira et al., 2013). The repeated evolution of high shattering in de-domesticated weedy rice stemming from different cultivated ancestors provides an excellent framework for exploring mechanisms of convergent evolution. Our previous examinations of AZ histology suggest that AZ-type cells converge in their cellular characteristics, irrespective of the evolutionary lineage; in contrast, the arrangement of AZ-type cells, in particular the length of the AZ cell layer, is associated with different levels of seed shattering across different weedy and cultivated rice populations (Li et al., 2024). To date, our understanding of the loci responsible for the evolution of high seed shattering in weedy rice from different domesticated lineages remains limited. Multiple loci are known to contribute to the shattering trait, and most of de-domesticated weedy genotypes carry the same alleles of known seed shattering genes as those found in cultivated rice (Table S13) (Li et al., 2022(Li et al., , 2023)). By examining patterns of gene expression in the AZ region of related weedy-cultivated rice pairs we have been able to identify the genetic associations with AZ morphology and shattering levels in the aus and temperate japonica lineages, and the extent of genetic convergence. Our findings indicate constraints on the genetic mechanisms determining AZ-type cell identity, but flexibility in the genetic networks shaping AZ morphology, particularly concerning AZ-type cell distribution, underscoring heterogeneity in the levels of genetic convergence underlying the crucial agricultural trait of shattering. Genetic convergence underlying the maintenance of AZ identity across lineages AZ-type cells exhibit uniform cellular characteristics, featuring the absence of lignin, small size and thin walls, across our four examined rice genotypes, though the number and the distribution of AZ-type cells varied (Figure 1). For example, cells are few and distributed sparsely in the temperate japonica cultivar, in contrast to the other genotypes. These AZ-type cells characteristics appear very early and persist throughout floral development (Figure 1), emphasizing the necessity in maintaining AZ-type cell identity. Some genes related to lignin and secondary cell Significant Gene Ontology (GO) terms enriched in the DE genes in the AZ between BHA and the aus cultivar as well as SWR and the temperate japonica cultivar at the BFT, FT and LFT developmental stages. Red and blue represent smaller and larger FDR values, respectively. The size of the dots represents the number of genes in each GO term over the total number of DE genes in the AZ. AZ, abscission zone. walls, especially transcription factors, consistently display comparable AZ-exclusive or tissue-specific expression patterns in both weedy and cultivated genotypes across lineages. These genes are likely instrumental in determining the distinct AZ-type cell identity, as evidence of genetic convergence in both lineages.

One such example is a gene encoding an AP2-like transcription factor detected in the AZ-region module for weedy and cultivated genotypes in both lineage during the BFT stage (Table S11). This gene exhibits high similarity to SHAT1 and SNB, both involved in AZ formation, the latter of which can inhibit lignin deposition in the AZ (Jiang et al., 2019). Another example is OsWRKY102, a repressor for lignin biosynthesis, which not only appears in AZ-region associated modules from both lineages during the FT stage (Figure 4; Table S11), but also preserves this AZ-exclusive expression pattern in BHA and the aus cultivar into the LFT stage (Figure 5; Table S12). Additionally, a higher expression of this gene is observed in the AZ of BHA compared to the aus cultivar (Table S19), correlated with the greater number of AZ-type cells in BHA, an expected outcome of more repression of lignin deposition. These findings suggest that lignin deposition remains regulated consistently throughout AZ development. Nevertheless, this repression may undergo some changes, as evidenced by the observed decrease of expression of OsWRKY102 during the LFT stage compared to the FT stage in BHA and the aus cultivar. Repression of secondary cell wall deposition is another process consistently regulated across AZ development. One transcription factor, OsXND1-like, associated with repressing secondary cell wall deposition, occurred in the AZ-region associated modules in both lineages during the FT stage (Table S11) and in BHA and the aus cultivar at the LFT stage (Table S12). Several cell wall modification enzymes also exhibited AZexclusive expression pattern in both lineages (Table S10) supporting genetic convergence.

Certain known seed shattering genes interact, and the outcomes of their interactions include the alteration of lignin or secondary cell wall deposition. Although these genes did not exclusively express in the AZ region, they exhibited a tissue-specific expression pattern and possibly interacted with each other if they were in the same module. For example, OsSh1 functions downstream of qSH1 (Li et al., 2020) while qSH1 is directly regulated by SNB (Jiang et al., 2019). The direct impact of their interaction is the suppression of lignin deposition in the AZ region, and mutations in SNB lead to altered development of AZ (Jiang et al., 2019). These three genes were found in the same module, exhibiting higher expression in the AZ region and floret tissues in BHA and aus cultivar during the BFT stage (Figure S7); additionally, qSH1 and OsSh1 are in the same module in SWR and the temperate japonica cultivar during BFT, and OsSh1 and SNB are in the same module in these two genotypes during LFT (Figure S7). NPC1 can impact the secondary cell wall deposition through modulating silicon distribution and resulting in changes in the content of cellulose and hemicellulose in the cell wall and the thickness of cell wall (Cao et al., 2016). This gene is included in the module exhibiting higher expression in the pedicel and AZ region during the BFT stage in BHA, the aus cultivar and SWR (Figure 6c), and potentially interacts with SH5 in BHA and the aus cultivar (Figure 6b). Additionally, relatively large number of genes overlapped between the qSH1 and OsSh1-containing and NPC1-containing modules from the two lineages (Figure 6a,c).

Although indirect, these results suggest that a shared group of genes, particularly transcription factors, are involved in determining AZ region identity in different genotypes of rice, indicating possible genetic constraints to making AZ-type cells in rice.

Genetic divergence underlying AZ morphological variation between weedy and cultivated rice in different lineages

How the AZ morphology varies between cultivated and weedy rice is lineage-dependent, suggesting that distinct genetic mechanisms explain these weed-crop differences. The primary AZ morphological distinction between BHA and aus cultivars lies in the length of the AZ (Figure 1), indicating that the genetic basis for AZ region formation is likely similar, with differences possibly in allowing or repressing AZ formation close to the vascular bundle (Figure 1g-i). This aligns with our finding that all AZ-region associated modules were consistently correlated with the AZ region in both BHA and the aus cultivar (Figure 3). Additionally, the modules including previously known seed shattering genes largely exhibited similar tissue-specific expression patterns between these two genotypes (Figure S7). Thus, the same networks appear to be utilized in both genotypes to maintain AZ cell identity, but these genes might exhibit varied expression levels and result in the AZ morphological difference and seed shattering difference between them. One possible causal gene might be SHAT1, a master transcription factor determining AZ identity (Zhou et al., 2012). Although not significant, this gene exhibits higher expression in the AZ in BHA than in the aus cultivar. This gene is also in a seed-shattering QTL detected through a cross of the same BHA and aus cultivar genotypes from this study and this QTL overlaps with a selective sweep (Li et al., 2023). One single nucleotide mutation in the 3 0 UTR (untranslated region) of this gene exhibited association with the change of seed shattering levels (Li et al., 2023).

In contrast to weeds and crops in the aus lineage, substantial morphological distinctions are present in the AZ region between SWR and temperate japonica cultivars. The discontinuous AZ in temperate japonica cultivars (Figure 1) suggests that AZ-type cells are randomly scattered in the AZ region, lacking the well-defined AZ cell layer observed in aus cultivars. Consistent with this disrupted or underdeveloped AZ region in the temperate japonica cultivar, fewer AZ-region associated modules were shared by both genotypes in the temperate japonica lineage compared to the aus lineage. Additionally, modules with AZregion exclusive expression in weeds but no expression in crop were only detected in the temperate japonica lineage (Figure 3). The random and scattered distribution of AZ cells in the temperate japonica cultivar suggests inconsistent regulation of AZ cell formation in this genotype. For example, the module that included OsSh1 and qSH1 in the BFT stage, displayed comparable tissue-specific expression pattern in BHA, SWR and the aus cultivar, while an altered expression pattern in the temperate japonica cultivar (Figure 6a).

That different genes have been involved in the reevolution of easy shattering in different de-domestication events of weedy rice had also been previously suggested by the lack of overlapping shattering QTLs detected in mapping populations from weed 9 crop crosses involving BHA weedy rice, and SH weedy ricea group descended from indica domesticated varieties (Li et al., 2023;Qi et al., 2015;Thurber et al., 2013). Different sets of genes exhibiting differential expression between crop and weed within each lineage further illustrate that different genes involved in the reacquisition of high shattering in these two independently evolved weedy rice populations.

Overall, our results shed light on the genetic convergence underlying phenotypic convergence: consistent repression of lignin and secondary cell wall deposition throughout development is crucial for maintaining the identity of the AZ region, regardless of weedy or cultivated rice genotype; this phenotypic convergence is supported by the observation that some genes associated with these modifications exhibit exclusive expression in the AZ region across both weedy and cultivated genotypes. However, only a few of these genes overlap between the two lineages, suggesting that the specific genes responsible for modifying lignin and secondary cell walls may differ between them. This variation maybe attributed to distinct evolutionary histories.

Our conclusions are derived from contrasting two independently evolved weed-crop pairs. However, more than 10 independent origins of de-domesticated weedy rice have been documented so far (Hoyos et al., 2020;Huang et al., 2017;Li et al., 2022;Presotto et al., 2024;Qiu et al., 2020;Reagon et al., 2010;Zhu et al., 2024), the great majority of which involve modifications to shattering levels compared to their cultivated ancestors. This raises further intriguing questions about limits to the diversity of genetic mechanisms underlying a convergent trait. Additionally, the wild ancestor of cultivated rice, O. rufipogon, also has high levels of shattering and an AZ morphology similar to that of weedy rice (Thurber et al., 2011). However, the levels of convergence in the genetic mechanisms regulating AZ development and seed shattering in wild and weedy rice have not yet been examined. Future comparisons incorporating the diversity of weedy, wild, and cultivated rice groups will be valuable in revealing the extent of genetic flexibility in the control of this agronomically important trait.

EXPERIMENTAL PROCEDURES

Plant materials

Two cultivated rice lines, Aus 196 and AN-09-13 were chosen to represent the cultivated rice variety groups of aus and temperate japonica. Lines 10A and ASF-09-54, were chosen to represent the two weedy rice populations Black Hull Awned (BHA) and Spanish Weedy Rice (SWR), which evolved from aus and temperate japonica cultivars, respectively. Aus 196 is a breeding and research material developed in Bangladesh and acquired from the International Rice Research Institute; 10A is a weedy rice variety collected in Arkansas, U.S. and maintained by the United States Department of Agriculture; AN-09-13 and ASF-09-54 were obtained through field collections and sampling of germplasm repositories in the Iberian Peninsula (Li et al., 2022). At least two seeds of the same genotype were planted in one pot; each pot served as one replicate for each genotype; and three replicates were used for each genotype in the RNA-seq experiments. The plants were grown under their optimal growth conditions (11 h light and 13 h dark with the day temperature of 29°C and the night temperature of 24°C) in a growth chamber at the University of Massachusetts Amherst.

Histology and seed shattering measurements

Histological analyses were conducted following established procedures as described in Li et al. (2024). Briefly, florets were collected at three distinct developmental stages for each genotype: approximately 1 week before the heading date (BFT), the heading date (FT) and 1 week after the heading date (LFT). Heading date refers the stage when the first panicle has emerged halfway from the boot. One week before the heading date, most florets at the base of the panicle reach approximately 0.5 cm in length and exhibit a pale or transparent green color. These florets were fixed in an FAA solution (37% formaldehyde:ethanol:H 2 O:acetic acid = 10:50:35:5), dehydrated in a graduated ethanol series, cleared with Histoclear, embedded in paraplast, and eventually sectioned into 10 lm longitudinal slices with a Leica RM2125microtome (Leica Biosystems, Dear Park, IL, USA). The resulting sections were stained with 1% toluidine blue and observed using a Leica DM750 LED biological microscope (Leica Microsystems, Wetlzlar, Germany). High-quality images were captured using an AmScopre MU1000-HS camera.

Seed shattering levels were evaluated by measuring the breaking tensile strength (BTS), the force required to separate the seeds from the plants. This measurement was performed 30 days after heading (LFT30), a point at which the BTS values stabilize (Thurber et al., 2010(Thurber et al., , 2011)). For each panicle, 10 florets were assessed and at least two panicles per genotype were included in the measurement. In addition to measuring BTS at LFT30, we also conducted BTS measurements at FT and LFT to capture the trends in BTS changes during development. Notably, during the BFT stage, the spikelet is exceptionally fragile, preventing the measure of BTS from this stage.

RNA isolation, RNA-seq library construction and sequencing

Florets were harvested from the four rice genotypes at the three developmental stages -BFT, FT and LFTduring the afternoons between 2:00 and 4:00 p.m. For each developmental stage, we collected three replicates of three different tissues: AZ, the tissue toward floret side (F), and the tissue toward side (P) (Figure S2). In each replicate, we combined tissues dissected from around 40-60 florets, typically sourced from at least two different plants grown in the same pot. The manual tissue dissection was carried out in RNAlater stabilization solution (Invitrogen, Waltham, MA, USA) using a scalpel under a Leica dissecting microscope (Leica Microsystems, Wetzlar, Germany), following the protocol described in Yu et al. (2020). The dissected tissues were promptly transferred into RNAlater stabilization solution, stored at 4°C overnight, and subsequently preserved at À80°C before extraction. Dissected tissues were manually ground into powder using pellet pestles in liquid nitrogen. We performed RNA extraction using the PicoPure RNA Isolation Kit (Thermo Fisher Scientific, Waltham, MA, USA) and included in-column DNase I (QIAGEN, Hilden, Germany) treatment. The concentration and integrity of the RNA samples were assessed using the Nanodrop2000c (Thermo Fisher Scientific, Waltham, MA, USA) and the Fragment Analyzer System with the DNF-471 RNA kit (15 nt) (Agilent, Santa Clara, CA, USA), respectively. For library preparation, we followed the manufacturer's instructions for the QuantSeq 3 0 mRNA-Seq FWD Library Prep Kit (Lexogen, Vienna, Austria) (Moll et al., 2014). A total of 108 libraries (4 genotypes 9 3 tissues 9 3 developmental stages 9 3 replicates) were single-end sequenced on a NextSeq 500 instrument at Lexogen. The raw sequencing data was deposited in the Sequence Read Archive under the BioProjectID of PRJNA1182773.

Reads alignment and principal component analyses

The raw reads for the four rice genotypes were cleaned with Trimmomatic using a sliding window trimming method with an average base quality threshold set at a minimum of 20 (Bolger et al., 2014). The resulting clean reads from the aus cultivar and its derived weedy rice BHA were mapped to the indica Shuhui 498 (R498) reference genome (Du et al., 2017) using Hisat2 (Pertea et al., 2016). The mapped reads were subsequently used as input for HTSeq-count (Anders et al., 2015) to calculate gene counts. The clean reads from the temperate japonica cultivar and its derived weedy rice SWR were mapped to the Nipponbare japonica reference genome (Kawahara et al., 2013), following the same method as described above, to obtain gene counts for subsequent analyses.

Additionally, we also mapped the clean reads from the temperate japonica cultivar and its derived weedy rice SWR to the R498 reference genome and conducted Principal Component Analysis (PCA) using the gene counts mapped to the same reference genome for all four genotypes to gain a comprehensive understanding of the expression pattern across the four genotypes.

Weighted gene co-expression network analyses

We performed weighted gene co-expression network analyses (WGCNA) for each cultivar and its corresponding derived weed pair (BHA and the aus cultivar; SWR and the temperate japonica cultivar) using all three tissues at each developmental stage to detect modules. Following the application of a variance-stabilizing transformation to the gene counts using the vst function in DESeq2 (Love et al., 2014), genes exhibiting low variance (interquartile range <0.5) were filtered out using the VarFilter function in the geneFilter package (Gentleman et al., 2023). The resulting gene counts were used for sample clustering and the construction of gene co-expression networks using the WGCNA package (Langfelder & Horvath, 2008). We employed the pickSoft-Threshold function to choose the optimal power for the network and blockwiseModules function to construct a signed network with a mergeCutHeight of 0.25, a corType of "bicor" with max-POutliers = 0.05 and robustY = FALSE, and a minModuleSize of 30. The Pearson correlation coefficient between each module and AZ region in both genotypes or in the weedy genotype exclusively, as well as the associated statistical significance, were calculated using the cor and corPvalueStudent functions. Intramoduler connectivity and signed eigengene-based connectivity, also known as module membership, for each gene within their respective modules were determined using intramodular and signedKME functions.

Differential expression analyses

We performed differential gene expression analyses between cultivated rice and its corresponding weedy rice in each tissue at each developmental stage using DESeq2 (Love et al., 2014). Genes with an average gene count greater than 10 in any specific tissue type at any developmental stage and in either genotype, cultivated or weedy, were included in the analyses. We applied a minimum two-fold change and an adjusted P-value below 0.05 as thresholds to determine genes with significantly differentiated expression.

Gene function and GO enrichment analyses

Gene Ontology (GO) enrichment analyses were performed using the genes within each AZ-associated module and the genes exhibiting differential expression in the AZ between cultivated rice and its corresponding derived weedy rice within each lineage using agriGO v2.0 with the setting of Hypergeometric as statistical test method, Hochberg (FDR) as multi-test adjustment method, 0.05 as significance level and one as minimum number of mapping entries (Tian et al., 2017). We referred to GO terms sourced from PLAZA 5.0 (Van Bel et al., 2022) to elucidate the functions of these genes.