Phytophthora cinnamomi isolates Pc2113 (A2 clade I) and Pc2109 (A2 clade II) were obtained from the Manosalva's laboratory pathogen collection, the CBS 144.22 isolate from the World Phytophthora Genetic Resource Collection (WPC) (University of California, Riverside), and the P. infestans isolate 1306-C (Pi1306-C) from Dr. Howard Judelson, UCR (Pan et al., 2018). Phytophthora isolates in this study were maintained as water agar plugs (Boesewinkel, 1976) and stored at 22°C.

Raphanus sativus cv. Saxa (radish), Solanum lycopersicum cv. M82 (tomato), and Nicotiana benthamiana seeds were germinated and grown under a 16-h light/8-h dark photoperiod at 22°C with 70% relative humidity. Two-week-old tomato seedlings were transferred to a greenhouse with an average air temperature between 25 and 28°C and an average humidity of 56%. Clonal Dusa^® avocado rootstock liners (Westfalia Technological Services, Tzaneen, South Africa) were obtained from Brokaw Nursery LLC (Santa Paula, Ventura, CA). The 6-month-old liners were removed from their bags and transplanted into pots after the nurse seed was removed. Plants were grown in a greenhouse under the conditions described above. Arabidopsis thaliana ecotype Landsberg erecta (Ler-0) were germinated and grown in a Conviron Adaptis A1000 growth chamber (Controlled Environments Limited, Winnipeg, Manitoba, Canada) at 23°C, 60% relative humidity, with a 16-h light/8-h dark photoperiod.

All Phytophthora cinnamomi isolates used in this study were grown in a liquid Plich media (Van der Lee et al., 1997) for 5–7 days at 22°C in the dark. The mycelia were harvested, washed with sterile Milli-Q water (MilliporeSigma, Burlington, MA, USA), flash-frozen in liquid nitrogen, and ground to a fine powder using a mortar and pestle. Approximately 100 mg of mycelia were transferred to microcentrifuge tubes and stored at −80°C until needed. Pc2113 and Pc2109 high-molecular-weight genomic DNA (gDNA) was extracted using the OmniPrep™ for Plants kit (G-Biosciences, St. Louis, MO, USA) according to the manufacturer's protocol, except for an additional RNAse treatment which involved adding 1.5 μl of LongLife™ RNAse (G-Biosciences, St. Louis, MO, USA) to the DNA solution and incubation at 37°C for 30 min. CBS 144.22 gDNA was extracted using the DNeasy Plant Mini kit (Qiagen, Germantown, MD) according to the manufacturer's instructions.

Pc2109 and Pc2113 gDNA were subjected to two independent Pacific Biosciences (PacBio) RS II and Sequel (PacBio, Menlo Park, CA) library preparations and sequencing runs at both Icahn Medical Institute at Mount Sinai (New York, NY) and Novogene Corporation (Sacramento, CA), respectively. Paired-end gDNA libraries with an insert size of 350 bp were generated for each isolate and sequenced using the Illumina NovaSeq 6000 PE150 platform (Illumina, San Diego, CA) by Novogene Corporation.

Flow cytometry (FCM) was conducted before genome assembly to determine the nuclear content. Phytophthora isolates were grown in 10% clarified V8 broth (cV8) (Belisle et al., 2019b) under the conditions indicated above for 3 days. Mycelium harvest, nuclei extraction, and staining were performed using the Cystain^® PI Absolute P kit (Sysmex Co., Kobe, Hyogo, Japan) as described by Bertier et al. (2013). Raphanus sativus cv. Saxa and/or Solanum lycopersicum cv. M82 leaf tissue served as reference standards. Approximately 10,000 nuclei per sample were measured in triplicate using MoFlow^® Astrios EQ Cell Sorter (Beckman Coulter, Indianapolis, IN) equipped with a 200 mW laser (488 nm). The data was analyzed using FlowJo™ software v10.7 (Tree Star Inc., Ashland, OR), and the genome size estimations were calculated as described by Doležel et al. (2007). The diploid P. infestans isolate Pi1306-C nuclear content was estimated and served as an FCM control. Genome size estimation was also performed using k-mer counting. K-mer histograms were generated with trimmed Illumina reads using Jellyfish (v2.3.0) (k-mer range 17–99 increments of 7) (Marçais and Kingsford, 2011), and the respective k-mer histograms served as the input for GenomeScope (accessed 20 November 2018) (Vurture et al., 2017) profiling.

PacBio subreads were combined and used for each de novo genome assembly using Canu v1.7.1 (Koren et al., 2017). The resulting assemblies were corrected with PacBio reads iteratively three times using Racon v1.3.2 (Vaser et al., 2017). The Illumina reads were filtered using Trimmomatic v0.36 (Bolger et al., 2014) and used to polish the Racon-corrected assemblies iteratively three times using Pilon v1.22 (Walker et al., 2014). Purge Haplotigs (Roach et al., 2018) were used to obtain a haploid genome assembly size that was consistent with our FCM estimations. Benchmark Universal Single-Copy Orthologs (BUSCO) v5.7.1 (Manni et al., 2021) using the Stramenopile database (orthoDB v.10) (Kriventseva et al., 2019) was used to assess assembly completeness. Assembly statistics were generated using QUAST v4.6.3 (Gurevich et al., 2013), and their quality and heterozygosity were assessed with KAT v2.3.4 (Mapleson et al., 2017). Pc2113 and Pc2109 whole genome alignments were performed using Nucmer (MUMmer v4.0), and differences were assessed using DNAdiff (Marçais et al., 2018). Dot plots were generated using Minimap2 v2.10 (Li, 2018) and visualized using D-GENIES v1.3.1 (Cabanettes and Klopp, 2018). Cross-read mapping was performed with Pc2113 and Pc2109 raw Illumina reads to the GKB4 P. cinnamomi reference genome (NCBI accession number PRJNA675400) using BWA-mem v0.7.17 (Li, 2013), and mapping percentage was acquired using Samtools v1.18 flagstats (Li et al., 2009).

First, repeat annotation and transposable element identification were performed using RepeatModeler 2 v2.0.1 (http://www.repeatmasker.org) for both assemblies. Assemblies for Pc2113 and Pc2109 were merged and used to build the database for each isolate. The respective database was used to annotate repeats, and a soft mask was performed for each genome using RepeatMasker v4.0.9 (http://www.repeatmasker.org) before conducting the gene annotation process using RNA-seq data as supporting gene evidence. The RNA sequences were trimmed and filtered using fastq-mcf v1.05 (Aronesty, 2011) (-q 30 -l 50), and alignments were performed with STAR v2.7.5a (Dobin et al., 2013) and assembled using StringTie v2.1.5 (Pertea et al., 2015). GFF utilities (Pertea and Pertea, 2020) were used to convert the GTF from StringTie to GFF to extract mRNA sequences. Peptide prediction was made with TransDecoder v5.5.0 (https://github.com/TransDecoder/), DIAMOND v0.9.24 (Buchfink et al., 2014), and HMMER3 (Eddy, 2009) (hmmscan v3.1b2). BRAKER v2.1.5 (Hoff et al., 2019) pipeline, GeneMark v4.61 (Stanke et al., 2006), and Augustus v3.3.3 (Stanke and Morgenstern, 2005) were employed to predict genes. Gene annotation was also conducted using MAKER pipeline v3.01.02, for which the ab initio gene prediction program SNAP v2006-07-28 (Korf, 2004) was trained with the mRNA sequences and peptides described above. The final round of MAKER was completed using the proteins from TransDecoder, including P. infestans v2 (NCBI accession number PRJNA17665), P. ramorum v1 (JGI Project ID 16273), P. sojae (JGI Project ID 16274), P. parasitica v2 (NCBI accession number PRJNA73155), P. capsici v11 (JGI Project ID 16272), P. cinnamomi v1 (JGI Project ID 1003774), and P. nicotianae v.1 (NCBI accession number PRJNA294216). BUSCO v4.1.4 (Simão et al., 2015) was used to assess the completeness of the annotations (BRAKER & MAKER) as described above, and MAKER and BRAKER annotation comparison was performed using Mikado v2.0rc2 (Venturini et al., 2018). BRAKER was chosen as the “base” annotations, and the MAKER annotations served as the “donor”. Unique genes from the maker were joined with the BRAKER annotations and the GFF file was sorted (https://github.com/billzt/gff3sort, v2017). DIAMOND was conducted for functional annotations using the SwissProt and TREMBL datasets (accessed on 3 September 2019). Descriptions were added to the annotations using AHRD v3.3.3 (https://github.com/groupschoof/AHRD).

Phytophthora cinnamomi and P. infestans ploidies were assessed using nQuire (Weiß et al., 2018) and the vcfR package (Knaus and Grünwald, 2017) in R v4.1.2 (R Core Team, 2010). The Illumina reads were trimmed, filtered using Fastq-mcf v1.1.2-537 (Aronesty, 2011), and the reads were mapped to their corresponding reference genomes as indicated above. The respective BAM files were processed using Samtools and assessed with Qualimap v2.2.1 (García-Alcalde et al., 2012). Single-nucleotide polymorphisms (SNPs) were called from the processed BAM files using FreeBayes v1.2.0 (Garrison and Marth, 2012). The respective variant call format file (VCF) was filtered with VCFtools v0.1.16-18 and imported into R for allele balance analysis following the methods described in Knaus and Grünwald (2018). The P. infestans isolate 1306-C was subjected to the same analysis as a diploid control using the T30-4 reference genome (NCBI accession number ASM14294v1).

Phytophthora infestans (NCBI accession number PRJNA17665), P. ramorum (NCBI accession number PRJNA12571), P. sojae (JGI Project ID 16274), P. fragariae (NCBI accession number PRJNA396163), P. parasitica (NCBI accession number PRJNA73155), and P. capsici (NCBI accession number PRJNA481983) proteomes were used to conduct orthology analysis using Orthofinder v2.5.2 (Emms and Kelly, 2019). Orthogroup intersections were visualized using UpsetR v1.4.0 (Conway et al., 2017) and extracted using ComplexHeatmap v2.6.2 (Gu et al., 2016) for downstream analysis. To identify proteins only present in P. cinnamomi, a tBLASTn v2.9.0+ analysis was performed (Altschul et al., 1997). tBLASTn results were filtered based on ≥80% identity and ≥80% of the query alignment length. The respective proteins were then subjected to Exonerate v2.4.0 (Slater and Birney, 2005) and filtered based on ≥80% identity and similarity. In addition, separate orthology analyses using GKB4 (Engelbrecht et al., 2021), Pc2113 and Pc2109 proteomes, predicted carbohydrate-active enzymes (CAZymes), and RXLRs were performed.

To determine the genome architecture for each P. cinnamomi isolate, two-dimensional gene binning based on their flanking intergenic regions' (FIRs) length was performed using methods adopted by Saunders et al. (2014) and visualized using the density-Mapr R script (https://github.com/Adamtaranto/density-Mapr). Gene-dense regions (GDRs) and gene-sparse regions (GSRs) of the genomes were delimited using methods described by Raffaele et al. (2010) and Rojas-Estevez et al. (2020) (Supplementary Method S1). The proteome for each P. cinnamomi isolate was scanned for signal peptide (SP) presence using SignalP (v5.0) (Armenteros et al., 2019). The resulting proteins containing SP sequence were subjected to TMHMM (v.2.0) analyses (Möller et al., 2001) to identify transmembrane domains (TMDs). TMD-containing proteins were removed from the secretome datasets. Secretomes were subjected to EffectorP v3.0 (Sperschneider and Dodds, 2022) to predict apoplastic, cytoplasmic, and dual-localized effectors. In addition, RXLR effectors were also predicted using regular expression (REGEX) searches following methods described by Win et al. (2007) and Boutemy et al. (2011) (Supplementary Method S2). Multiple sequence alignments of the predicted Pc2113 and Pc2109 RXLRs, 61 GKB4 RXLRs, and 41 RXLRs from other Phytophthora spp. and Hyaloperonospora arabidopsis described in Joubert et al. (2021) were performed with MUSCLE v3.8.425 (Edgar, 2004) using default settings. A maximum likelihood phylogenetic tree was generated using IQ-Tree2 v2.1.3 (Minh et al., 2020) with the VT+R7 model determined by ModelFinder (Kalyaanamoorthy et al., 2017). The phylogenetic tree was visualized using the ggtree v3.2.1 (Yu et al., 2017) and ComplexHeatmap R packages. CAZymes were identified using the dbCAN2 web server (Yin et al., 2012; Zhang et al., 2018) using the CAZyme database (Cantarel et al., 2009), as described by Engelbrecht et al. (2021). Proteins were considered CAZyme if they were positive for two or more programs (HMMER, Hotpep & Diamond). Proteins were putatively classified as cell wall-degrading enzymes (CWDEs) if they contained polysaccharide lyase (PL), carbohydrate esterase (CE), glycoside hydrolase (GH), auxiliary activity (AA), or carbohydrate-binding module (CBM) modules, which are associated with plant cell wall degradation (Blackman et al., 2014). Necrosis-inducing proteins (NIPs), carbonic anhydrases (CAs), crinkling and necrosis (CRNs), elicitins, and elicitin-like proteins were identified from the functional annotations generated in this study.

Leaves of 5-week-old N. benthamiana, 9-month-old Dusa^® avocado, and three-and-a-half-week-old A. thaliana Ler-0 were used for Pc2113 detached leaf inoculations as described in Belisle et al. (2019b). At 16- and 24-h post-inoculation (hpi), three infected leaf discs were harvested from three individual plants as one biological replicate per time point. Three biological replicates were used for this RNA-seq experiment, apart from A. thaliana samples at 90 hpi, which consisted only of two biological replicates. Leaf tissue was frozen in liquid nitrogen, ground to a fine powder using a mortar and a pestle, and stored at −80°C. Pc2113 mycelia bearing sporangia were harvested, vacuum-filtered, ground to a fine powder in liquid nitrogen with a mortar and pestle, and stored at −80°C. Pc2113 zoospore suspension (2 × 10⁵ zoospore/ml) was centrifuged at 362 × g for 5 min and the resulting pellet was resuspended in 1 ml of sterile water, ground in liquid nitrogen using a mortar and a pestle, and stored at −80°C. Mycelia-bearing sporangia and zoospore samples were harvested in three biological replicates. Total plant and pathogen RNA were extracted using the RNeasy Plant Mini Kit (Qiagen, Germantown, MD, USA) according to the manufacturer's instructions, followed by DNAse treatment with the Invitrogen TURBO DNA-free™ Kit (Thermo Fisher Scientific, Waltham, MA). RNA cleanup and purification were performed using the Zymo RNA Clean and Concentrator Kit (Zymo Research Corp., Irvine, CA). RNA pellets were resuspended in nuclease-free water, and the total RNA concentration and absorbance 260/280 ratio, indicative of RNA purity, were measured using a NanoDrop (Thermo Fisher Scientific, Waltham, MA). To confirm successful infection and pathogen accumulation before submitting samples for RNA sequencing, P. cinnamomi biomass at the respective time points post-inoculation was calculated using TaqMan real-time PCR as described by Belisle et al. (2019b). Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) using the Plant RNA Nano Assay was used to assess the RNA quality and purity before sending samples to Novogene for library preparation and 250–300 bp paired-end sequencing using the Illumina HiSeq2 platform. RNA sequencing reads were mapped to the Pc2113 assembly with HISAT2 v2.1.0 (Kim et al., 2019) using default settings. Samtools v1.12 was used for alignment mapping file processing, and reads were counted with HTseq-count v0.11.1 (Anders et al., 2015). Differential gene expression analysis was conducted using the R package DEseq2 v1.30.1 (Love et al., 2014). Genes with < 10 counts across the rows of the matrix were removed, and differential expression analysis was performed using the DEseq function. Comparisons within DEseq2 were performed by comparing the in planta timepoints to the Pc2113 mycelium (bearing sporangia) using the contrast function with the following settings: alpha = 0.05 and lfcThreshold = 1. Genes were characterized as differentially expressed (DE) if the false discovery rate adjusted p values (Benjamini and Hochberg, 1995) were of < 0.05 and had a |log-fold change (LFC)| of ≥ 2. Genes with an LFC of ≥ 2 were considered upregulated, and genes with an LFC of ≤ -2 were considered downregulated. The same method was applied to compare Pc2113 sporulating mycelium and Pc2113 zoospore samples to identify DE genes in the zoospore sample. Gene expression in terms of transcript per million (TPM) was generated using the same bam files described previously using StringTie v2.2.1 (Pertea et al., 2015). The transcript TPM values were extracted from the output GTF files and averaged across each of the sample replicates at each timepoint, respectively. Genes were considered expressed if the TPM > 5. The respective TPM values were log₂ transformed for visualization using ComplexHeatmap v.2.6.2 in R.

qRT-PCR was used to validate expression patterns identified by the RNA-seq data with the primers in Supplementary Table S1. qRT-PCR analysis of three highly upregulated P. cinnamomi effector genes was performed to validate RNA-seq results. Detached leaf inoculations with Pc2113 in N. benthamiana were conducted as described above. Leaf discs were harvested for three replicate samples, with three plants per replicate at each time point. cDNA synthesis was performed using the Invitrogen Superscript III Kit (Thermo Fisher Scientific, Waltham, MA) following the manufacturer's protocol. Each 10-μl reaction included 4 μl of a cDNA (1:20 dilution), 5 μl of iQ SYBR Green Supermix (Bio-Rad Laboratories, Inc., Hercules, CA), and 0.5 μM of each forward and reverse qRT-PCR primers. Reactions were amplified using a CFX384 Touch Bio-Rad Real-Time PCR system (Bio-Rad Laboratories) under the following conditions: 95°C for 3 min (initial denaturation), followed by 40 cycles of denaturation at 95°C for 10 s, annealing/extension at 60°C for 30 s, and a melt curve of 5 s at 65–95°C. Gene expression levels were calculated using the 2^−ΔΔCt (Ct = cycle threshold) comparative method (Livak and Schmittgen, 2001). The results were reported as the mean ± standard error of three biological replications and three technical replications. Pearson's correlation test was used to evaluate significance. All primers in this study were designed using Primer3 v0.4.0 software. The three candidate gene primer pairs used (sequences not shown) were Pc2113T1C00002279g0000010.1 (RXLR052), Pc2113T1C00002373g0000080.1 (Elicitin2373), and Pc2113T1C00002266g0000100.1 (RXLR047) (Supplementary Table S1). Phytophthora cinnamomi 40s ribosomal gene Ws21 (86101) was used as an endogenous reference gene for normalization. Primer efficiency for each pair was calculated via the slope of a standard curve constructed by the amplification of a 10-fold dilution series of DNA [covering five dilution points, 1–10,000 picograms (pg)]. All primer sets were between 95% and 99% efficient. This experiment was conducted two times.

Raw Illumina reads post-filtering resulted in a total of ~156 and ~130 million reads, representing 23.1 and 19.1 GB of filtered data for Pc2113 and Pc2109, respectively (Supplementary Table S2). Genome size estimations by GenomeScope yielded a size range from 144.5 Mb (17-mer) to 178.6 Mb (99-mer) for Pc2113 and from 140.4 Mb (17-mer) to 183.6 Mb (86-mer) for isolate Pc2109 (data not shown). Recent studies reported the predominance of triploidy in P. cinnamomi populations collected from avocado orchards in South Africa (Engelbrecht et al., 2021, 2022). Considering that Pc2113 and Pc2109 genome size estimations were close to the South African triploid P. cinnamomi isolate GKB4 (Engelbrecht et al., 2021), we used a triploid model to calculate the haploid genome size of the P. cinnamomi isolates by flow cytometry. Consistent with Van Poucke et al. (2021), all the P. cinnamomi isolates from this study exhibited similar DNA content, with an average ranging from 0.41 to 0.43 pg (Supplementary Figure S1). The use of R. sativus (2C = 1.09 pg) (Arumuganathan and Earle, 1991) as the reference standard resulted in a slight overlap of the P. cinnamomi G2_P and R. sativus G1_R nuclear populations (Supplementary Figures S1A–C). To resolve the overlap, S. lycopersicum (2C = 2.01 pg) (Marie and Brown, 1993) was substituted as a reference standard (Supplementary Figure S1D). No significant differences (p = 0.7187, Kruskal–Wallis rank sum test) were observed in the estimated DNA content (pg) for the Pc2113 isolate using either reference standard (Supplementary Figure S1E). As expected, we obtained an estimated genome size of 0.54 pg for P. infestans isolate Pi1306-C (Catal et al., 2010; Wang et al., 2016) (Supplementary Figure S1E). These results indicate a larger genome size for California P. cinnamomi isolates compared to the other publicly available P. cinnamomi genome assemblies, and these estimations were used to guide the respective de novo genome assemblies.

Post-filtering of the PacBio reads resulted in approximately 1.8 million subreads for each isolate, representing ~20 Gb and ~1 Gb for Pc2113 and Pc2109, respectively (Supplementary Table S2). The initial raw Canu assemblies yielded 259.9 Mb and 232.9 Mb for Pc2113 and Pc2109, respectively. Three iterations of Racon followed by three iterations of Pilon increased the quality of the assemblies but resulted in larger genome sizes than those estimated by the K-mer and FCM analyses. Purge Haplotigs successfully reduced the number of heterozygous contigs observed within the assemblies to match the FCM estimations under triploidy. The resulting assemblies representing the haploid state yielded a length of 140.7 Mb and 141 Mb for the Pc2113 and Pc2109 genomes, respectively, in 735 and 897 contigs with N50 values ranging from 289 to 255 kb. BUSCO analyses indicated a completeness of 100% for each isolate, and 77–86% were single BUSCOs with duplications of 23% and 14% for Pc2113 and Pc2109, respectively. Genome annotation for Pc2113 and Pc2109 yielded 21,402 and 21,111 genes, respectively. In addition, we found 375 and 399 alternatively spliced variants and a total of 21,777 and 21,510 proteins for Pc2113 and Pc2109, respectively. The secretome of each isolate is comprised of 1,400 and 1,342 proteins for Pc2109 and Pc2113, respectively (Table 1).

The ploidy of P. cinnamomi isolates Pc2113, Pc2109, and CBS 144.22 was confirmed as triploid based on two different approaches. Phytophthora cinnamomi allele balance histograms generated from vcfR displayed a binomial distribution with two dominant peaks at 1/3 and 2/3 typical for triploid isolates (Supplementary Figures S2A–C). As expected, the diploid Pi1306-C control isolate showed an allele balance histogram displaying a dominant peak at 12, indicative of diploidy (Supplementary Figure S2D). nQuire ploidy models for P. cinnamomi isolates generated the lowest delta log-likelihood values for the triploid model, while the lowest value for P. infestans was for the diploid model (Supplementary Figure S2E).

High macrosynteny and sequence identity (>75%) were observed in the dot plot between Pc2109 and Pc2113, and each assembly with GKB4 (Figure 1A, Supplementary Figures S3A, B). Cross-read mapping of the California P. cinnamomi Illumina reads to the GKB4 assembly yielded a mapping percentage of >97%, supporting the observations in the respective dot plots. Pc2113 and Pc2109 genome comparison using DNAdiff also showed that 99.7% of their total contigs aligned with ~88% of sequence identity. No major differences were detected regarding breakpoints between both genomes; however, Pc2113 showed 539 more events than Pc2109, including translocations and insertions. Interestingly, DNAdiff revealed that the Pc2109 genome exhibited more tandem insertions (105) than Pc2113 (89), with an average size of 1.6 kb (Figure 1B). KAT k-mer density plots were generated to assess genome complexity. Both P. cinnamomi isolates exhibited similar distributions, with a strong heterozygous k-mer signal at a quarter coverage and two homozygous signals, one at half coverage (weak signal) and another at full coverage (strong signal) (Figures 2A, B). KAT spectra-cn plots were used to investigate the k-mer frequencies in common between the assemblies of Pc2113 and Pc2109 and their respective read sets. Consistent with the KAT density plots, each isolate displayed three distinct peaks (one heterozygous and two homozygous) (Figures 2C, D). The first homozygous peak is well represented in the respective assemblies, with a small portion of the read set (black) not captured in the assembly (red), while the second homozygous peak is almost completely represented in the assembly. In addition, both homozygous peaks from each isolate exhibited different levels of k-mer duplications (2 × and 3 × ) (Figures 2C, D). Together, these results indicated completeness and high genome complexity, with different levels of duplications supporting polyploidy and/or aneuploidy events.

Orthology analysis revealed that Pc2113 and Pc2109 shared 12,777 orthogroups, representing 88% of their respective proteomes. In Pc2113, there were 754 orthogroups containing 773 proteins that were absent in Pc2109. Interestingly, Pc2109, the most virulent isolate, exhibited more orthogroups (n = 836) comprised of 929 proteins not shared with Pc2113 (Figure 3A). All Phytophthora spp. compared shared a total of 7,457 orthogroups, which represent ~51% of the total proteins in Pc2113 (n = 11,030) and Pc2109 (n = 10,904) (Figure 3B). In addition, 2,540 of these orthogroups were single-copy “core” orthologs. All the proteins and their corresponding encoding genes shared among all genomes analyzed will be referred to as “OG7457” for the remainder of this article. A total of 322 orthogroups are shared among all the Phytophthora spp. analyzed, except for the California P. cinnamomi isolates. Interestingly, Pc2109 and Pc2113 contain 904 orthogroups that are not present in the other Phytophthora spp. analyzed, representing ~6.2% of their proteomes (1,349 and 1,308 for Pc2113 and Pc2109, respectively). Proteins that belong to these orthogroups will be referred to as “OG904” for the remainder of this article. A total of 148 and 136 orthogroups were also identified only in Pc2113 (386 proteins) and in Pc2109 (362 proteins), respectively. To determine the total number of proteins only present in each isolate, we combined these proteins with the Pc2113 and Pc2109 unassigned proteins. After removing positive hits using tBLASTn and Exonerate, we calculated a total of 1,019 and 1,165 proteins present only in Pc2109 (OG-Pc2109) and Pc2113 (OG-Pc2113), respectively.

The optimal L value for delimiting the genome into gene-dense and gene-sparse regions was determined to be 1.1 Kb (Supplementary Figures S4A, B). Delimitation of the genomes with an FIR L value of 1.1 Kb yielded 8,281 (38.7%) and 7,849 (37.2%) genes in GDRs for Pc2113 and Pc2109, respectively. Approximately 37% of the genes from Pc2113 (8,030) and Pc2109 (7,853) resided in the in-between regions, and 18% of Pc2113 (3,938) and Pc2109 (3,970) genes were in GSRs. The positions of < 7% of Pc2113 (1,153) and Pc2109 (1,439) genes could not be determined (ND) (Supplementary Figures S4C, D). The genome delimitation allowed us to determine the distribution of predicted genes associated with pathogenicity and virulence within these regions. Functional annotation revealed that Pc2109 encoded slightly more NIPs and elicitins than Pc2113, but fewer CAs and CRNs. Overall, regions where these effectors resided were consistent among the two isolates; however, Pc2109 displayed slightly higher effector counts within GSRs than Pc2113, except for the elicitins. A higher proportion of predicted NIPs (~50%), CAs (~50%), elicitins (~46%), and CRN (~34%) shared orthology with other Phytophthora spp. (OG7457) (Table 2). A total of 632 and 583 effectors were predicted in Pc2109 and Pc2113 genomes, respectively, using EffectorP. The predicted subcellular localization of these effectors was similar for both isolates; however, Pc2109 exhibited more apoplastic elicitins. After excluding annotated/predicted elicitins, NLPs, CWDEs, CAs, RXLRs, and CRNs from the EffectorP predictions, a total of 332 and 333 effectors were identified for Pc2109 and Pc2113, respectively (Supplementary Table S3, Table 2).

CAZyme prediction yielded 378 and 372 putative CAZymes for Pc2113 and Pc2109, respectively (Supplementary Table S4). Isolate Pc2113 encoded more CWDEs (299) compared to Pc2109 (276). For both isolates, the GH family made up the majority of the predicted CWDEs, followed by the PL family. In both isolates, ~60% of CWDEs were common to other Phytophthora spp. (OG7457), and ~30% resided in GSRs. We identified several CWDEs that were present in the California P. cinnamomi isolates (OG904), with more than 50% residing in the GSRs. Isolate Pc2109 contained a higher proportion of private CWDEs (OG-Pc2109) in GSRs than Pc2113 (OG-Pc2113).

Isolate Pc2109 harbors 222 predicted RXLRs, 49 more than Pc2113 (Supplementary Table S5). Approximately 22% of the predicted RXLR effectors in both isolates contain a putative WY-domain based on the HMMER search (E-value of < 0.05). For Pc2109, more effectors resided in GSRs (102) and in-between regions (84) when compared with Pc2113 encoding 74 and 64 in those regions, respectively. On the contrary, more Pc2113 effectors (23) resided in GDRs compared with Pc2109 (14). Less than 10% of the effector distributions for each isolate were ND (Figures 4A, C). As a control, the majority of the selected Pc2113 and Pc2109 housekeeping genes (the 40s and 60s ribosomal proteins) resided in GDRs, with ~1% in GSRs (Figures 4B, D). Approximately 22% and 6% of the RXLRs predicted from both isolates were also present in the OG7457 and the OG904 clusters, respectively (Supplementary Table S5). Orthology analysis of Pc2113, Pc2109, and GKB4 proteomes revealed several protein clusters specific to each isolate and 12,061 orthogroups shared among all isolates, representing < 10% and 76% of their proteomes, respectively (Supplementary Figure S3C). The clustering of the predicted CAZymes and RXLRs from the two California isolates and the GKB4 also revealed common and isolate-specific orthogroups (Supplementary Figures S3D, E). We identified a cluster of 184 common orthogroups containing 362, 293, and 279 predicted CAZymes in GKB4, Pc2113, and Pc2109, respectively, representing ~75% of the total repertoire in each isolate. A total of 12 GKB4 orthogroups, corresponding to ~15% of the CAZyme repertoire, were not present in the California isolates. Only two Pc2113 orthogroups, corresponding to ~1% of the CAZyme repertoire, were not present in GKB4 or Pc2109 (Supplementary Figure S3D). Similarly, we found 45 orthogroups containing 89, 60, and 83 common RXLRs in GKB4, Pc2113, and Pc2109, respectively, representing ~34–49% of the total repertoire in each isolate. A total of nine, seven, and two orthogroups were only found in GKB4, Pc2109, and Pc2113 isolates, respectively. GKB4 and Pc2109 exhibited more private RXLRs (13%) when compared with Pc2113 (8%) (Supplementary Figure S3E).

Prior to RNA-seq, we performed detached leaf inoculations with Pc2113 infecting A. thaliana, N. benthamiana, and avocado and used trypan blue staining at 16, 24, and 90 hpi (A. thaliana only) to track Pc2113 infection (data not shown). The Pc2113 infection in A. thaliana at 90 hpi was consistent with the Pc2113 infection at 24 hpi in N. benthamiana and avocado, which led us to include the 90 hpi time point for A. thaliana. Following, we conducted an RNA-seq experiment to evaluate the expression patterns of Pc2113 when infecting A. thaliana, N. benthamiana, and avocado at 16, 24, and 90 hpi. Pathogen reads mapped ~91% and < 8% in Pc2113 growing in vitro and in planta, respectively. Interestingly, Pc2113 inoculation produced different lesions in each host evaluated, with more chlorotic and less necrotic lesions observed in the most susceptible A. thaliana ecotype (Robinson and Cahill, 2003), followed by N. benthamiana, which displayed water-soaking-like necrotic lesions, and avocado, which exhibited darker necrotic lesions. Lesion progression in each host was correlated with the increase of Pc2113 biomass in the infected tissue. In addition, the mapping percentage in each host increased over time and was correlated with the lesion phenotypes in each host (Supplementary Figure S5, Supplementary Table S6).

A total of 872 DEGs were identified across all timepoints in A. thaliana (546 up- and 326 down-regulated), 4,224 DEGs in N. benthamiana (1,162 up- and 3,048 downregulated), and 5,693 DEGs in avocado (2,048 up- and 3,645 downregulated). In avocado, 58% of the Pc2113 DEGs shared orthology with other Phytophthora spp. (OG7457), compared with 66% of the Pc2113 DEGs in A. thaliana and N. benthamiana. Less than 5% of the Pc2113 DEGs in each host belonged to the OG904 or OG-Pc2113 clusters. The distribution of Pc2113 DEGs that resided in GSRs was similar in all hosts analyzed (~17%) (Supplementary Table S7). A total of 249, 1,820, and 5,113 Pc2113 DEGs were identified in A. thaliana, N. benthamiana, and avocado, respectively, at 16 hpi, increasing over time, reflecting the disease development and lesion appearance when Pc2113 infected each host (Supplementary Figure S5). We observed a slight decrease in the number of Pc2113-upregulated genes in avocado and A. thaliana at the latest time point analyzed (Supplementary Table S8). RNA-seq analyses revealed Pc2113 DEGs in all hosts, two hosts, and specific hosts at 16 and 24 hpi (Figure 5). We found nine and eight Pc2113 genes significantly up- and down-regulated at 16 hpi (Figures 5A, C) and increased to 130 up- and 21 down-regulated at 24 hpi in all hosts, respectively (Figures 5B, D). There were more common subsets of upregulated genes between the two hosts, while only N. benthamiana and avocado shared overlapping sets of downregulated genes. Several subsets of upregulated genes decreased over time in A. thaliana and avocado but increased in N. benthamiana (Figures 5A, B). Similarly, the amount of downregulated genes decreased over time in N. benthamiana and avocado; however, no downregulated genes were found only in A. thaliana (Figures 5C, D).

Necrosis-inducing proteins, elicitins, and CAs have been shown to play important pathogenicity roles in several Phytophthora–plant interactions (Raffaele et al., 2010; Reitmann et al., 2017; Hardham and Blackman, 2018). The differential expression of these effectors varied among the hosts analyzed, with the least exhibited in A. thaliana, followed by N. benthamiana and avocado. Approximately 12, 8, and 1 of the predicted NIPs were upregulated when infecting avocado, N. benthamiana, and A. thaliana, respectively. The NIP upregulated in A. thaliana (Pc2113T1C00002230g0000120.1) was also upregulated in the other hosts. Interestingly, four of the NIPs were specifically upregulated in avocados. Only two NIPs were downregulated in avocado and no other hosts. Similarly, 17, 9, and 6 of the elicitins were upregulated and 7, 6, and 2 were downregulated in avocado, N. benthamiana, and A. thaliana, respectively. Four elicitins were found to be exclusively upregulated in planta (Pc2113T1C00002373g0000080.1, Pc2113T1C00002392g0000050.1, Pc2113T1C00002418g0002680.1, and Pc2113T1C00002607g0000450.1). Several elicitins exhibited host-specific DE in avocado, including the upregulated Pc2113T1C00002583g0000030.1 (OG904) and the downregulated Pc2113T1C00002220g0000590.1. The Pc2113T1C00002220g0000580.1 elicitin was only downregulated in N. benthamiana. Different subsets of CAs were DE, with 7, 4, and 2 upregulated and 5, 5, and 1 downregulated in avocado, N. benthamiana, and A. thaliana, respectively. The CA, Pc2113T1C00003123g0000220.1, was upregulated exclusively in avocados. Furthermore, none of the DE NIPs, elicitins, or CAs were A. thaliana-specific. The number of up- and down-regulated genes reported as host-specific is limited to those DE in planta (not DE in the zoospore stage) (Supplementary Table S7).

A total of 141, 80, and 23 of the predicted CWDEs (299) were upregulated in avocado, N. benthamiana, and A. thaliana, respectively. Additionally, 57 CWDEs were upregulated in zoospores. The majority of the upregulated CWDEs belonged to the GH CAZy family, followed by the PL family, except for A. thaliana, with only one PL upregulated. Only one member of the CBM CAZy family was upregulated in avocado and N. benthamiana (Supplementary Table S7, Figure 6A). A total of 14 CWDEs were upregulated in all hosts (13 GH and one CE). Several CWDEs were upregulated between two hosts, except for A. thaliana and N. benthamiana. A total of 45 CWDE were commonly upregulated in avocado and N. benthamiana (33 GH, 5 CE, 6 PL, and 1 CBM) and one GH in avocado and A. thaliana. A total of 51 CWDEs were exclusively upregulated in avocado, with 42 upregulated at both timepoints and 17 located in GSRs. Several CWDEs belonged to the OG7457 clusters (23) and two PLs, Pc2113T1C00000567g0000150.1 and Pc2113T1C00001098g0000110.1, to the OG904 and OG-Pc2113 clusters, respectively. Five CWDEs were exclusively upregulated at 24 hpi when infecting N. benthamiana, with all but one CE (Pc2113T1C00002636g0002400.1) located in GSRs. Only two, Pc2113T1C00002589g0001670.1 (PL) and Pc2113T1C00002397g0000460.1 (GH), shared orthology with Phytophthora spp. (OG7457). We identified four CWDEs exclusively upregulated in A. thaliana, including three GHs and one PL. Only two A. thaliana-specific CWDEs (GHs) resided in GSRs and shared orthology with Phytophthora spp. (OG7457) (Supplementary Table S7, Figure 6B).

In contrast to the upregulated CWDEs, there were fewer downregulated CWDEs in avocado (22), N. benthamiana (21), and A. thaliana (4), with the majority in the GH CAZy family. No downregulated CWDEs in the CBM family were detected, and members of the CE family were only found downregulated in zoospores (Supplementary Table S7, Figure 6C). A GH CWDE (Pc2113T1C00002345g0000380.1) was downregulated in all hosts. Only avocado and N. benthamiana shared seven downregulated CWDEs. Of these seven, two of the GHs shared orthology with Phytophthora spp. (OG7457), one PL (Pc2113T1C00001147g0000260.1) belonged to the OG904 cluster, and another PL (Pc2113T1C00002233g0000020.1) to the OG-Pc2113 cluster. Five CWDEs were exclusively downregulated in N. benthamiana. Of these five, only two of the GHs shared orthology with Phytophthora spp. (OG7457). Interestingly, the GH (Pc2113T1C00002958g0000060.1) and PL (Pc2113T1C00002426g0000100.1) enzymes were found to be exclusively upregulated in avocado and A. thaliana, respectively, but downregulated in N. benthamiana. Six CWDEs were exclusively downregulated in avocado, and four belonged to the OG7457 cluster. We did not detect any downregulated CWDEs exclusive to A. thaliana (Supplementary Table S7, Figure 6D).

Differential RXLR expression was detected when infecting different hosts. Approximately 27% (47) and 6.9% (12) of the predicted Pc2113 RXLRs (173) were up- and down-regulated across all samples and timepoints, respectively. In A. thaliana, the number of upregulated RXLRs decreased over the time course of the infection. Three out of the seven upregulated RXLRs in this host belonged to the OG7457 (Pc2113RXLR061, Pc2113RXLR098, and Pc2113RXLR092) and one (Pc2113RXLR047) to the OG904 clusters. No downregulated RXLRs were detected in this interaction. When infecting N. benthamiana, 10.4% (18) and 3.4% (6) RXLRs were up- and down-regulated during the infection, respectively. Of these upregulated RXLRs, seven resided in GSRs, seven shared orthology with Phytophthora spp. (OG7457), and two (Pc2113RXLR047 and Pc2113RXLR052) belonged to the OG904 clusters. There were two downregulated RXLRs in GSRs; one belonged to the OG7457 clusters, and two (Pc2113RXLR149 and Pc2113RXLR150) to the OG904 clusters. A total of 37 and 12 RXLRs were up- and down-regulated in avocados, respectively. In total, 15 of these upregulated RXLRs resided in GSRs; 13 belonged to the OG7457 clusters, and two (Pc2113RXLR047 and Pc2113RXLR052) to the OG904 clusters. Four downregulated RXLRs in avocado resided in GSRs, four shared orthology with Phytophthora spp. (OG7457), and two (Pc2113RXLR149 and Pc2113RXLR150) belonged to the OG904 clusters. We found 39 DE RXLRs (16 up- and 23 down-regulated) in zoospores. A total of 4 up- and 10 down-regulated RXLRs shared orthology with Phytophthora spp. (OG7457). In addition, three up- and two down-regulated RXLRs belonged to the OG904 cluster, and one downregulated RXLR (Pc2113RXLR077) belonged to the OG-Pc2113 clusters. A total of 31 and 3 RXLRs were exclusively up- and down-regulated in planta. There were RXLR subsets that were DE in all hosts, and some were host-specific. In total, 3 RXLRs were upregulated in all hosts (Pc2113RXLR047, Pc2113RXLR096, and Pc2113RXLR092), 11 were common to N. benthamiana and avocado, and 1 RXLR was overlapping between A. thaliana and avocado. We found one RXLR upregulated only in N. benthamiana and 15 in avocado. Interestingly, none of the downregulated RXLRs were commonly expressed in all hosts or between two hosts, and only three RXLRs were significantly downregulated in avocados (Supplementary Table S7).

Phylogenetic analysis revealed the evolutionary relatedness of the Pc2113 and Pc2109 RXLRs, along with previously characterized RXLRs from various Phytophthora spp. (Figure 7, Supplementary Figure S6). Several P. cinnamomi RXLRs effectors from this study clustered with a P. infestans suppressor of necrosis 1 (SNE1) effector (NCBI accession number ABI74673.1) (Kelley et al., 2010). This cluster contained three Pc2113 RXLRs (Pc2113RXLR058, Pc2113RXLR057, and Pc2113RXLR062) annotated as SNE1, with two that resided in the GSR, and all Pc2109 RXLRs in this cluster resided in GSRs and were annotated as SNE1, except for Pc2109RXLR187. Pc2113RXLR058 was upregulated only in N. benthamiana, Pc2113RXLR057 was upregulated in avocado and N. benthamiana, and Pc2113RXLR062 was exclusively upregulated in avocado. This SNE1 clade is closely related to another clade containing several Pc2113 and Pc2109 RXLRs annotated as “unknown proteins”. In this clade, we observed an expansion of RXLRs in both Pc2113 and Pc2109. For Pc2109RXLR125, we observed an expansion of RXLRs in Pc2113 (7). Similarly, we observed a larger cluster of Pc2109 RXLRs (12) with the majority in GSRs (8) and the remaining four in the in-between regions (Figure 7A). Another RXLR clade grouped several Pc2113 and Pc2109 RXLRs with a characterized P. infestans suppressor of plant immunity RXLR (NCBI accession number D0N6D2.1) (Zheng et al., 2014; He et al., 2019). The subclade contained D0N6D2.1 and was composed of one putative homolog in Pc2113 and three homologs in Pc2109. Pc2113RXLR171 from this cluster was upregulated in avocado and N. benthamiana. Associated with D0N6D2.1, we identified another subclade composed of two RXLRs (RXLR004 and RXLR169) annotated as “Avirulence protein 1b” and three RXLRs (RXLR096s and RXLR131) annotated as “Avr1b Avirulence-like protein”, all located in GSRs in both genomes. Pc2113RXLR004 was exclusively upregulated in avocados, and Pc2113RXLR096 was upregulated in all hosts (Figure 7B). We found a subclade that contained RXLRs related to previously annotated P. cinnamomi (NCBI accession numbers QVE55524.1 and QVE55553.1) (Joubert et al., 2021) and P. infestans Avr4 (NCBI accession number ABV66276.1) (van Poppel et al., 2008). This clade was further divided into two subclades: one subclade on top contained three Pc2113 RXLRs and two Pc2109 RXLRs related to QVE55524.1. Two Pc2113 RXLRs (Pc2113RXLR135 and Pc2113RXLR080) were upregulated in zoospores, and Pc2113RXLR135 was also upregulated in avocado at 16 hpi. The bottom subclade was made up of P. cinnamomi RXLRs annotated as ‘Avr4′ with one member, Pc2113RXLR065, upregulated only in avocado and N. benthamiana (Figure 7C). We identified clades that contained previously annotated RXLRs from P. cinnamomi (NCBI accession numbers QVE55566.1, QVE55565.1, and QVE55538.1) (Joubert et al., 2021) and P. infestans (NCBI accession number EEY65678.1) (McLellan et al., 2013). Interestingly, no homologs or few homologs to those RXLRs were found in Pc2113 or Pc2109. One QVE55538.1 homolog was identified in Pc2113 (Pc2113RXLR107) and was upregulated in all samples (Figure 7D). A clade that contained several Pc2113 and Pc2109 RXLRs related to P. sojae Avr1a (NCBI accession number ABO47652.1) was identified; however, none of these RXLRs were annotated as Avr1a homologs. Pc2113RXLR124 and Pc2109RXLR041 were annotated as “Avr1b-1 avirulence-like protein” and were grouped with a previously annotated P. cinnamomi RXLR (NCBI accession number QVE55549.1). Pc2113RXLR124 was found to be exclusively upregulated in avocados. Pc2113RXLR101 and Pc2113RXLR115 were annotated as “avh87”, and only Pc2113RXLR101 exhibited upregulation in avocado and N. benthamiana. Pc2109RXLR35, Pc2109RXLR34, and Pc2113RXLR161 were annotated as “Avh242”, and only Pc2113RXLR161 exhibited upregulation in avocado (Figure 7E).

To confirm the host-specific patterns of the P. cinnamomi effectors exhibiting significant up- or down-regulation in this study, we assessed the overall gene expression profiles in terms of TPM of 18 CAs, 220 CWDEs, 34 elicitins, 22 NIPs, and 81 RXLRs in all samples. Consistent with the DE analyses, global expression indicated that P. cinnamomi expressed (TPM > 5) fewer effectors when infecting A. thaliana and N. benthamiana at 16 hpi when compared with N. benthamiana at 24 hpi and avocado. Interestingly, more effectors were expressed in planta and located in the GSRs when compared with their expression in P. cinnamomi growing in vitro, especially for CWDEs, RXLRs, and elicitins (Supplementary Figure S7 and Supplementary Table S9). These results are consistent with the lesion phenotypes of the pathogen infecting each host (Supplementary Figure S5). We selected three effectors and validated their expression in N. benthamiana leaves infected with P. cinnamomi using qPCR analyses. Consistent with our RNA-seq expression data, Pc2113RXLR047 and Pc2113RXLR52 were highly upregulated at 16 and 24 hpi when compared with the mycelia and sporangia samples from the pathogen growing in vitro. Similar results were found for the elicitin (Pc2113T1C00002373g0000080.1) expression assessed by qPCR, showing significant upregulation at 16 and 24 hpi (Supplementary Figure S8).