A host‐free transcriptome for haustoriogenesis in Cuscuta campestris: Signature gene expression identifies markers of successive development stages

Abstract The development of the infection organ of the parasitic angiosperm genus Cuscuta is a dynamic process that is normally obscured from view as it happens endophytically in its host. We artificially induced haustoriogenesis in Cuscuta campestris by far‐red light to define specific morphologically different stages and analyze their transcriptional patterns. This information enabled us to extract sets of high‐confidence housekeeping and marker genes for the different stages, validated in a natural infection setting on a compatible host. This study provides a framework for more reproducible investigations of haustoriogenesis and the processes governing host–parasite interactions in shoot parasites, with C. campestris as a model species.

. All are root-and leafless, and either completely lacks the ability to perform photosynthesis or show minimal photosynthetic activity of unsustainable levels (van der Kooij et al., 2000;Vurro et al., 2017). These parasites are therefore considered to be obligate stem parasites.
Cuscuta has a thread-like stem that grows while rotating in a counter-clockwise motion. As soon as it encounters a host stem, it twines around it and effectively develops lateral parasitizing organs termed haustoria on the side of the stem facing the host (Kokla & Melnyk, 2018). Early stages of haustoria promote the attachment of the parasitic shoot to the host surface. This is mediated by the production of sticky substances that allow the parasite to adhere to the host (Galloway et al., 2020;Vaughn, 2002) and coincides with a swelling of the Cuscuta stem and a re-shaping of the epidermal cells into club-shaped cells. During a successful infection on a compatible host plant, the haustorium evolves endophytic structures that penetrate the host surface, growing into the host (Johnsen et al., 2015), and enable the parasite to sequester water, inorganic salts and organic compounds, among them RNAs, proteins, hormones and metabolites (Kim & Westwood, 2015). Specialized so-called feeding hyphae at the tip and the sides of this mature haustorium serve to connect to different cell types of the host (Shimizu & Aoki, 2019;Vaughn, 2003Vaughn, , 2006. Host molecules have been shown to particularly influence the formation of these feeding connections between the endophytic mature haustorium and the host (Liu et al., 2020;Narukawa et al., 2021), so that this organ with its unique cell specializations is a product of internal (i.e., Cuscuta-derived) and external (i.e., host-derived) regulating factors.
Attacks by Cuscuta parasites appear to go unnoticed by susceptible hosts, ostensibly because their surface chemistry and architecture are quite similar. Attempts to decipher action and reaction in the host-Cuscuta interaction have proven difficult, mainly since most of the infection organ is hidden from view. Early penetrating stages of the haustorium cannot be visually distinguished from feeding mature haustoria without dissecting and analyzing the infection sites. This, however, may lead to artifacts in gene expression analysis due to wound responses and degradation of RNAs in the cut tissues. In addition, the assignment of the harvested tissue samples to a development stage may differ, depending on the experimenter's judgment and morphological observation skills. The publication of the genome sequence of Cuscuta campestris (Vogel et al., 2018) and Cuscuta australis (Sun et al., 2018), the possibility to study fluorescent fusion proteins in calli (Švubová & Blehová, 2013) and the adhesive disks of the infection organs (Lachner et al., 2020), and the development of an advanced artificial host system providing control of the parasite environment throughout its life cycle (Bernal-Galeano & Westwood, 2021) are just some of the recent milestones in parasitic plant research that helped to ascertain Cuscuta a role as a shoot parasitic model. In order to be able to draw conclusions across datasets from past (Ranjan et al., 2014), present (Jhu et al., 2021;this study) and future transcriptomic studies, it is highly desirable to overcome potential shortcomings in morphological sample assignment and to identify robust molecular markers that are highly specific for the different developmental stages.
Initiation and progression of haustorium formation seem to rely on several signals that are not necessarily host-dependent and can be conveniently used to stimulate haustorium development without a living host. As Bernal-Galeano and Westwood (2021) have recently shown, Cuscuta can even be grown entirely on an artificial feeding support system. The possibility to replace one highly variable factor (the host) can provide points of reference for future studies and may help to decipher how different aspects of haustorium development are controlled. The most commonly used way of host-free haustorium induction is a combination of far-red (FR) light and tactile stimuli (Lachner et al., 2020;Olsen et al., 2016;Tada et al., 1996), although also blue light was shown to induce haustoria (Haidar, 2003;Kaga et al., 2020). Host-free induced haustoria are easy to monitor, and they exhibit a more uniform and predictable development compared to host-induced haustoria. In the present study, FR induction was employed to produce non-hostinduced haustoria of C. campestris, provide easy-to-identify hallmarks of early haustorium development and correlate gene expression patterns with characteristic morphological traits during the respective stages. Genes whose expression changed with the transition from one stage to the next were identified, and their use as haustorial stagedefining markers was verified using independently collected on-host samples. We present a list of robust marker genes that can be used to supplement or substitute visual stage assignment. With this study, we thus provide a framework for more reproducible investigations of haustoriogenesis and the processes governing host-parasite interactions in shoot parasites, with C. campestris as a model species.

| Plant material
Plants were maintained in a greenhouse at the Phytotron of the Arctic University of Norway, Tromsø, under 24 h daylight and approximately 21 C. Cuscuta campestris was originally obtained from the Botanical Garden of the University of Kiel (Germany) and propagated on Pelargonium zonale as compatible host. Solanum lycopersicum cv. M82 was grown in sphagnum peat (Veksttorv, Tjerbo, Norway) mixed at a 2:1 (v/v) ratio with perlite (Agra-perlite, PULL Rhenen, The Netherlands).
Solanum lycopersicum cv. M82 is tolerant to C. campestris attacks.

| Host-free haustorium induction
Distal portions of C. campestris shoots, including tips (approximately 10 cm), were harvested from individuals feeding on P. zonale. Host-free induction of haustoriogenesis was carried out as described by Olsen et al. (2016) by placing them between two reversed plastic Petri dish halves (ø = 13.5 cm). Gentle pressure was applied, and halves were taped together. Shoots were placed in an upright position with their bottom into water, and irradiated with FR light (740 nm) for 2 h. They were then kept in the dark for 6 days. Shoots were finally inspected under a SteREO Lumar V12 (Zeiss) microscope equipped with an AxioCam MRc5 (Zeiss) camera. Consecutive stages of haustorium development were defined based on morphological characteristics ( Figure 1A-H). Ten visually similar sites were cut from as many FR-induced stems as needed and pooled together before frozen in liquid nitrogen. Since each FR-treated stem contained variable numbers of haustoria in each stage, the number of stems represented in each pool differed but was at least three stems. Further, stem sections below and above areas with developing infection sites were collected using the same routine. Three biological replicates were harvested for each sample type and used for sequencing.

| Read mapping and quantification
Quality assessment of the read sequences was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (I) Representation of Venn intersections relative to niS. Numbers are overlapping differentially accumulated transcripts between pairwise comparisons of sample types. Only transcripts that have an FDR-corrected P-value ≤0.05 and a jlog 2 (FC)j ≥1.5 were retained. (J) Volcano plots for differentially accumulated transcripts between infective (SWE, ATT and PEN) and non-infective (niS) tissues. The x-axis shows the log 2 fold change (FC) in transcript accumulation between sample types, and the y-axis shows the statistical significance of the differences. Horizontal and vertical dashed lines are P-value ≤0.05 and jlog 2 (FC)j ≥1.5 thresholds, respectively. Up-and downregulated transcripts that meet both criteria are highlighted in red. ar, attachment ring with large expanded cells; c, cortex; dec, digitate epidermal cells; e, epidermis; ecz, elongated cell zone of the haustorial primordium; hbc, haustorial body cells; hi, haustorial initials; htc, hyphae-like digitate haustorial tip cells; mcz, meristem cell zone of the haustorial primordium; p, pith; v, vascular tissue ). Quality trimming, adapter clipping and mapping of processed reads on the coding sequences derived from the reference genome of C. campestris (v. r0.31) (Vogel et al., 2018) were performed using CLC Genomics Workbench (v. 11) (Qiagen Bioinformatics). For trimming, a base-calling error probability of 0.05, a maximum of 2 ambiguous nucleotides and a minimum length of 15 nucleotides were set as thresholds. The build-in CLC read mapper was set to a mismatch cost of 2, an insertion/deletion cost of 3, a length fraction of 0.8 and a similarity fraction of 0.8. Raw read counts per gene were normalized to reads per kilobase million (RPKM).

| Differential expression analysis
Differentially expressed gene transcripts were identified between sample groups using the Bioconductor software package EdgeR (v. 3.28.1) (Robinson et al., 2010). Raw read counts were used as input data, following the user's guide recommendation. Gene transcripts that had (1) a minimum of one read count in a "worthwhile" number of samples (which is defined by the software) and (2) a total of at least three counts per million across all samples were considered expressed and kept for analysis. Normalization factors were obtained using the trimmed mean of M values (TMM) method. Pairwise comparisons were based on an exact test, using the quantile-adjusted conditional maximum likelihood (qCML) method and allowing both common and tagwise dispersion approaches. P-values were adjusted using the Benjamini-Hochberg (BH) procedure (Benjamini & Hochberg, 1995).
A false discovery rate (FDR) cut-off of 0.05 and a minimum log 2 foldchange of 1.5 were used to retain transcripts of interest.

| Hierarchical clustering
Two-dimensional hierarchical clustering was performed on transcripts that were differentially regulated in at least one of any pairwise comparison of stages. Counts per million (CPM) were calculated from raw read counts, log-transformed and converted into z-score values. Clustering of both samples and transcripts was obtained by applying the complete-linkage method to the respective Euclidean distance matrices.

| Correlation of expression pattern to stages
Low count (noisy) transcripts were discarded by applying a minimum mean RPKM cut-off of 2 across all samples. For each transcript, gene significance (GS) was assessed after Langfelder and Horvath (2008).
Pearson correlation between transcript accumulation and simulated stage-specific expression patterns (with 1 indicating expression in the biological samples of one stage and 0 indicating no expression in the others), was calculated. Transcripts of interest were defined as those having a minimum GS of 0.7 (on a scale of 0-1) and a significant Pvalue (P ≤ 0.05).

| Gene annotation and enrichment analysis
Assignment of C. campestris transcripts to MapMan4 (v.3.0) functional categories was performed using the online Mercator annotation tool (Schwacke et al., 2019). Gene set enrichment analyses were performed by applying a hypergeometric test. P-values were adjusted using the BH procedure. Bins with an FDR ≤0.05 were considered significantly enriched.

| Definition of housekeeper candidates
The following criteria were adopted: (1) expression in all samples; (2) a minimum mean log 2 (RPKM) of 5; (3) low variance, defined as lying within the 10% lower coefficient of variation (CV) and median absolute variation (MAD); (4) no significant fold change, defined as a log 2 (FC) <1.5 and an FDR >0.05 in any of the all versus all pairwise comparisons of host-free haustorium development stages.

| Definition of marker candidates
Differential expression and estimation of gene significance were carried out independently. By applying a reductionist approach (which ensured both high specific fold-changes in expression and a consistent pattern among sample replicates), candidate markers of the different stages of haustorium development were defined as those that (1)

| Host plant parasitization
Distal portions of C. campestris shoots, including tips (approximately 10 cm), were harvested from individuals feeding on P. zonale and attached to S. lycopersicum cv. M82 stems. A continuous 16 h light, 2 h FR light and 6 h dark regime was applied. Coiled regions containing tissues of both parasite and host plants were randomly sampled. Cross-sections were made on both sides of individual coils using razor blades. Interface regions were visually inspected under a SteREO Lumar V12 (Zeiss) microscope equipped with an AxioCam MRc5 (Zeiss) camera and assigned to one of the aforementioned stages of haustorium development. Samples containing parasites in contact with host tissues were frozen in liquid nitrogen. Cross-sectioned non-infective dodder stems (that were exposed to the same light regime) were further sampled in a similar way. Tissue homogenization, total RNA extraction and DNase treatment were performed as described above. RNA concentration and purity were measured with a NanoDrop ND-1000 (Thermo Fisher Scientific) spectrophotometer.
2.13 | Reverse transcription quantitative real-time PCR (RT-qPCR) The expression of selected reference and marker candidates was measured by RT-qPCR in both host-free and host-dependent samples. Gene-specific forward and reverse primer pairs were designed using Primer3 (v. 2.4.0) (Untergasser et al., 2012). Target specificity was checked using the blastn-short tool from the NCBI BLAST+ suite (v. 2.6.0) (Altschul et al., 1990). Both software were used as standalone to enable batch computing. The SuperScript II Reverse Transcriptase (Invitrogen) was used with anchored oligo(dT) 18 primers to reverse transcribe DNase-treated RNA

| Evaluation of housekeeper stability
Three widely used algorithms were applied to rank the expression stability of the housekeeper candidates in both the hostfree and host-induced systems based on their Cq values. geNorm (Vandesompele et al., 2002) was used to calculate the expression stability values M. Candidates with the lowest M-value were the most stably expressed ones. Pairwise variations Vn/Vn + 1 were then used to determine the minimum number of reference genes required for normalization of transcript relative abundances.
NormFinder (Andersen et al., 2004) was used to calculate stability values by taking both intra-and inter-group variation into account, with the lowest values indicating the most appropriate candidates to be used. BestKeeper (Pfaffl et al., 2004)

| RESULTS
3.1 | A host-free transcriptome for early haustorium development When exposed to FR light in the presence of a tactile stimulus, apical portions of Cuscuta stems develop haustoria that in their early stages bear many of the morphological and molecular characteristics of naturally developing haustoria, as shown here for the sequenced species C. campestris ( Figure 1A-H). These morphological signs are easily detectable under low magnification (e.g., with a stereo microscope). The first visual sign appears after one to two days with a slight bump where the otherwise flat epidermal cells appear rounded on the surface ( Figure 1A,C). Cross-sections of these sites show that in comparison to the non-induced stems, these epidermal cells undergo a strong elongation perpendicular to the surface and haustorium initials appear in the region between the cortex and the vascular tissue ( Figure 1B,D). This early infective "swelling stage" (SWE), soon becomes more pronounced and culminates in a macroscopically visible structure of 1-2 mm ( Figure 1E) that begins to stick to the surface it faces (Vaughn, 2002) ("attaching stage" [ATT]). In cross-sections, the newly formed haustorial primordium with its meristem and elongation zones (Lee, 2007(Lee, , 2008 is visible upon staining with the polychromatic stain Toluidine Blue O by its bright blue/turquoise color (indicating the presence of lignin or other polyphenolic compounds in the cell wall) ( Figure 1F). The epidermal cells facing the surface are palisade-formed and dark purple stained cells owing to their high pectin content (Vaughn, 2002). In the center of this lateral structure, the haustorium finally emerges and the adhesive surface forms a ring around it ( Figure 1G,H). The emerging haustorium has digitate cells at its tip that sometimes protrude from the surface ( Figure 1H) and are occasionally termed "search hyphae" (Kaga et al., 2020) although they more likely serve host penetration rather than feeding purposes. On a host, the split in the sticky adhesive ring together with the action of degrading enzymes causes a rupture in the infected tissue and the intrusion of the haustorium that will then begin to feed. In the host-free system used here to induce haustoria, this is the final stage that can be observed and will be referred to as the "penetrating stage" (PEN). FR-exposed stem sections below and above the region where haustoria developed (non-  (Table S1, Appendix S1). Hierarchical clustering analysis of the sequencing data confirmed the consistency between biological replicates ( Figure S1).
3.2 | Host-free haustorium development is accompanied by major transcriptional reprogramming Clusters 7 and 10 grouped DEGs that were mostly upregulated in ATT. Cluster 8 was mainly associated with PEN. The transition from niS to SWE was marked by a shift toward a strong representation of MapMan4 functional bins related to phytohormones action, chromatin and cell cycle organization as well as RNA biosynthesis and processing, protein biosynthesis, modification and homeostasis, cytoskeleton and cell wall organization, and solute transport ( Figure 2C). It should also be noted that almost all the metabolism-related categories contained bins significantly enriched with DEGs. The transition from SWE to ATT showed a rather moderate number of functional categories (including RNA biosynthesis) with few representatives. Finally, the transition from ATT to PEN was associated with a distinct number of DEGs in bins related to phytohormones action, RNA biosynthesis, cell wall organization, solute transport, and enzyme classification.
3.3 | Stage-specific patterns of gene expression during host-free haustorium development allow for the identification of marker genes Gene transcripts that are expressed at a constant level during haustoriogenesis and can be used as reliable internal controls to correct for sample-to-sample variations in further experimental studies were first determined. Applying filtering criteria (see Section 2.10) resulted in a list of the 716 most stable transcripts from 688 gene models, which were considered as candidate housekeepers for in vivo investigation and normalization of gene expression. Among these transcripts, 5 from the top 20 that were assigned to different MapMan4 categories were selected for further investigation (Tables 1 and S2).
Marker genes that can be used to assign samples to a specific stage in later analyses of parasite-host interactions were then identified. Applying filtering criteria (see Section 2) resulted in a list of 585 transcripts from 556 gene models that were considered as candidate markers for the different haustorium development stages, including 68 transcripts (67 genes) for niS, 154 (149) for SWE, 7 (7) for ATT and 356 (333) for PEN. For each stage, three marker candidates from different gene models were selected for validation (Tables 2 and S3).
Due to the limited number of marker candidates for ATT, whose low accumulation level (mean RPKM <5 in most cases) in this stage made them difficult to detect by quantitative RT-PCR techniques, the filtering criteria had to be relaxed by allowing that not all possible differential comparisons between this stage and the others were significant.
This yielded Cc010463.t1 as a candidate.

| Stage-specific markers provide for the classification of natural infection sites
The expression of reference and marker candidates was investigated by quantitative real-time PCR (RT-qPCR) (Appendixes S2 and S3). To this end, the host-free induction experiment was reproduced so that a new set of RNA samples from three biological replicates for each haustorial stage was obtained. Moreover, infection sites of C. campestris on the tomato (S. lycopersicum) cultivar M82 were collected, their morphology documented by microscopy after sectioning and their RNAs extracted.
The stability of the five selected reference candidates was first evaluated using three different statistical approaches: geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004) and BestKeeper (Pfaffl et al., 2004). Rankings provided by all three methods were integrated by calculating the geometric mean for each accession (Table 3). In both the host-free and host-induced systems, only two genes were found to be sufficient for proper normalization, as indicated by the low average expression stability M (<1) and pairwise variation V2/3 (<0.2) values provided by the geNorm algorithm (Table 3, Figure S3). However, ranking integration provided contrasting results, with Cc028808.t1 and Cc006757.t1 being the most suitable pair in the host-free system, and Cc028378.t1 and Note: Expression values are average log 2 (RPKM) ± SE across all samples (four types, three biological replicates each). "CV" and "MAD" are respectively coefficient of variation (commonly used) and median absolute deviation (more robust to outliers), with a low value indicating a more stable expression. "Bincode" and "Description" refer to MapMan4 v.3.0 functional annotation. After normalization of transcript abundance, the selected markers mostly showed a consistent expression pattern in the parallel hostfree experiment compared with the RPKM values observed in the sequenced samples ( Figures 3A,B and S4). Altogether, the marker triples allowed for a proper clustering of the parallel host-free biological replicates ( Figure 3C). Finally, for a proof of principle, 11 infection sites were randomly selected on S. lycopersicum M82, which is tolerant to C. campestris, and cross-sectioned ( Figure 4A  Note: Expression values are average log 2 (RPKM) ± SE across all samples (four types, three biological replicates each). "Max" and "Min" are respectively maximum and minimum expression. "Stage" refers to the development stage toward which a marker is directed. "Zscore" refers to the average z-score in the corresponding stage. "GS" and "p.GS" are gene significance and corresponding P-value for that stage. "Bincode" and "Description" refer to MapMan4 v.3.0 functional annotation. Abbreviations: ATT, attaching stage; niS, non-infective stem; PEN, penetrating-like stage; SWE, swelling stage. Note: "M" refers to geNorm's average expression stability value; SV refers to NormFinder's stability value; "r" refers to BestKeeper's Pearson correlation coefficient; "GM" refers to geometric mean of the rankings from the three methods.

T A B L E 3 Stability ranking of the selected reference candidates in host-free and host-induced systems
interface region ( Figure 4C), highlighting the suitability of the markers also in a host-induced system.

| DISCUSSION
One major breakthrough in functional genomic studies in dodders was the almost simultaneous release of two reference genomes (Sun et al., 2018;Vogel et al., 2018). Here, we built on this opportunity to generate a comprehensive and continuous host-free transcriptome in C. campestris using RNA sequencing and explore the dynamic changes in gene expression. Unlike previous transcriptomic studies (Kaga et al., 2020;Olsen et al., 2016;Ranjan et al., 2014), we went beyond the simple distinction between young and mature haustoria, and defined three successive and macroscopically distinctive stages of development that can be induced without a host present: swelling, attaching and penetrating. The rationale behind the use of a host-free induction system as used in this study, is that it greatly facilitated the sampling of morphologically and transcriptionally uniform samples with high reproducibility. The representation of mature haustoria that are in a feeding and interacting modus with the host is perhaps more limited, but with the refinement of the host-free support system for . Infection sites with several coils (H1 and H2) were divided into two or three subsamples (À1, À2, À3). Swelling (SWE) stages appeared as smooth to slightly bumped parasite stems (H1-1, H2-1, H4, H5, H7 and H8). Attaching (ATT) stages appeared as clearly bumped parasite stems with no visibly emerged haustoria (H1-2, H2-2 and H6). Penetrating (PEN) stages appeared with visible endophytic haustoria (H2-3 and H3). (C) RT-qPCR transcript accumulation transformed into z-scores. Expression values were normalized against Cc028378. t1 and Cc002986.t1. Hierarchical clustering was based on Euclidean distance. A z-score value is positive (negative) if the transcript accumulation in a sample is larger (smaller) than the mean row accumulation. niS, non-infective stem expansion and division in the cortical layers of the stem and resulting in the swelling at the initiation sites (Kokla & Melnyk, 2018). Here, protein-and hormone-related categories were highlighted in SWE.
That such pathways were activated indicates that the parasite perceived the external stimuli to which it was exposed and responded to them by inducing haustoriogenesis. Protein phosphorylation plays a major role in signal transduction. In plants, protein phosphorylation has been implicated in responses to many signals, including light, thanks to a large array of different protein kinases (Stone & Walker, 1995). The current data suggest a large contribution of leucine-rich repeat (LRR) receptor-like kinases, which are known to be involved in growth and development (Liu et al., 2017). Similarly, phytohormones are key components in the development of infective structures in parasitic plants (Kokla & Melnyk, 2018 (Tada et al., 1996). The validity of ATT as a discernable stage, although relatively transient in our host-free system and with remarkably few indicative DEGs, is sustained by the sequencing data. Functional categories related to pectin modification and degradation, representing genes that could promote host surface attachment, were enriched at this stage, consistent with the fact that the parasite is getting prepared for invasion (Yoshida et al., 2016). Because a host is lacking, there is no further incentive for the haustoria to differentiate beyond the penetrating stage, although feeding-hyphae-like structures can sometimes be observed when the haustorium is prevented from drying out (data not shown). Recent observations that a fluorescent dye was taken up with amazing speed by the haustoria and the surrounding adhesive ring (Lachner et al., 2020), support the notion that some form of uptake competence may already exist at this point. More investigations are needed in the future to explore this possibility. Haustoria that are induced in the absence of a host do not contain xylem bridges that are required to form vascular connections (Kaga et al., 2020;Yoshida et al., 2016). It should nevertheless be noted that among the hormone-related functional categories that were enriched, TDIF peptide receptors (TDR) can be found, which are involved in cellular differentiation into tracheary elements (Morita et al., 2016). This could indicate that the parasite has already developed the competence to differentiate these structures that are crucial for its survival at this point and is awaiting the necessary cues from a host to proceed.
The analysis of expression patterns in the different stages revealed both constitutively expressed genes and genes with dramatic shifts in their expression. Further analyses were carried out to identify genes that can facilitate in vivo investigations of haustoriogenesis.
Among common methods for gene expression analysis, RT-qPCR has emerged as a fast, reliable and easy to use technique to measure transcript abundance in response to developmental variations and experimental conditions (Pabinger et al., 2014). Here, a new set of reference genes for RT-qPCR was identified along with stage-specific markers and validated in both a host-free and host-induced system. Some of the markers had no known function based on MapMan4 predictions, while the function of others reflects their involvement in haustoriogenesis and/or host invasion. For instance, Cc008389.t1 (PEN) is an alpha-class cell wall expansin which (as mentioned earlier) is expected to play a role in the regulation of cell wall expansion, while Cc004177.t1 (PEN) is a berberine bridge enzyme-like protein that possibly inactivates oligogalacturonides released upon pectin alteration and prevents them from triggering either the parasite or the host immunity (Benedetti et al., 2018). We considered three markers as a strict minimum for the molecular characterization of a sample and its assignment to one of the stages. This number can be seen as a balance between cost effectiveness and accuracy. It should be noted that the markers, although highly correlated with one of the development stages, often show weaker expression in others. It is therefore important that the entire set of markers is tested on a sample for its proper characterization. Also, infection sites can vary in length and a gradient in development stages among the haustoria can be observed. As a result, longer sites consisting of several coils must be split so that the less extensive the sites or the less haustoria per site or sample, the better for the marker pattern interpretation.

ACKNOWLEDGMENTS
We thank Leidulf Lund for plant care and maintenance, and Prof.

DATA AVAILABILITY STATEMENT
The generated raw reads from this article can be found in the NCBI Sequence Read Archive (SRA) under the accession number PRJNA666991.