Tight genetic linkage of genes causing hybrid necrosis and pollinator isolation between young species

The mechanisms of reproductive isolation that cause phenotypic diversification and eventually speciation are a major topic of evolutionary research. Hybrid necrosis is a post-zygotic isolation mechanism in which cell death develops in the absence of pathogens. It is often due to the incompatibility between proteins from two parents. Here we describe a unique case of hybrid necrosis due to an incompatibility between loci on chromosomes 2 and 7 between two pollinator-isolated Petunia species. Typical immune responses as well as endoplasmic reticulum stress responses are induced in the necrotic line. The locus on chromosome 2 encodes ChiA1, a bifunctional GH18 chitinase/lysozyme. The enzymatic activity of ChiA1 is dispensable for the development of necrosis. We propose that the extremely high expression of ChiA1 involves a positive feedback loop between the loci on chromosomes 2 and 7. ChiA1 is tightly linked to major genes involved in the adaptation to different pollinators, a form of pre-zygotic isolation. This linkage of pre- and post-zygotic barriers strengthens reproductive isolation and probably contributes to rapid diversification and speciation.

The mechanisms of reproductive isolation that cause phenotypic diversification and eventually speciation are a major topic of evolutionary research. Hybrid necrosis is a post-zygotic isolation mechanism in which cell death develops in the absence of pathogens. It is often due to the incompatibility between proteins from two parents. Here we describe a unique case of hybrid necrosis due to an incompatibility between loci on chromosomes 2 and 7 between two pollinator-isolated Petunia species. Typical immune responses as well as endoplasmic reticulum stress responses are induced in the necrotic line. The locus on chromosome 2 encodes ChiA1, a bifunctional GH18 chitinase/lysozyme. The enzymatic activity of ChiA1 is dispensable for the development of necrosis. We propose that the extremely high expression of ChiA1 involves a positive feedback loop between the loci on chromosomes 2 and 7. ChiA1 is tightly linked to major genes involved in the adaptation to different pollinators, a form of pre-zygotic isolation. This linkage of pre-and post-zygotic barriers strengthens reproductive isolation and probably contributes to rapid diversification and speciation.
Speciation depends on the evolution of reproductive barriers that reduce gene flow between previously interbreeding populations [1][2][3] . In plants, pollinator-mediated isolation constitutes a crucial pre-zygotic isolation barrier 4 . The preferences of distinct classes of pollinators for sets of floral traits such as colour, morphology, scent and nectar production are considered a driving force in the rapid diversification of the angiosperms 5 . These matching sets of floral traits are called pollination syndromes. However, pollinator preference for a flower type is rarely absolute, and hybridization between (incipient) species is often observed. Thus, plant speciation generally involves the accumulation of multiple reproductive barriers 3 .
In plants, hybrid necrosis (HN) is a post-zygotic isolation barrier in which hybrids show necrotic leaves and poor growth 6,7 . Many of the reported cases were caused by deleterious epistatic interactions between alleles at two or more loci, creating a post-mating barrier that reduces fitness in the hybrids but not in the parents 8 . HN often involves genes encoding components of the immune system 9,10 . For instance, many reported cases of HN in Arabidopsis are due to improper combinations of nucleotide-binding domain leucine-rich repeat (NLR) proteins, which are key players in the immune systems of both plants and vertebrates 11-14 . Article https://doi.org/10.1038/s41477-023-01354-8

An introgression line shows necrotic symptoms
The introgression line IL5 is derived from a cross between hawkmothpollinated P. axillaris N and hummingbird-pollinated P. exserta 28 . It segregates for a P. axillaris region of low recombi nation on chromosome 2 in an otherwise P. exserta background (Fig. 1a). Several major genes affecting pollination syndrome traits were previously identified in this so-called supergene region [19][20][21][22] . Unexpectedly, the progenies of IL5 that are homozygous for P. axillaris N at the supergene region (hereafter called IL5-Ax) displayed necrotic symptoms that strongly affect plant growth (Fig. 1b) and flower production (Fig. 1c). IL5-Ex (homozygous for P. exserta at the supergene region) plants were healthy. The necrotic symptoms in IL5-Ax can be observed as early as 38 days after sowing (Extended Data Fig. 1a). The symptoms were also observed in plants heterozygous for the introgression (IL5-Het) but were milder than in the IL5-Ax plants, suggesting semi-dominance (Fig. 1d). Necrosis was associated with the production of reactive oxygen species (ROS) and increased levels of expression of key pathogenesis-related genes ( Fig. 1e and Extended Data Fig. 1b,c). These data are consistent with HN, in which an immune response is triggered in the absence of pathogens 10 . Endoplasmic reticulum (ER) stress marker genes 29 , for example, two luminal binding protein coding genes BiP4 and BiP5, and transcription factor bZIP60, were highly induced in the necrotic line (Extended Data Fig. 1d), suggesting ER-stress-induced cell death may play a role in HN development 30 .
We have previously studied the genetic basis of the evolution of pollination syndromes in the South American genus Petunia (Solanaceae). Petunia axillaris (P. axillaris) is widely distributed in southern South America, while Petunia exserta (P. exserta) is a highly endemic species that is found exclusively in shaded shelters in the Guaritas Region of Rio Grande do Sul, Brazil 15 . P. axillaris has white, ultraviolet (UV)-absorbent, fragrant flowers and is pollinated by hawkmoths 16,17 , while P. exserta has red, scentless flowers with exserted stamens and stigma and is pollinated by hummingbirds 15,18 . Single genes underlying major quantitative trait loci (QTLs) for UV absorption (MYB-FL), scent production (CNL1) and pistil length (EOBII) were found to be tightly linked in a so-called supergene region on chromosome 2 [19][20][21][22] . Such tight genetic linkage is thought to impede the dissolution of pollination syndromes by recombination. The two species are thought to have evolved recently from a common ancestor that was most likely hawkmoth pollinated [23][24][25] . Natural hybrids are present in a few locations in the Guaritas Region where the two species are sympatric, showing that reproductive barriers are incomplete 18,25,26 .
In this Article, we describe a case of HN in a cross between P. axillaris and P. exserta. The necrotic symptoms are associated with a deleterious combination between two loci located on chromosome 2 and chromosome 7. We show that the locus on chromosome 2 encodes ChiA1, a bifunctional chitinase/lysozyme with a major role in pattern-triggered immunity against fungi and bacteria 27 . We discuss the potential relevance of the tight linkage of ChiA1 to pre-zygotic isolation barriers.  Article https://doi.org/10.1038/s41477-023-01354-8

HN is caused by genetic interaction of two loci
HN often results from incompatibilities between genetic loci that have diverged over time and are brought together in hybrids 10 . Such interactions can involve pairs of nuclear genes or nucleo-cytoplasmic interactions 31,32 . Therefore, we generated reciprocal F2 populations derived from P. axillaris N and P. exserta to dissect the genetic architecture of HN and identify the loci interacting with the chromosome 2 region. Both populations displayed comparable distributions of HN scores, that is, approximately three-quarters of the plants did not show symptoms and one-quarter of the plants displayed symptoms. The similar distribution of HN scores in the two populations indicates that the loci controlling HN are nuclear, as cytoplasmic inheritance is unidirectional (Fig. 2a,b). Bulked segregant RNA sequencing (BSR-seq) was used to locate the genomic regions causing HN. Strong signals associated with the necrotic phenotype were detected on chromosomes 2 and 7 ( Fig. 2c and Extended Data Fig. 2). We named these loci HNe2 and HNe7. The necrosis is triggered by a combination of P. axillaris N alleles at HNe2 and P. exserta alleles at HNe7 (Fig. 2c and Extended Data Fig. 2) in accordance with the genetic make-up of IL5-Ax (Fig. 1a).

Map-based cloning of HNe2
To reduce the interval carrying HNe2, a fine-mapping approach was undertaken through the search of recombinant progenies of IL5-Het. The relationship between the phenotypic segregation of HN and the genotype of the recombinant lines led to a first interval of 8.7 Mb (Fig. 3a and Supplementary Table 1). Then, the recombination breakpoints of ten of the most informative lines were precisely characterized through whole-genome sequencing, resulting in a 1.74 Mb interval ( Fig. 3a and Extended Data Fig. 3). As this region still contains 60 genes (Supplementary Table 2), an RNA sequencing (RNA-seq) analysis of leaf tissues from IL5-Ax and IL5-Ex plants was performed to investigate their functional relevance on the basis of expression. Among the 34 expressed genes in the candidate region (read counts >100 in at least one of the samples; Supplementary Table 3), 8 were differentially expressed (q < 0.001) including 4 genes known to play a role in plant-pathogen interactions: ChiA1, ChiA2, FMO1 and CNGC1 (Table 1  and Supplementary Table 4). Chitinase A1 (ChiA1) is one of the most highly expressed genes in IL5-Ax (Table 1 and Extended Data Fig. 4a) and is annotated as a bifunctional endochitinase 33 (Extended Data Fig. 4b). The protein encoded by ChiA1 is homologous to AtLYS1/ChiA (AT5G24090), which has a central role in triggering immune responses in Arabidopsis 27 (Extended Data Fig. 4c). Importantly, ChiA1 carries a nonsense mutation in P. exserta, which leads to a truncated protein of 158 amino acids (Extended Data Fig. 4d,e). ChiA2 is closely related to ChiA1 but shows a much lower expression level than ChiA1 in IL5-Ax (Table 1). Flavin-containing monooxygenase 1 (FMO1) plays key roles in systemic acquired resistance by synthesizing N-OH-Pip from pipecolic acid 34 . It has two copies in the P. exserta genome, whereas only one copy was present in P. axillaris N. Cyclic nucleotide-gated ion channel 1 (CNGC1) belongs to a gene family involved in plant immunity by triggering calcium signalling and hypersensitive response 35,36 .
To functionally validate these four candidates for HNe2, we knocked down their transcript levels in IL5-Ax by virus-induced gene silencing (VIGS) 37 . VIGS of ChiA1 caused a robust decrease of necrosis in IL5-Ax plants compared with the empty vector or untreated controls (Fig. 3b,c). In contrast, VIGS of ChiA2, FMO1 and CNGC1 showed    Tables 1-4. b,c, VIGS of the four selected candidate genes for HNe2: representative leaf phenotypes from one branch arising after the ChiA1 VIGS treatment and empty control VIGS treatment in IL5-Ax; leaves are shown from a top-to-bottom order; scale bar, 1 cm (b). Yellow area ratio of IL5-Ax plant leaves after VIGS of ChiA1, ChiA2, FMO1, CNCG1 and empty control; five biological replicates were used for each treatment (n = 5 plants). Symptoms were quantified on the basis of the ratio of the area with yellow colour compared with the whole leaf area. A higher value indicates a higher proportion of yellow surface. Different letters indicate significant differences (P < 0.05, one-way ANOVA, Tukey's honestly significant difference (HSD) test); for P values, see source data (c). d,e, Transient overexpression of ChiA1 Ax in IL5-Ex plant leaves: A representative phenotype of IL5-Ex leaf after agro-infiltration with 35S::ChiA1 Ax or 35S::ChiA1 Ex as a control; scale bar, 1 cm (d). Yellow area ratio of IL5-Ex plant leaves after agro-infiltration of 35S::ChiA1 Ax , 35S::ChiA1 Ex , 35S::GFP and empty infiltration buffer. Three individual plants were used for each treatment. ChiA1 Ex , GFP and infiltration buffer were used as negative controls. The phenotype was analysed 2 weeks after the infiltration. Different letters indicate significant differences (P < 0.05, one-way ANOVA, Tukey's HSD test); for P values, see source data (e). f, Diagram showing the predicted protein products of ChiA1 gene in IL5-Ax, IL5-Ex and chia1-cas9-1. Red box indicates the mismatched protein sequence caused by frameshift. The numbers indicate the amino acid sites. SP, signal peptide; GH18, glycosyl hydrolase 18 domain. g,h, Plant phenotype (g) and leaf yellow area ratio (h) of chia1-cas9-1 in IL5-Ax (middle) compared with IL5-Ax (left) and IL5-Ex (right). For the leaf yellow area ratio, three individual plants were analysed for each genotype. Different letters indicate significant differences (P < 0.05, one-way ANOVA, Tukey's HSD test); for P values, see source data. Article https://doi.org/10.1038/s41477-023-01354-8 no significant reductions of symptoms (Fig. 3c). To confirm the role of ChiA1 in HN, we transiently overexpressed ChiA1 Ax (P. axillaris N allele) as well as ChiA1 Ex (P. exserta allele) in IL5-Ex leaves driven by the 35S promoter. ChiA1 Ax was highly expressed, while ChiA1 Ex was barely detectable (Extended Data Fig. 4f). Only leaves overexpressing ChiA1 Ax turned yellow in the region of infiltration (Fig. 3d,e). Thus, the lack of messenger RNA accumulation is post-transcriptional, most likely due to nonsense-mediated decay (NMD) 38 . Furthermore, we generated a ChiA1 loss-of-function mutant in the IL5-Ax background by clustered regularly interspaced short palindromic repeats (CRISPR)-Cas9. The mutant (chia1-cas9-1) has a 4 bp deletion, leading to a frameshift and a premature stop codon, resulting in a truncated protein of 173 amino acids (instead of the 292 amino acid full length), similar to the protein encoded by ChiA1 Ex (Fig. 3f). The necrotic symptoms were strongly reduced, and plant morphology was similar to IL5-Ex (Fig. 3g,h and Extended Data Fig. 5a,b). The combined results prove beyond reasonable doubt that ChiA1 is the gene underlying HNe2.

HN does not require ChiA1 enzymatic activity
ChiA1 harbours the glycoside hydrolase 18 (GH18) domain and a secretion signal peptide (Extended Data Fig. 4b) and is probably a bifunctional enzyme digesting both chitin and peptidoglycan as the homologous proteins in rubber tree (Hevea brasiliensis) and Arabidopsis 27,39 . We observed that chitinase and lysozyme activities were extremely high in IL5-Ax leaves compared with the other lines ( Fig. 4a). In the chia1-cas9-1 mutant, these activities were reduced to a level similar to that in IL5-Ex. P. axillaris N plants overexpressing ChiA1 (35S::ChiA1) showed higher activity levels compared with wild-type P. axillaris N. Moreover, the enzyme activities were in line with the transcription levels of ChiA1 (Fig. 4b). These results indicate that ChiA1 is responsible for the high chitinase and lysozyme activities in IL5-Ax leaves.
As ChiA1 degrades microbial cell wall components, we asked whether the necrotic symptoms might be caused by conserved degradation products of microbial cell walls that could act as microbe-associated molecular patterns 40 . When grown under axenic conditions, IL5-Ax plants still developed strong necrotic symptoms (Fig. 4c,d; and Extended Data Fig. 6). This suggests that necrotic symptoms develop in the absence of a microbial substrate for ChiA1. To test whether symptom development relied on ChiA1 enzymatic activity, we generated a ChiA1 mutant in which the critical amino acids D148 and E150 in the glycoside hydrolase active site DXDXE motif are replaced by alanine (Fig. 4e). Mutation of these two amino acids abolishes both chitinase and lysozyme activities [41][42][43] (Extended Data Fig. 7a). IL5-Ex leaves overexpressing either the catalytically inactive ChiA1 D148A, E150A or ChiA1 Ax showed strong necrosis at the site of infiltration, whereas the overexpression of ChiA1 Ex and GFP did not induce symptoms ( Fig. 4e and Extended Data Fig. 7b,c). These results show that the necrotic symptoms represent a case of autoimmunity that does not rely on the enzymatic activities of ChiA1.

ChiA1-dependent induction of WRKY18 induces necrotic symptoms
Transcriptional reprogramming is an integral part of plant immunity 44 . To better understand the mechanism underlying the activation of the immune system in IL5-Ax, we screened the promoters of the 1,717 differentially expressed genes (DEGs; q < 0.001) between IL5-Ax and IL5-Ex for conserved cis-elements (Supplementary Table 5). A WRKY binding site (Fig. 5a) was present in 60% of the DEG promoters (Supplementary Table 6). WRKY transcription factors are involved in pathogen responses, abiotic stress responses and senescence 45,46 . Of all the DEGs encoding WRKY transcription factors, WRKY18 (homologous to AtWRKY18) showed the highest expression level in IL5-Ax, while its expression was negligible in IL5-Ex or P. axillaris N wild-type plant leaves (Fig. 5b). Moreover, AtWRKY18 binds to the typical binding motif in its target genes 47 (Fig. 5a). Overexpression of WRKY18 induced necrotic symptoms in leaves of IL5-Ex as well as Nicotiana benthamiana (Fig. 5c,d). Importantly, agro-infiltration of 35S::ChiA1 Ax but not 35S::ChiA1 Ex in IL5-Ex leaves induced WRKY18 expression (Figs. 5e and 3d and Extended Data Fig. 4f). We conclude that the massive transcriptional reprogramming observed in IL5-Ax is mediated through ChiA1-dependent induction, most likely indirectly, of one or more WRKY18-type transcription factors.

The synteny of ChiA1 and MYB-FL is specific to Petunia
ChiA1 is located at the edge of a region of chromosome 2, where recombination in crosses between P. axillaris N and P. exserta is strongly reduced 19 (Fig. 6a). Indeed, ChiA1 is tightly linked (1.63 cM) to MYB-FL, a well-characterized gene inside the supergene that is responsible for the gain and loss of floral UV absorption during evolutionary shifts in pollinator preference 21 (Fig. 6a and Supple mentary Table 7).
Analysis of ChiA1 and its neighbouring genes showed clear microsynteny between Petunia and the related Solanaceae species Solanum lycopersicum (tomato) and S. tuberosum (potato). However, in both these Solanum species, MYB-FL is not linked to ChiA1 but resides on a different chromosome (Extended Data Fig. 8 and Supplementary Table 8). The exclusive genetic linkage between ChiA1 and MYB-FL in Petunia is consistent with the notion that the ChiA1-mediated HN constitutes a post-zygotic barrier that acts in concert with the pre-zygotic barrier to limit gene flow.

The ChiA1 Ex allele is widespread in nature
To analyse the frequencies of the two alleles of ChiA1, we genotyped ChiA1 across the range of the two species (for details, see Supplementary Tables 9 and 10). P. exserta was sampled from 18 sites representing its highly restricted distribution in the Guaritas Region of Brazil. Of the 75 P. exserta accessions genotyped, 58 were homozygous for the nonsense mutation (ChiA1 Ex ), 13 were homozygous for the P. axillaris genotype (ChiA1 Ax ) and 4 were heterozygous (Fig. 6b, right, Extended Data Fig. 9 and Supplementary Table 9). The overall ChiA1 Ex allele frequency among these accessions was 80% (Fig. 6b, right). P. axillaris accessions were sampled from multiple locations in Argentina, Uruguay and Brazil, including the regions of sympatry with P. exserta (Fig. 6b, left, red-dashed square). Of the 238 P. axillaris accessions, 173 were homozygous for the P. axillaris allele (ChiA1 Ax ), 61 were homozygous for the P. exserta allele (ChiA1 Ex ) and 4 were heterozygous (Fig. 6b, left, and Supplementary Table 10). The overall ChiA1 Ax allele frequency among these accessions was 74% (Fig. 6b, left). Sixty-two P. axillaris accessions with the ChiA1 Ex allele are located in the Corrientes and Entre Ríos regions in Argentina (Fig. 6b, solid square in red), which is ~500 km away from the Guaritas Region where P. exserta accessions are exclusively found. These accessions belong to P. axillaris subspecies parodii 48 , which is distinguished from P. axillaris ssp. axillaris mainly by its longer corolla tube 49 .
These results confirm the presence of both ChiA1 alleles in the natural P. axillaris and P. exserta populations. ChiA1 Ex is more present in P. exserta, whereas ChiA1 Ax is more frequent in P. axillaris.

Discussion
We discovered a case of HN in crosses between P. axillaris and P. exserta relying on a deleterious epistatic interaction between the HNe2 and HNe7 loci (Fig. 6c). The gene underlying HNe2 is ChiA1, which encodes a protein with both chitinase and lysozyme activity. ChiA1 Ax allele encodes a functional enzyme, whereas ChiA1 Ex has a nonsense mutation leading to a premature stop codon.
Cases of HN characterized in multiple species often have commonalities: activation of immune responses and the direct participation of immune system components 10 . In a well-studied case in tomato, a secreted cysteine protease causes autoimmunity when the allele from a related wild species is present in combination with an avirulence (Avr) receptor of the cultivated tomato 50,51 . In this example, however, the protease activity of the enzyme is thought to be essential, whereas replacement of the critical amino acids in the active site of ChiA1 does not interfere with its ability to cause necrosis (Fig. 4e).
If it is not the activity of the protein, by what mechanism does ChiA1 cause necrosis? The level of expression is high for ChiA1 Ax in IL5-Ax but low for the ChiA1 Ex and chia1-cas9-1 alleles. The low expression of ChiA1 Ex and chia1-cas9-1 is mostly probably caused by the mRNA surveillance pathway through the mRNA NMD mechanisms 38 . In contrast, ChiA1 Ax is one of the most highly expressed genes in the necrotic IL5-Ax background. The high expression of ER stress marker genes in the same background indicates ER stress. Many key pathogenesis-related proteins rely on ER and ER quality control for proper folding and secretion [52][53][54] . We propose that the high and chronic demand for the processing of ChiA1 in the ER exceeds the ER quality control working capacity and causes prolonged ER stress and cell death [55][56][57]  Article https://doi.org/10.1038/s41477-023-01354-8 seems plausible (Fig. 6d). The basal expression of ChiA1 Ax in the necrotic line may trigger the activation of HNe7 Ex (probably via WRKY18), and thus HNe7 Ex induces autoimmunity, which further enhances the expression of ChiA1. In this case, HNe7 may encode an auto-active immune system component (for example, NLR 10 ). Alternatively, HNe7 Ex activates the expression of ChiA1 Ax , and ChiA1 Ax induces the autoimmunity. In that case, HNe7 might be a transcriptional or post-transcriptional regulator. HNe7 is located in a large region of very low recombination containing close to 1,000 genes, and its identity is not yet known. Thus, the molecular details of the interaction of ChiA1 and HNe7 remain to be determined. A unique aspect of the work presented here is that ChiA1 is tightly linked to a region of low recombination containing major genes involved in the evolution of pollination syndromes. Among them is MYB-FL, encoding a transcription factor that induces the synthesis of UV-absorbing pigments in the moth-pollinated P. axillaris but was inactivated in the hummingbird pollinated P. exserta 19,21 . This genetic linkage is relatively recent, as ChiA1 and MYB-FL are located on different chromosomes in tomato and potato (Extended Data Fig. 8). This association may have served to enhance reproductive isolation during the process of speciation.
Both the ChiA1 Ax and ChiA1 Ex alleles were found in P. axillaris as well as in P. exserta wild accessions collected from multiple sites in South America (Fig. 6b). The presence of ChiA1 Ax in wild P. exserta could be due to introgression 58 . P. axillaris ssp. parodii does not reach the P. exserta region today, and in the past, there were no conditions suitable for P. parodii in Serra do Sudeste, even considering Last Glacial Maximum (~ 22 kya) or mid-Holocene (6 kya) 59 . The presence of ChiA1 Ex in P. axillaris ssp. parodii strongly suggests incomplete lineage sorting 60 , consistent with a function in reproductive isolation during speciation. As chitinases are important in pattern-triggered immunity against pathogens, the loss of ChiA1 function is likely to compromise defence. The potential benefit of the ChiA1 Ex allele in allopatric P. axillaris ssp. parodii populations will need further study. The presence of ChiA1 Ax in P. exserta suggests conditions in which the benefits of pathogen resistance outweigh its effect on reproductive isolation.
Speciation occurs when reproductive barriers accumulate and substantially reduce gene flow between lineages 3 . Pollinator preference is a strong barrier against gene flow, but it is rarely absolute 61 . This is also true for Petunia 62 . Multiple reproductive barriers acting in concert are therefore often needed to complete reproductive isolation 63,64 . The genetic linkage of distinct isolation mechanisms would further enhance reproductive isolation and thereby enhance the tempo of diversification and speciation. Such linkage may be a more general phenomenon that can help explain the rapid and successful diversification of the angiosperms.

Plant materials
P. axillaris N is from the Rostock Botanical Garden (Germany), P. exserta is from R.J. Griesbach (Beltsville, United States). The accessions are maintained by self-fertilization. The introgression line IL5 has been described previously 28 . Wild accessions of P. axillaris and P. exserta have been described previously 21,23,25,48 and are described further in Supplementary Tables 9 and 10.
Plants were grown in a greenhouse or a growth chamber. Greenhouse plants were grown with additional lighting resulting in 14 h day at 18-25 °C in pots. Plants grown in a growth chamber are under a light:dark regime of 15 h:9 h, at 22 °C:17 °C at 60%-80% relative humidity, in commercial soil (70% Klasman substrate, 15% Seramis clay granules and 15% quartz sand) and fertilized once a week (Plantaktiv, 16+6+26 type K fertilizer, 0.1% concentration).

DAB staining
Diaminobenzidine (DAB) staining was performed according to a previous protocol 65 . After the treatment, the leaf was used for immediate observation or kept at 4 °C.   Tables 9 and 10. c, Schematic showing the predicted interaction of pre-zygotic and post-zygotic barriers between P. axillaris and P. exserta. ChiA1 is tightly linked with several genes controlling pollination syndromes in the supergene. The different pollination syndromes helped the species to adapt to their correspondent pollinators. The HN helps to fix the right combinations of floral traits for its pollinator and thus enhances isolation. d, Model of the molecular mechanism of HN. Left: interaction of full-length P. axillaris ChiA1 and active P. exserta HNe7 induces necrosis. This interaction proceeds through a positive feedback loop involving ER stress and WRKY transcription factor(s). Middle: truncated P. exserta or chia1-cas9-1 ChiA1 fails to induce active P. exserta HNe7 because ChiA1 mRNA does not accumulate due to NMD. Right: full-length P. axillaris ChiA1 does not cause symptoms in combination with inactive P. exserta HNe7.  Table 11.

Leaf yellowness analysis
The ratio of leaf yellow area was calculated on the basis of a previous method with modifications 67 . Briefly, leaves from each plant were removed and photographed on a black background under uniform artificial light conditions with a colour reference card included in each photo. The software Fiji (ImageJ) was used to process the pictures 68 .
Total leaf area as well as yellow and green leaf area were measured on the basis of colour thresholds for hue, saturation and brightness values. Only pictures taken in the same session were compared to guarantee light uniformity. To define the hue parameters that identify a region as yellow, we first applied the parameters to the colour reference card in each photo and verified that the yellow areas identified were consistent. The yellow area measured on the colour cards is shown in Extended Data Fig. 10. We then used the hue thresholds to analyse the leaf photos. The Fiji macro outputs a painted version of the photos, where the yellow area measured is shown. We manually verified that this area correctly overlapped the yellow parts of the leaves. Detailed scripts and instructions can be found here: https://github.com/Kuhlemeier-lab/fiji_macros.

BSR-seq and data analysis
We used the method for BSR-seq (a bulk-segregant analysis based on RNA-seq data) published by Soyk et al. 69 to map the loci responsible for HN. P. axillaris N (mother plant) and P. exserta (pollen donor) were used to produce an F2 population. A total of 384 F2 plants were sown for sampling. On the basis of the phenotypic scale shown in Fig. 2a, we selected plants 13 weeks after sowing showing a strong necrotic leaf phenotype (19 individuals, score 4) or showing a healthy leaf phenotype (89 individuals, score 0). A leaf disc (Ø = 9 mm) from a mature and healthy leaf of each plant was collected, flash frozen in liquid nitrogen and stored at −80 °C. Total RNA was extracted from pools of leaf discs either from the necrotic plants or the healthy plants. RNA-seq was performed by NovoGene with a polyA-enriched library prep, to obtain paired-end reads of 150 bp. Alignment of quality-controlled reads was performed with STAR 70 (2.6.0c) on the reference genome of P. axillaris N version 4.03. Variants were called using GATK 71 (4.0.4.0). The variants were filtered to keep only high-quality, biallelic single-nucleotide polymorphisms (SNPs) with a minimum read depth of 100 in each sample (using minDP argument in vcftools), and with minimum 100 bp between each variant position. We calculated the difference in alternate allele frequency between the two pools (necrotic -healthy) and used this ∆SNP frequency to define genome-wide thresholds at 0.01 and 0.05 (lower thresholds) and 0.95 and 0.99 (higher thresholds) quantiles. The SNPs were then grouped in stepping windows of 100 SNPs along the genome and the proportion of SNPs with value outside of the genome-wide thresholds in each window was used to produce the main figure. Allele frequency calculations and plotting was performed in R 72 . Detailed parameters, software versions and scripts are deposited on Github: https://github.com/Kuhlemeier-lab/Petunia_hybrid_necrosis#bsr-seq.

Fine mapping using recombinant lines of IL5-Het
As IL5-Het is a near-isogenic line segregating for a chromosome 2 region comprising HNe2 while the rest of the genome is homo zygous, we use it as a starting material to search for recombinant lines among its progenies and reduce the size of the interval of interest. We screened 3,800 progenies of IL5-Het with the markers mHIND-cn9140 and mHIND-MYB58 that are flanking the heterozygous region, allowing the identification of 37 recombinant lines (Supplementary Table 1). These lines harbour a smaller heterozygous region compared with IL5-Het that were genotypically characterized with the addition of more genetic markers (Supplementary Tables 1 and 11). Then, phenotypic characterization of homozygous progenies of these recombinant lines allowed the fine mapping of HNe2 down to an 8.7 Mb region (Supplementary Table 1).
To get a better resolution of the recombination breakpoints and reduce the HNe2 locus, ten informative recombinant lines and four controls homozygous for the HNe2 region (two P. axillaris and two P. exserta) were selected for whole genome sequencing (Supplementary Table 1). Genomic DNA extraction from leaf tissue of the selected lines was performed with a modified CTAB method as previously reported 73 . DNA was sequenced by the Next Generation Sequencing platform of the University of Bern using Illumina whole genome sequencing library preparation to obtain paired-end 150 bp reads. Quality-controlled reads were aligned with BWA MEM default parameters 74 to the genome of P. axillaris N v. 4.03, providing an average coverage of 3.3 × along the genome. Variants were called and quality filtered with GATK 71 . Observation of the genotypes of the variants associated with the HN phenotype in the recombinant lines allowed reduction of the size of the region containing HNe2. The detailed software versions, parameters and scripts are deposited on Github: https://github.com/Kuhlemeier-lab/ Petunia_hybrid_necrosis#il-shallow-sequencing.

RNA-seq analysis
IL5-Ax and IL5-Ex plants coming from the same IL5-Het parent were grown for leaf sampling. The leaf displaying the onset of necrosis in IL5-Ax was harvested, and equivalent leaves were sampled for the other genotypes. Three biological replicates from three different plants were set for each genotype. Total RNA was extracted from leaf tissues using an innuPREP DNA/RNA Mini Kit (Analytik Jena; code: 845-KS-20800250). Quality controls were performed using a Nanodrop (NanoDrop 1000, Thermo Fisher) and a fragment analyser (2100 Bioanalyzer Instrument, Agilent). RNA-seq was performed by Lausanne Genomic Technologies Facility (University of Lausanne) with TruSeq stranded RNA library prep to obtain single-end reads of 125 bp. Quality-controlled reads were aligned to the reference genome of P. axillaris version 4.03 using STAR 70 . Read counts were generated with Subread 75 . Differential expression analysis was performed with DESeq2 in R 76 between IL5-Ax and IL5-Ex (q < 0.001). The detailed software versions, parameters and scripts are deposited on Github: https://github.com/Kuhlemeier-lab/Petunia_hybrid_necrosis#rnaseq.

Phylogenetic analysis
ChiA1 and its homologous sequences from species listed in Extended Data Fig. 4c were identified in the NCBI nr (nonredundant protein sequences) database by performing a BLASTp search with ChiA1 protein sequences (Supplementary Table 12). Amino acid sequences were aligned with MUSCLE in the MEGA-X software package using the default settings for protein multiple alignments. Evolutionary distances were computed using Poisson correction analysis. The bootstrap method with 1,000 replicates for phylogeny testing was used.

Protein structure prediction
The structures of ChiA1 variants were predicted by AlphaFold 77 . The structures were visualized and aligned by PyMOL (Version 2.5.3, Schrödinger LLC).
Article https://doi.org/10.1038/s41477-023-01354-8 VIGS VIGS was performed as described previously 23 . For each candidate gene, we targeted a specific part of the coding region to avoid off-targeting (Supplementary Table 11). After amplification with primers containing BamHI and EcoRI restriction sites, these fragments were cloned into the pTRV2-MCS plasmid (ABRC accessions CD3-1040) and transformed into Agrobacterium tumefaciens strain GV3101. Seven weeks after the infection, only leaves from the branches arising from infected meristems were phenotyped.

Transient and stable overexpression
ChiA1 D148A, E150A coding sequence was generated by PCR mutagenesis method with ChiA1 Ax as a template. The coding sequences of ChiA1 Ax , ChiA1 D148A, E150A and ChiA1 Ex were amplified with primers containing AttB sites (Supplementary Table 11) and cloned into the Gateway-compatible binary vector pGWB402 (Addgene plasmid #74796) containing a CaMV35S promoter. For transient overexpression, the constructs were transformed into Agrobacterium tumefaciens strain GV3101. The Agrobacterium cells were grown to OD 600 of 0.8, pelleted, and resuspended in infiltration buffer (10 mM methylester sulfonate, 10 mM MgCl 2 and 150 µM acetosyringone, pH 5.7) and infiltrated into 7-or 8-week-old P. exserta or 5-week-old Nicotiana benthamiana plant leaves using a needleless syringe. The phenotype was observed 2 weeks (Petunia) or 5 days (Nicotiana) after the infiltration. Stable transgenic P. axillaris N lines were generated by leaf disc transformation with Agrobacterium tumefaciens strain GV3101 following an adapted protocol 23 based on Conner et al. 78 .

CRISPR-Cas9 editing
was ordered from Vector-Builder (https://en.vectorbuilder.com/). gRNA 1: CTCACGTCCACTAG GAGATG, gRNA 2: AGGGAGGGACAGCAGAACAT. The construct was transformed into Agrobacterium tumefaciens strain LBA4404. IL5-Ax CRISPR-Cas9 line was generated with the Agrobacterium by leaf disc transformation following the same protocol for stable transgenic P. axillaris N lines. Editing in the T0 generation was detected by PCR performed on genomic DNA with primers targeting the gRNA sites (Supplementary Table 11) followed by Sanger sequencing.

Protein extraction, chitinase and lysozyme activity measurement
Leaf apoplastic proteins were obtained by vacuum infiltration of 7-week-old Petunia plant leaves (the stage when IL5-Ax plants start to show necrotic symptoms) with activity buffer from Chitinase Assay Kit (abbexa, catalogue: abx298854) or Lysozyme Activity Assay Kit (Fluorometric) (Abcam, catalogue: ab211113). The leaves were dried on a paper towel. Then leaf discs (Ø = 9 mm) were obtained by a hole puncher so that the total area of leaves used for protein extraction could be determined. Twenty-two leaf discs from three plants were pooled as one biological replicate. Afterward, the pooled leaf discs were placed in a 50 ml Falcon tube and spun at 1,000 g for 5 min at 4 °C. Three biological replicates (3 × 22 leaf discs from nine plants) were set for each genotype tested.
Chitinase activity was detected with a Chitinase Assay Kit (abbexa, code: abx298854) following its standard protocol. One unit of chitinase activity was defined as the amount of enzyme required to produce 1 µg of N-acetylglucosamine per hour at 37 °C. Lysozyme activity was detected with Lysozyme Activity Assay Kit (Fluorometric) (Abcam, code: ab211113) following its standard protocol. One unit of lysozyme activity was defined as the amount of enzyme that generates 1.0 µmol of 4-MU per minute at pH 5.0 at 37 °C.

Sterile culture
Transparent 5-litre beakers were used as containers for the sterile culture. The beakers with 1.5 litre commercial soil (70% Klasman substrate, 15% Seramis clay granules and 15% quartz sand) were sealed with aluminium foil and autoclaved at 121 °C for 20 min. Then 500 ml autoclaved water was added to the beaker in a sterile environment under laminar flow. Seeds with the genotype of IL5-Ax and IL5-Ex were sterilized with 1% bleach for 10 min and rinsed in sterile water five times. Then the seeds were sowed in the beakers in a sterile environment. The beakers were sealed with several layers of transparent plastic film and put in the growth chamber. The phenotypes of the leaves were analysed 12 weeks after sowing. Four biological replicates for each genotype were observed with similar results.

Conserved cis-elements screening
One-kilobase promoter sequences (1 kb upstream of translation start codons) of the DEGs were extracted and subjected to the motif discovery by MEME-ChIP 79 (http://meme-suite.org/tools/meme-chip) with default parameters.

Genotyping
Genotypes of ChiA1 in the wild accessions of P. axillaris and P. exserta were determined by the CAPS marker designed on the nonsense mutation site. The G-to-T mutation in P. exserta allele of ChiA1 impairs the FokI (GGATGN7) restriction enzyme digestion site. CAPS marker primer sequences can be found in Supplementary Table 11. Genotypes of MYB-FL were determined as previously described 21 .

Microsynteny analysis
Sequences of MYB-FL, ChiA1 and its four nearby genes in the genome were BLASTed against the ITAG release 4.0 cDNA library for S. lycopersicum and PGSC DM4.03 cDNA library for S. tuberosum. For detailed gene IDs and coordinates of the genes shown in the microsynteny analysis, see Supplementary Table 8.

Statistical analysis
GraphPad Prism v.6.0.7 and Microsoft Excel 2016 were used for the statistical analyses (one-way analysis of variance (ANOVA) and two-sided Student's t-test).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
BSR-seq reads have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA708139. Shallow whole-genome sequencing reads data have been deposited under BioProject PRJNA705072. RNA-seq reads have been deposited under BioProject PRJNA705649. The P. axillaris N 4.03 genome assembly has been deposited at NCBI GenBank under the accession JANRMM000000000 (https://www.ncbi. nlm.nih.gov/bioproject/?term=JANRMM000000000).

Code availability
All scripts used in this paper have been deposited on Github: https:// github.com/Kuhlemeier-lab/Petunia_hybrid_necrosis. Major computations were performed on a SLURM computing cluster (UBELIX http:// www.id.unibe.ch/hpc, the HPC cluster at the University of Bern).
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.  Table 1). Red color indicates the genotype of the sample at the specified SNP is homozygous for P. exserta; yellow color indicates the genotype of a sample at the specified SNP is homozygous for P. axillaris N. The black frame indicates the region associated with the HN phenotype in the fine mapping.