The genome and lifestage-specific transcriptomes of a plant-parasitic nematode and its host reveal susceptibility genes involved in trans-kingdom synthesis of vitamin B5

Sequencing and assembly of the H. schachtii genome

We measured (Supplemental Fig. 1), sequenced (BioProject PRJNA722882), and assembled the genome of H. schachtii (population Bonn) using a combination of flow cytometry, Pacific Biosciences sequencing, and Illumina sequencing. H. schachtii has the largest genome (160–170 Mb) of any cyst nematode measured/sequenced to date (Supplementary Table 1). It was sequenced to 192-fold coverage using Pacific Biosciences sequencing (fragment n50 of 16 kb), and 144-fold coverage using Illumina sequencing (150 bp Paired-end reads). The final, polished, contamination-free (Supplemental Fig. 2), assembly (v1.2) included ~179 Mbp contained within 395 scaffolds: 90% of the sequence is contained on scaffolds longer than 281,463 bp (n = 154). The assembly is a largely complete haploid representation of the diploid genome, as evidenced by core eukaryotic genes being largely present, complete and single copy (CEGMA 93.15% complete with an average of 1.12 copies each, and BUSCO (Eukaryota odb9) 79% complete with 8.2% duplicated—Supplementary Table 2). Over three million variants were phased into haplotypes (2029 blocks, N50 239.5 kb, covering 94.7% of the reference) which can be used to predict true protein variants (Supplementary data 1), and 601 larger structural variants were identified (Supplementary data 2).

The trans-kingdom, lifestage-specific, transcriptomes of H. schachtii and A. thaliana provide a holistic view of parasitism

We devised a sampling procedure to cover all major life stages/transitions of the parasitic life cycle to generate a simultaneous, chronological, and comprehensive picture of nematode gene expression, and infection-site-specific plant gene expression patterns. We sampled cysts and pre-infective second-stage juveniles (J2s), as well as infected segments of A. thaliana root and uninfected adjacent control segments of root at 10 hours post infection (hpi – migratory J2s, pre-establishment of the feeding site), 48 hpi (post establishment of the feeding site), 12 days post infection females (dpi – virgin), 12 dpi males (differentiated, pre-emergence, most if not all stopped feeding), and 24 dpi females (post mating), each in biological triplicate (Fig. 1A). We generated approximately nine billion pairs of 150 bp strand-specific RNAseq reads (Supplementary data 3) covering each stage in biological triplicate (for the parasite and the host): in the early stages of infection we generated over 400 million reads per replicate, to provide sufficient coverage of each kingdom.

Strand-specific RNAseq reads originating from host and parasite were deconvoluted by mapping to their respective genome assemblies (H. schachtii v.1.2 and TAIR10). For the parasite, ~500 million Illumina RNAseq read pairs uniquely mapping to the H. schachtii genome were used to generate a set of 26,739 gene annotations (32,624 transcripts – detailed further in the next section), ~77% of which have good evidence of transcription in at least one lifestage (≥10 reads in at least one rep). Similarly for the host, ~2.8 billion Illumina RNAseq read pairs uniquely mapping to the A. thaliana genome show that ~77% of the 32,548 gene models have good evidence of transcription in at least one stage (≥10 reads in at least one rep, even though we only sampled roots). A principal component analysis of the host and parasite gene expression data offers several insights into the parasitic process. Principle component 1 (60% of the variance) and 2 (19% of the variance) of the parasite recapitulate the life cycle in PCA space (Fig. 1B). The 12 dpi female transcriptome is more similar to the 24 dpi female transcriptome than to the 12 dpi male transcriptome. Principle components 1 (75% of the variance) and 2 (10% of the variance) of the host show that the greatest difference between infected and uninfected plant tissue is at the early time points (10 hpi), and that the transcriptomes of infected and uninfected plant material converge over time, possibly due to systemic effects of infection. A 12 dpi male syncytium transcriptome is roughly intermediate between a control root transcriptome and a 12 dpi female syncytium transcriptome. Given that at this stage most if not all of the males will have ceased feeding, this could be due to inadequate formation of the feeding site, or regression of the tissue. In any case, by comparing both principal component analyses, we can see that what is a relatively small difference in the transcriptomes of the feeding sites of males and females is amplified to a relatively large difference in the transcriptomes of the males and females themselves (Fig. 1B).

The consequences, and possible causes, of large-scale segmental duplication in the Heterodera lineage

To understand the evolutionary origin(s) of the relatively large number of genes in H. schachtii in particular, and Heterodera spp. in general, we analysed the abundance and categories of gene duplication in the predicted exome. Compared to a related cyst nematode, Globodera pallida (derived using comparable methodology and of comparable contiguity) the exomes of H. schachtii and H. glycines are characterised by a relatively smaller proportion of single-copy genes (as classified by MCSanX toolkit¹⁷, and a relatively greater proportion of segmental duplications (at least five co-linear genes with no >25 genes between them), with relatively similar proportions of dispersed duplications (two similar genes with >20 other genes between them), proximal duplications (two similar genes with <20 other genes between them), and tandem duplications (two similar genes that are adjacent) (Fig. 2A). Genes classified as segmentally duplicated are clustered into islands (Fig. 2B) of varying gene duplication depth from 2 to 131. Setting an arbitrary duplication depth threshold of 10, approximately half of all segmentally duplicated genes (2273/4547) are grouped into 58 islands across the genome.

To understand when this large-scale segmental duplication occurred, and what genetic capital it has acted on, we cross-referenced the gene duplication analyses with a recent orthologous gene clustering of 61 species across the phylum Nematoda¹⁸. Interestingly, most (434/833) of the orthologous gene clusters that contain H. schachtii genes classified as segmentally duplicated also contain H. glycines genes classified as segmentally duplicated. We, therefore, conclude many are likely bona fide, and that at least half of the modern-day segmentally duplicated regions of the H. schachtii genome were present in the last common ancestor with H. glycines.

Interestingly, gene families that have arisen in the Heterodera lineage (i.e., contain more than one H. schachtii or H. glycines gene, with no representatives from all 61 other proteomes in the analysis) are numerous (1877 orthogroups, accounting for 2072 genes in H. schachtii) and are also enriched in segmental duplications (28% vs 17% for H. schachtii; Fig. 2C). Finally, H. schachtii specific gene families (i.e., those that only contain more than one H. schachtii gene, with no representatives from all 61 other proteomes in the analysis, including H. glycines) are also preferentially located in regions of the genome that have been segmentally duplicated (26% vs 17% of all genes; Fig. 2C). We reason that some of these segmental duplications must have arisen relatively recently (i.e., since divergence from H. glycines) because although it is not certain whether H. schachtii-specific genes arose in the H. schachtii lineage or were lost in the H. glycines lineage, if these gene families were segmentally duplicated before the divergence of H. schachtii and H. glycines, then H. glycines would need to have independently lost both segments to not be in the orthogroup, and this seems unlikely for all cases. Taken together, we conclude that substantial segmental duplication started in the Heterodera lineage after the split from Globodera, was still active in the H. schachtii lineage after the split from H. glycines, and has acted on old and new genetic capital.

To understand the consequences, and possible causes, of these large-scale segmental duplications, we adopted two complementary approaches: (i) a targeted approach to examine what impact this may have had on genes known to be involved in parasitism, and (ii) a non-targeted approach to examine Gene Ontology (GO) terms that are enriched in segmental duplications and/or gene families that are segmentally duplicated. We therefore additionally annotated the predicted proteome of H. schachtii to identify putatively secreted proteins (2669 genes, 3138 transcripts—Supplementary data 4), putative orthologues of known effectors (38 families, 209 genes, 245 transcripts—Supplementary data 5), genes putatively acquired by horizontal gene transfer (HGTs, 263 genes, 319 transcripts—Supplementary data 6), putative cell wall degrading enzymes (CWDEs, 9 families, 39 genes, 45 transcripts–Supplementary data 5), and PFAM domains (to assign Gene Ontology terms, 11,836 genes, 14,705 transcripts—Supplementary data 4).

In the targeted approach, we find preferential duplication of proper subsets of the secretome and the HGTs that are directly implicated at the plant-nematode interface (i.e., effectors and CWDEs, respectively; Fig. 3A). For example, secreted proteins are not generally more duplicated than the rest of the genes in the genome, but effector genes are—including a substantial proportion in segmental duplications (~1/5th). Similarly, putative HGT events are not generally more duplicated than the rest of the genes in the genome, but the cell wall degrading enzymes are (zero single copy genes; Fig. 3A). In the non-targeted approach, we find that genes in segmental duplications are enriched in a very small set of highly-related GO terms involved in nucleic acid synthesis/manipulation and proteolysis (Fig. 3B). Many of the PFAM domains underlying these enriched annotations are associated with DNA replication/integration of viruses/transposons and, remarkably, they tend to be located towards one or both edges of islands of segmental duplication (Fig. 3C). The highly statistically significant functional enrichment of GO terms involved in DNA integration (adj. p value = 7.46 × 10⁻⁹²), DNA biosynthetic process (adj. p value = 1.53 × 10⁻¹¹), viral DNA genome packing (adj. p value = 3.86 × 10⁻²⁰), and DNA replication (adj. p value = 4.67 × 10⁻⁰⁶), coupled with their tendency to be positioned at the edge of segmental duplications, leads us to hypothesise that these functions may have played a role in the large-scale segmental duplication in this lineage.

Differential expression analysis reveals contrasting evolutionary histories of host and parasite genes deployed at specific times

To understand how host and parasite genes are expressed during the parasitic life cycle, only those RNAseq reads that are uniquely mapped to the host genome, not the parasite genome, and vice versa, were used to identify host and parasite genes that were significantly differentially expressed between conditions (Log2FC > +0.5 or <−0.5 and adj. p value < 0.01, Supplementary data 4 and Supplementary data 7). Under these criteria, 69% (18,575) of the nematode genes, and 59% (19,071) of A. thaliana genes, are differentially regulated during the course of infection. To focus on the local infection-specific host response (i.e., removing changes associated with time/systemic changes), we selected only the subset of the A. thaliana transcriptome, 25% (8351 genes), that is differentially expressed between infected and uninfected tissues at each timepoint (hereafter referred to as differentially expressed genes). The most differentially expressed gene at each stage 10 hpi, 48 hpi, 12 dpi female, 12 dpi male, and 24 dpi female respectively) are potential marker genes for future study (Supplementary Table 3).

A majority of host and parasite differentially expressed genes (81% and 94%, respectively) were assigned to superclusters that describe different stages, or groups of stages, of the infection cycle (29 and 28 superclusters for host and parasite, respectively—e.g., Supplemental Fig. 3). For the host, the three largest clusters are of comparable size (10hpi_48hpi, 12 dpi fem._24dpi, and 10hpi_48hpi_12dpi male), and together account for over 2/5ths of all differentially expressed genes (Supplementary Table 4). For the parasite, the cluster with the most differentially expressed genes has a peak at 12 dpi males (Supplementary Table 5). It contains nearly 1/5th of all clustered differentially expressed genes, and is over twice as large as the next largest cluster (cyst). This is noteworthy because males and cysts have generally been less often considered in previous differential expression analyses in the literature compared to pre-parasitic/biotrophic stages, and yet these two life stages alone account for nearly 1/3rd of all differentially expressed genes across the life cycle. Interestingly, genes in the male-specific cluster are non-randomly distributed across the genome (for example, >31% of genes on scaffold 23 (115/364) are in the male-specific cluster, hypergeometric test, adj. p value = 1.51 × 10⁻²⁵ (Supplementary Data 8)). Male-specific genes are often located in segmentally duplicated islands, even when the scaffold as a whole is not enriched, because 30% of all male-specific genes have sequence similarity to ADI82807.1 (860/2,881), and 72% of genes with sequence similarity to ADI82807.1 are in segmental duplications (2235/3094). This gene family was noted as large (474), and male-expressed, in G. pallida¹⁹. This would be an exceptionally large family of related sequences, even for a parasitic organism, and while DNA transposition may be implicated in its formation, the details of that link are at present unclear.

To analyse the evolutionary histories of clusters of genes that peak at specific life stages, we cross-referenced these data with existing analyses of orthologous gene clusters of plants and nematodes. We focused on a subset of eight plant species (Amborella trichopoda (outgroup), the monocots Hordeum vulgare and Zea mays, the solanaceous Solanum lycopersicum and S. tuberosum, the fabaceous Glycine max, and the brassicaceous Brassica rapa (ssp. rapa) and A. thaliana – GreenPhylDB) and a subset of eight nematodes species (the free-living nematode Caenorhabditis elegans (outgroup), the pine wilt nematode Bursaphelenchus xylophilus, the root-knot nematodes Meloidogyne graminicola and M. hapla, the potato cyst nematodes G. pallida and G. rostochiensis, and H. glycines and H. schachtii¹⁸). Six categories of orthologous clusters were considered for the parasite (nematodes, plant-parasites, root-knot nematodes and cyst nematodes, cyst nematodes, Heterodera, and H. schachtii), and six categories of orthologous clusters were considered for the plant (Magnoliopsida, Mesangiospermae (monocots and dicots), Pentapetalae, Rosids, Brassicaceae, and A. thaliana, Supplemental Fig. 4).

We correlated the putative annotation, orthologue definition, and transcriptional clustering data to explore the relatedness and functions of subsets of differentially expressed genes. Initially, we analysed the six differential expression clusters that describe the transition to biotrophy (J2, J2_10hpi, 10 hpi, 10_48hpi, 48 hpi, and J2_10_48hpi; Fig. 4A). Of these six, five are significantly enriched in genes encoding putatively secreted proteins (hypergeometric test, false discovery rate (FDR) adj. p value < 0.01), the exception being the J2 specific cluster, which also contains zero known effectors. This does not mean that no effectors are expressed in the J2 stage, but rather that none are specific to the pre-infective J2 timepoint: all effectors expressed in the pre-pre-infective J2 timepoint are also expressed at some other timepoint such that they are clustered elsewhere. Importantly, at 10 hpi the nematode is still a motile J2 and, consistent with this, the only cluster significantly enriched in CWMEs is J2_10hpi (hypergeometric test, FDR adj. p value = 5.41 × 10⁻¹⁰) commensurate with the idea that CWMEs facilitate both entry into, but also movement within, the host. While 65% of all putative orthologues of known effectors are contained within five clusters that describe the transition to biotrophy, only one of these (48 hpi) is significantly enriched (hypergeometric test, FDR ad. p value = 1.08 × 10⁻⁴⁹) and contains 38% of all putative orthologues of known effectors. The only other cluster enriched for effectors contains exclusively life stages involved in biotrophic interaction with the host (48hpi_12dpifem_12dpimale_24dpi, hypergeometric test, FDR adj. p value = 1.15 × 10⁻¹⁰).

Cross-referencing the differential expression clusters that describe the transition to biotrophy with the orthogroup analysis suggests that the transition to biotrophy is generally characterised by a proportional decrease in genes shared among nematodes, and a proportional increase in various orthologous genes cluster categories specific to subsets of plant-parasitic nematodes. On a finer scale, how the distribution of genes across the orthologous gene cluster categories compares between differential gene expression clusters that describe the transition to biotrophy is very informative. To summarise the biggest differences: the beginning of the transition to biotrophy is characterised by the greatest proportional increases in plant-parasitic nematode-specific genes and cyst and root-knot nematode-specific genes (J2_10hpi), and the greatest proportional decrease in Heterodera-specific genes (J2, Fig. 4B); the middle of the transition (10 hpi) is characterised by a proportional increase in Heterodera-specific genes; and the end of the transition (10hpi_48hpi and to a greater extent 48 hpi) is characterised by a proportional increase in cyst nematode-specific genes and Heterodera-specific genes. Considering the functional annotation of genes in these clusters (Fig. 4A), we conclude that genes expressed exclusively in the pre-infective J2 are not involved in host manipulation (unless they are also expressed in another timepoint (e.g., J2_10hpi)). Given that the 10 hpi cluster is characterised by a larger proportional increase in Heterodera-specific genes than cyst nematode-specific genes, coupled with the fact that the greatest difference between infected and uninfected plant tissue is at 10 hpi and dominated by immunity-related functions, these data appear consistent with the idea that a greater degree of specialisation from the parasite is needed to overcome the early, immunity-related, host response. Finally, we conclude that H. schachtii is a model for the most agriculturally important cyst nematode, H. glycines, because H. schachtii-specific genes represent a tiny proportion of those expressed during the transition to biotrophy (generally <5%).

Subsequently, we analysed the five differential expression clusters, of the host and of the parasite, that describe discrete stages of infection (10 hpi, 48 hpi, 12dpi_male, 12dpi_female, and 24dpi_female; Fig. 5). Strikingly, host and parasite genes upregulated at the same time of infection have different, and contrasting, evolutionary histories within their respective phyla. Parasite genes in expression clusters that peak at individual times of infection are generally less well conserved across the phylum Nematoda when compared to all differentially expressed genes in the experiment (Fig. 5 and Supplemental Fig. 5). In contrast, host genes in expression clusters peak at those same times of infection are generally better conserved across the kingdom Plantae when compared to all differentially expressed genes in the experiment (Fig. 5 and Supplemental Fig. 5). Although these evolutionary time scales are not directly comparable, these data may suggest that infection-specific gene expression is characterised by lineage-specific nematode genes (including, for example, effectors) modulating the expression of widely conserved plant genes (including, for example, genes involved in the development of the feeding site).

Congruent differential expression of metabolic pathways during infection highlights nematode susceptibility genes

To understand the roles of widely conserved host genes, and lineage-specific parasite genes, that are differentially regulated during biotrophy (Fig. 5), we used the KEGG Automatic Annotation Server Ver. 2.1 (https://www.genome.jp/kaas-bin/kaas_main) to annotate metabolic pathways in the host (3,860 KO terms on 9,240 genes) and the parasite (3479 KO terms on 4990 genes (Supplementary Data 9)). Given the partial overlap of KO terms (1,998 in common between host and parasite) we examined pathways that were incomplete in the parasite, but complemented by the host.

The vitamin B5 (pantothenate) biosynthetic pathway (M00119) is an example of congruent differential expression of metabolism involving widely conserved host genes and lineage-specific parasite genes. The complete pathway is conserved in most plants, and all steps (even when more than one gene encodes a particular function) are upregulated during infection of A. thaliana with H. schachtii (Fig. 6A), with the exception of the last step pantothenate synthetase (AtPANC). The dominant peak of expression in all upregulated cases is at 12 dpi in the syncytia associated with the female. The expression of this pathway in syncytia associated with the 12 dpi male is either indistinguishable from the control (e.g., 2.2.1.6), or significantly upregulated compared to the control, but much less highly upregulated than in syncytia associated with the 12 dpi female and approximately at the same level as in the syncytia associated with the 24 dpi female. While all animals require vitamin B5 in order to complete the basal metabolism of proteins, carbohydrates, and lipids, the biosynthesis pathway is considered absent (vitamin B5 is an essential nutrient for animals). Unusually for an animal, the H. schachtii genome encodes two genes (Hsc_gene_23032 and Hsc_gene_23033) annotated to carry out only the last step (PANC) of the vitamin B5 biosynthesis pathway²⁰. These putative H. schachtii PANC-encoding genes (Hs-panc) are dissimilar in sequence to plant orthologues. Querying the non-redundant nucleotide archive, the most similar sequences to Hs-panc are from a few animals with obligate interactions with plants (primarily, although not exclusively, plant-parasitic nematodes) and the vast majority from bacteria (Fig. 6C). Phylogenetic analyses group the Hs-panc, along with many other plant-parasitic panc-like sequences, towards the base of a well-supported clade (100 bootstrap) containing sequences primarily from actinobacteria (the most closely clustered sequence being from Candidatus Rhodoluna planktonica, WP_070954084.1 [https://www.ncbi.nlm.nih.gov/protein/WP_070954084.1?report=genpept], a well characterised but as-yet uncultured organism). Importantly, both Hs-panc encode several well-supported introns and are in the middle of a large 1.3 Mb scaffold surrounded by classical nematode genes, disproving the possibility of contamination (Fig. 6B). The most parsimonious explanation for the existence of panc in the H. schachtii genome is, therefore, their acquisition by horizontal gene transfer from bacteria to the last common ancestor of the suborder Hoplolaimina—indicative of an important role in plant-nematode interactions across the group – and subsequent tandem duplication in H. schachtii or its progenitor. Remarkably, Hs-panc complement the upregulation of the first steps in the pathway in the host (and the absence of upregulation of the last step in the pathway in the host): Hs-pancs are upregulated in the 12 dpi female, to a lesser extent in the 12 dpi male, and to a greater extent in the 24 dpi female (Fig. 6A).

Taken together, we conclude that the vitamin B5 (pantothenate) biosynthetic pathway is congruently differentially expressed between kingdoms in this interaction. The presence of panc-like sequences in two other groups of animals outside the plant-parasitic nematodes (Tetranychus urticae – a species of plant-feeding mite generally considered to be a pest—XM_015927639.1²¹, and Apis mellifera—the western honey-bee—XM_001123024.1 is intriguing given their similarly obligate relationship with plants (Supplemental Fig. 6). In the case of A. mellifera, it is somewhat more ambiguous than for the plant-parasitic nematodes and the mite. The predicted protein sequences of the A. mellifera PANC are 65% identical with the most similar bacterial sequence in NCBI (with much less similar nucleotide sequences), the gene does not contain introns, however, and it is on a short scaffold. Other examples of congruently differentially expressed metabolic pathways that warrant further investigation are more complicated (due to not being exclusively linear pathways) but show similar trends (e.g. biotin (Supplemental Fig. 7 and Supplementary Tables 6 and 7 for parasite and host respectively), riboflavin, stilbenoid, and vitamin B6).

-These data support the hypothesis that the plant-encoded enzymes of the first part of the vitamin B5 biosynthesis pathway support nematode infection, and may thereby function as nematode susceptibility genes. To investigate whether AtPANB1 and AtPANB2 (the penultimate step in the pathway, and the last step upregulated upon infection) plays a role in cyst nematode parasitism, we ordered T-DNA insertion lines for both AtPANB1 and AtPANB2. After several unsuccessful attempts to generate a homozygous loss-of-function for AtPANB2, we focused our further analysis on AtPANB1 (atpanb1-1; Supplemental Fig. 8). With the exception of a slight delay in flowering time, atpanb1-1 was phenotypically indistinguishable from wild-type in the absence of infection: hypocotyl length, root and shoot fresh weight, and siliques sizes, were not altered compared to Col-0 (Supplemental Fig. 9). No morphological differences were observed in roots between Col-0 and atpanb1-1. During infection, however, the following marked differences were noted between the mutant and wild-type lines. We found that atpanb1-1 supported fewer females at 14 dpi (Fig. 7A; p value < 0.0001). Moreover, even those nematodes that were able to infect the mutant are less fit as evidenced by the formation of smaller syncytia (Fig. 7B, D, p value < 0.0001), smaller females (Fig. 7C, D, p value = 0.02), smaller cysts (Fig. 7E; p value = 0.049), and lower average number of eggs per cyst (Fig. 7F; p value = 0.049). We have not noted striking differences in egg morphology, size or hatching rate. However, we have not quantified these parameters. Importantly, the reduced susceptibility of atpanb1-1 to nematodes can be rescued by overexpressing AtPANB1 under the control of a 35 S promoter in the atpanb1-1 background (Fig. 8). These experiments were performed multiple times (at least three) independently, with the same outcome (Data for additional replicates are available in the Source Data file)

Next, we conducted microscopic examinations of anatomical and ultrastructural features of syncytia induced in roots of atpanb1-1 plants and compared them with syncytia induced in wild-type Col-0 plants (Fig. 7G–J and Supplemental Fig. 10). In both genotypes, syncytia were induced in the vascular cylinder and they were composed of vascular cylinder cells only (Fig. 7G–J)⁹. In wild-type Col-0, syncytia were surrounded by a continuous layer of periderm (Fig. 7G; Supplemental Fig. 10B, D)²², excluding regions next to nematode heads (Supplemental Fig. 10A, C). However, in the case of syncytia induced in atpanb1-1 roots, the periderm layer was not yet differentiated next to the nematode head (Supplemental Fig. 10E) and syncytia developed poorly at 5 dpi (Fig. 7H and Supplemental Fig. 10F). Overall, syncytia induced in atpanb1-1 were less well developed than in Col-0 plants as indicated by reduced size and number of cells incorporated into syncytia. This was especially evident at 5 dpi (Fig. 7G, H; Supplemental Fig. 10A, B, E, F). At 10 dpi, syncytia induced in atpanb1 roots were still smaller on cross sections than syncytia induced in wild-type plants, but this difference was less apparent (Fig. 7I, J and Supplemental Fig. 10D, H). At the ultrastructural level, the organisation of syncytial protoplasts in Col-0 plants was typical for syncytia induced in Arabidopsis roots. Syncytial cytoplasm is marked by a high density of organelles (Supplemental Fig. 10I, J)^{9, 10}. A contrasting situation was found in syncytia induced in atpanb1-1 roots (Supplemental Fig. 10K, L). At the widest parts of syncytia, some distance from the nematode heads, the syncytial cytoplasm exhibits large regions almost free from plastids, mitochondria and structures of endoplasmic reticulum both at 5 and 10 dpi. Syncytia are metabolically hyperactive sinks from which nematodes feed throughout their life cycle. The ultrastructural features developed in atpanb1-1 are unusual, and might indicate that syncytia induced in these plants are metabolically less active and therefore less efficient in supporting nematode development leading to/stemming from the formation of smaller and undernourished females. These data suggest that a reduction in pantoate supply in atpanb1-1 mutant compromises nematodes’ ability to synthesise vitamin B5, leading to smaller females, which are unable to fully support the establishment and/or maintenance of a syncytium. It is therefore likely that the syncytial phenotype observed is an indirect effect caused by the nematode’s general lack of fitness.

To investigate whether AtPANB1 plays a similar role in other plant-interacting organisms, we analysed changes in susceptibility of atpanb1-1 towards the necrotrophic fungus Botrytis cinerea, and the beneficial endophyte Serendipita indica. In both cases, no differences were detected between atpanb1-1 and the wild-type control (Supplemental Fig. 11A, B). Next, we investigated whether the decrease in susceptibility of atpanb1-1 to cyst nematodes is due to changes in classical immune responses. To this end, we tested ROS production, a hallmark of basal defence in plants, upon treatment of an immunogenic peptide, flg22. We did not observe any changes in ROS burst, both over time or cumulative (Supplemental Fig. 11C). Taken together, these data suggest that AtPANB1, the last step in the pathway upregulated during infection, is a specific non-immunity-related nematode susceptibility gene. Given that this pathway is linear, and that all previous steps are similarly upregulated, the implication is a series of nematode susceptibility genes for future study. The wide conservation of panc in plant-parasitic nematodes (Fig. 6B), and where we have data on their similar upregulation during parasitism (e.g., G. rostochiensis, Gr-panc—GROS_g05752—Supplemental Fig. 12), extends this intriguing possibility to at least the cyst nematodes, but perhaps most other plant-parasitic nematodes of global agricultural importance²³.

Balancing the cost vs reward of any S gene is an important consideration in the application of the technology. Generally, one would aim to achieve a partial knockout of function such that infection is impaired, but the useful biology of the plant is not impaired. Intuitively one might think that the gene with the highest expression is the most suitable for knockout—in the case of S genes, weak phenotypes (with no or minimal impact on “normal” plant physiology, but large impacts on parasitic physiology, as shown herein) are ideal. Redundancy in the pathway (as is the case for AtPANB) can help in this regard, and is anticipated to be particularly relevant to plant-parasitic nematodes given the extreme metabolic strain and relative fragility of the feeding site. With the exception of a slight delay in flowering time, atpanb1-1 was phenotypically indistinguishable from wild-type in the absence of infection. Generating a CRISPR-based allelic series of panb mutants in one of the major crop plants infected by cyst nematodes, for example, potato, would be an excellent use case because: the pathogen is closely related; the crop is polyploid (and so a wide variety of reduction of function mutants can be generated); flowering time is not required for crop production (and may not be similarly affected); and finally the corresponding homologue of HsPANC (GROS_g05752) is similarly and almost exclusively upregulated during biotrophy (14 days post-infection, Supplemental Fig. 12).

Trans-kingdom compartmentalisation of the vitamin B5 pathway

In contrast to AtPANB, AtPANC is neither upregulated in syncytium nor plays a role in cyst nematode infection: loss-of-function homozygous mutants for AtPANC (Fig. 6; Supplemental Figs. 6 and 7) are indistinguishable from wild-type plants during infection (Fig. 9A–C). Complementing this lack of upregulation, the expression of both Hs-panc1, but particularly Hs-panc2, is not only increased after onset of parasitism (48 hpi) but remains high as the nematode proceeds through the rest of its life cycle, pointing to a functional role for Hs-panc, in particular Hs-panc2, in parasitism (Fig. 6A). To test this hypothesis experimentally, we used in vitro RNAi targeting the more highly expressed paralogue (by a factor of nearly 10) Hs-panc2 (Supplemental Fig. 13), which caused a significant decrease in expression of Hs-panc2 (Fig. 9D) and resulted in fewer female nematodes infecting wild-type Col-0 (Fig. 9E). The average sizes of the females and syncytia were also significantly reduced on Col-0 infected with nematodes showing a reduced expression of Hs-panc (Fig. 9F, G). Inconsistent effects were seen on the number of males (Supplemental Fig. 14). Pre-infective juveniles of cyst nematodes (pre-J2s) are non-feeding, essentially arrested in development, and use their lipid reserves to locate and invade the host plants. They start feeding only after their feeding sites are established. It is therefore highly unlikely that knocking down HsPANC affects nematode behaviour pre-infection (indeed, HsPANC is lowly expressed when the nematode is outside the plant). Decreased expression of HsPANC is hypothesised to compromise the fitness of nematodes post infection, which will in turn affect nematode’s ability to establish and/or maintain their feeding sites, and therefore manifest in the decreased infection measured herein.

Importantly, Hs-panc functions in the nematode. Hs-panc does not encode an N-terminal secretion signal, and in situ hybridisation of digoxigenin-labelled probes to Hs-panc mRNA reveals expression throughout the nematode body with a particularly strong signal in the intestinal region of the nematode (Fig. 9H). Taken together, we conclude that the contribution of Hs-panc to parasitism success is due to its role in the nematode, during parasitism.

To understand how Hs-panc contributes to parasitism success in the nematode, we tested whether Hs-pancs encode a functional enzyme. Hs-panc1 cDNA was cloned into PUC18, and was tested for its ability to complement panc mutant from Escherichia coli (AT1371). The E. coli panc gene served as a positive control. All strains grew well in the presence of pantothenate in the growth medium. In absence of exogenous pantothenate, the panc mutant containing the empty vector did not grow. However, the strains transformed with either bacterial or nematode panc grew well with no significant differences detected between strains (Fig. 10). We, therefore, concluded that Hs-panc encodes a bona fide PANC.

The presence of panc homologues in the cyst nematode’s genome indicates an evolutionary benefit, but the molecular fundamentals of this peculiar separation remain elusive. Previous studies show that plant PANC is subjected to uncompetitive substrate inhibition at higher pantoate concentrations²⁴. Similarly, the activity of plant PANC is also inhibited by products, and in particular by the accumulation of pantothenate²⁵. E. coli, on the other hand, produces vitamin B5 in large amounts, and any excess is secreted into the medium, strongly suggesting that vitamin B5 biosynthesis, including PANC activity, is not very tightly regulated. In comparison to bacteria, plants are unlikely to depend on continuous excretion of vitamin B5, as this would require an elaborate transport system²⁶. The biosynthesis of vitamin B5 is therefore expected to be more tightly regulated in plants. In line with this idea, PANC from plants e.g., Lotus japonicus is shown to have strong substrate inhibition by pantoate, whereas enzymes from E. coli and Mycobacterium tuberculosis follow hyperbolic kinetics^{24, 27, 28}. Similarly, a comparative analysis of kinetic properties showed that although the two-step reaction mechanism is conserved between the AtPANC and EcPANC, AtPANC exhibits large allosteric effects, whereas EcPANC displays essentially non-allosteric behaviour²⁶. Based on these observations, we propose a hypothesis (Supplemental Fig. 15) that compartmentalisation of the vitamin B5 pathway between plants and nematodes helps to avoid feedback/feed-forward inhibition and ensures a consistent supply of vitamin B5 to rapidly developing nematodes. Clarifying further details of the nematode vitamin metabolism, including details of the activity of HsPANC and AtPANC may aid in the development of nematode-resistant germplasm through the targeted deletion of susceptibility genes.

To conclude, we have shown that H. schachtii has the largest genome, and second-largest exome, of any cyst nematode analysed to date. Substantial segmental duplication contributed to the large gene number, it started at some point in the Heterodera lineage after the split from Globodera, was still active in the H. schachtii lineage after the split from H. glycines, and has acted on old and new genetic capital. A reference trans-kingdom, lifestage-specific, transcriptome of H. schachtii and A. thaliana reveals that host and parasite genes upregulated at the same time of infection have different, and contrasting, evolutionary histories within their respective phyla. Nevertheless, congruent differential expression analysis of metabolic pathways in this time course transcriptome highlights nematode susceptibility genes involved in vitamin B5 biosynthesis. Taken together, we conclude that the vitamin b5 biosynthesis pathway of the plant-nematode hologenome starts in the plant, ends in the nematode, and is putatively enabled by a horizontal gene transfer from a bacterium.

Source: Ecology - nature.com