Genetic insights into developmental variations of spiny bracts among hazels through the pangenome construction

Abstract

Hazels (Corylus L. in birch family) are globally celebrated for their delectable nuts (Molnar, 2011). Corylus includes approximately 18 species that are widely distributed across the temperate Northern Hemisphere and is classified into four sections (Acanthochlamys, Siphonochlamys, Colurnae and Phyllochlamys) primarily based on variable bract (Figure 1a). Acanthochlamys, including two C. ferox varieties, is characterized by the spiny structure developed by bract lobes, resembling chestnut fruits, while bract lobes in remaining sections lack of this characteristic. At present, hazelnuts from different cultivars are consumed across 34 countries, with an annual yield of 1.1 million tons of in-shell nuts (FAOSTAT, 2023). The primary cultivars are domesticated from C. avellana, and hybridization between some local germplasms (e.g. C. heterophylla) and introduced C. avellana also breed some cultivars with strong adaptability to local environment. With advancements in sequencing, more crops are utilizing pangenome to capture crucial variations—especially genomic structural variations (SVs)—responsible for adaptation and agronomic trait enhancement in molecular breeding (Chen et al., 2023). However, only four hazel species currently possess high-quality reference genomes, leaving vast wild hazel species unconcerned (Brainard et al., 2023; Li et al., 2021). These unconcerned hazels also have some valuable traits. For example, C. ferox has extremely abundant nuts in each infructescence, while the spiny bract impedes its utilization for breeding (Figure 1a).

To demonstrate developmental variations of spiny bracts, we assembled high-quality chromosome-level genomes of eight wild hazels from seven species, including C. chinensis, C. fargesii, C. wulingensis, C. yunnanensis, C. kweichowensis, C. heterophylla and two C. ferox samples (Table S1). These hazels were collected from diverse geographical locations in China and assembled with the long sequencing reads (HiFi or ONT), NGS and Hi-C approaches. Except for one C. ferox sample with a contig N50 of 1.44 Mb, all of them exhibited a contig N50 longer than 9 Mb (Table S3). Eleven chromosomes were successfully anchored in each genome, revealing significant high collinearity and no large inter-chromosome rearrangement (Figure 1b). All the genomes showed a high quality by different assessments (Table S2). With six published genomes from C. avellana, C. americana, C. heterophylla and C. mandshurica (Table S3), 14 genomes cover four sections in Corylus. We predicted a range of 164–222 Mb for repeat elements and 22 137 to 28 267 for protein-coding genes across different genomes (Tables S4 and S5). Based on these genomes, we further revealed a highly supported topology, in which Acanthochlamys splits first at ~24.2 million years ago (Mya), close to splitting times between Siphonochlamys and the clade comprising Phyllochlamys and Colurnae (~21.8 Mya) and between Phyllochlamys and Colurnae (~18.7 Mya).

We used these 14 genomes to construct the Corylus pangenome (Figure 1c,d). A gene family based pangenome was firstly constructed, and a total of 10 583 core gene families were identified. and they mainly involved in maintenance the basic process, such as organic substance metabolic process and biological regulation (Figure S1). A total of 13 702 dispensable (shared by at least two samples) and 8048 private gene families were also identified (Table S6). The dispensable genes were enriched in stress resistance functions (Figure S2). Private gene number ranges from 73 to 1559 for different species (Table S7). The graph-based pangenome constructed by PanPop (Zheng et al., 2024) comprised 3 334 456 nodes, 6 050 688 edges, 848 361 non-reference nodes, and approximately 440 Mb of non-reference sequences. Full SVs were further identified within the constructed pangenome. Approximately 601 k non-redundant SVs were detected, totalling around 320 Mb (Figure 1e,f). Among these non-redundant SVs, 32.4% were shorter than 100 bp; but there were still 3650 SVs exceeding 10 kb in length (Figure S3a). Interestingly, SVs occurring in non-repeat regions were significantly more abundant than those in repeat regions (377 k vs. 166 k, P < 10⁻¹⁵, Figure S3b). Similarly, we found 90% (1035) of (resistance) R-genes with at least one SV located within gene region or in the up/downstream regions within a distance of less than 2 kb, and found 41.3% of these genes exhibited significant different expression (Figure S4). These results suggest that SVs should have high impacts on gene functions in hazels.

We then primarily demonstrated the SVs associated with bract development in hazels. Acanthochlamys showed a widely distribution range and bract development from simple lobes to spiny structure (Figure 1a,h). We firstly identified 3582 genes influenced by the C. ferox-specific SVs, while genes related to these SVs showed no significant enrichment of Gene Ontology (GO) terms. This could be attributed to variations responsible for spiny bracts and the other traits. So, we performed developmental transcriptome analyses of C. ferox, C. heterophylla and C. fargesii, to identify crucial genes associated with the development of spiny bracts. The transcriptome samples were collected from a sympatric location—resembling a natural common garden experiment, and another allopatric location. For early developmental stages (less than 90 days after flowering, DAF), we sampled every 20 days, while for the mature stages (>90 DAF), we sampled every month (Table S8). Using weighted gene co-expression network analysis (WGCNA) analyses, we identified 23 modules with sizes ranging from 44 to 2181 genes (Figure S5). Among them, magenta, yellow and blue modules were significantly overexpressed in C. ferox at 57, 92 and 158 DAF, respectively (Figure S6). These three modules correspond to three distinct and important developmental stages of bracts that meticulously documented through photography and were categorized into three distinct phases according to the bract development in C. ferox: branching (around 57 DAF), elongation (around 92 DAF) and hardening (128 ~ 158 DAF) (Figure 1h). The functional enrichment of genes in these models also corresponding the different development stage. In the magenta module (branching stage), functional genes are mainly associated with plant epidermis morphogenesis and biosynthesis of cuticle, cellulose and fatty acid (Table S9). Genes in yellow module (elongation stage) mainly involved in cell wall macromolecule biosynthesis, such like lignin, cellulose, pectin and D-xylose (Table S10), while those in the hardening stage (blue module) only related to obsolete vacuolar (Table S11).

As the initial branching stage determined the bract phenotype (Figure 1h), we further focused on this stage. We retrieved differentially expressed genes (DEGs) between C. ferox and other hazels, and 140 DEGs showed the significantly overexpressed in this stage for C. ferox, with 18 DEGs exhibiting a fold change >10 times. Three of these 18 DEGs contained the unique SVs of C. ferox in their 2-kb flanking regions and were also members of the magenta module (Table S12). CHE07516 and CHE10948 belong to the peroxidase superfamily and sulphite exporter TauE/SafE family, respectively. However, there is no detailed report about phenotypes of these two genes. Nevertheless, CHE05907 (designated as CheCUC3) was annotated as a homologue of CUC3, belonging to the NAC family. In Arabidopsis, CUC3 is known to regulate leaf morphology, and in cucumbers, it is associated with the spiny trait (Chen et al., 2014). In C. ferox, CheCUC3 showed a higher expression in the in the basal (represent the branching tissue) bract than in top regions by qRT-PCR with 57 DAF bracts, indicated that this gene likely plays an important role in the development of spiny bracts (Figure 1i). Two unique SVs—1.2-kb and 180-bp deletions—were together located at 890 bp upstream and 1.2-kb downstream of CheCUC3 in C. forex (Figure 1g). Then, we evaluated the impact of the SV in upstream promote regions on the gene function by isolating promoter fragments from the three species and performing a dual-luciferase reporter assay (LUC) in transiently transformed tobacco plants. LUC activity driven by the promoters of C. ferox was significantly higher than that from the promoters of C. heterophylla and C. fargesii (Figure 1j–l). Therefore, these results demonstrated that the unique deletion in C. ferox was accelerated the expression of CheCUC3 and it should facilitate the development of spiny bracts.

In summary, our comprehensive Corylus pangenome establishes a valuable resource for researchers and breeders to investigate and harness genes present in both cultivated and wild hazels. Through development transcriptome analysis, we revealed the function of each stage of spiny growth. Furthermore, we have identified and verified a significant structure variant that influences the development of the distinctive bracts in hazels. The genomic resources of R-gene and relative genes of spiny development could aid further breeding. Our findings contribute to a deeper understanding of hazelnut biology and provide a foundation for further research and breeding programmes aimed at improving hazelnut cultivation and production.