RXLR effector genes mediate regional adaptation of Phytophthora infestans

Abstract

Local adaptation has been a central theme of eco-evolutionary research for decades. It is generally assumed that plant pathogens are locally adapted due to their standing interactions with biotic and abiotic factors in the ecosystem. Effectors, secreted small proteins encoded by pathogens, play critical roles in host–pathogen interactions, by activating host genotype-specific resistance, suppressing plant immunity, and playing other functions. In this study, we investigated the potential involvement of RXLR effector genes in ecological adaptation by examining the simple sequence repeat (SSR), virulence, and effector profiles in Phytophthora infestans isolates collected from two geographic regions differing in ecological environments. Genotypic analyses with SSR markers and virulence assays showed that the pathogen from the two regions shared genetic background but differed in virulence spectrums. High-throughput sequencing and expression analysis of 24 selected P. infestans isolates further showed variations in the RXLR effector repertoire, ranging from 536 to 548 for each isolate and the expression of effector genes was highly associated with the accumulation of homologous sRNA. Regional specific alleles were detected at 94 RXLR effector genes, and a specific accumulation of homologous 25–26 nt sRNAs was found at 67 RXLR effector genes. Two of the regional specific RXLR effector genes were confirmed to be virulence factors. Taken together, these results suggest that genomic and epigenetic variations in RXLR effector genes contribute significantly to the ecological adaptation of P. infestans populations and that regional specific effector genes will help to understand the adaptive landscape of pathogens and efficient use of host resistance genes.

Background

Pathogens have widespread impacts on ecosystem structure including the community composition, biodiversity, and species distribution, and through such impacts regulate ecosystem processes including primary production, secondary production, biogeochemical cycles, and disturbance mechanisms (Fischhoff et al. 2020). The contribution of pathogens to ecosystems and their processes begins with their interactions with hosts, altering plant abundance, distribution, and behaviours, as well as community composition (Fischhoff et al. 2020; Mainwaring et al. 2023), which leads to rapid and reciprocal shifts in genetic structure and diversity of both host plants and pathogens and drives their co-evolution (Khavkin 2021).

Over this long-term antagonistic interaction, plant hosts continuously develop a resistance repertoire to defend against pathogens, while pathogens constantly sharpen their invasion capabilities to overcome the defense mechanisms of plants (Ngou et al. 2022). Effector, a group of small proteins secreted by pathogens, play critical roles in host–pathogen co-evolution, involved in overcoming host defense system, nutrient absorption, and ultimately the development of disease symptoms (Uhse and Djamei 2018; Chen et al. 2019). In addition to its role in host–pathogen interaction, effectors also contribute to other biological and ecological functions, such as regulating osmotic pressure (Kasmi et al. 2018), antimicrobial activities (Snelders et al. 2020), and microbiota composition (Cai et al. 2023).

In addition to sequence characteristics, the contribution of effectors to host–pathogen interactions also depends on their expression and successful delivery to host cells. Many effector genes are not constitutively expressed. Their expression and delivery are induced upon contact with the host and involve many biotic and abiotic environments. Abiotically, it has been documented that the expression and delivery of effectors are influenced by climatic conditions, such as air temperature, light spectrum, and humidity (Goel et al. 2008). For example, increased translocation of effectors into host plants was detected during Pseudomonas syringae pv. tomato DC3000 infection in Arabidopsis at elevated temperature (Huot et al. 2017), and enhanced UV-B radiation can down-regulate the expression of disease-causing genes in Magnaporthe grisea (Huang et al. 2019). In P. infestans, both UV radiation and temperature can regulate the expression and evolution of Pi02860 (Yang et al. 2022). The ability of effector proteins AvrE and HopM1 to cause water soaking requires high atmospheric humidity in P. syringae (Xin et al. 2016). Climatic conditions can also alter mutation and selection of effector genes, leading to spatial heterogenicity across different ecosystem (Velásquez et al. 2018). For example, population genetic analyses showed that P. infestans effector genes, such as Avr1 (Shen et al. 2021), Avr2 (Yang et al. 2020), Avr3a (Yang et al. 2018), and Avr4 (Waheed et al. 2021), were strongly associated with annual mean temperature. Under UV radiation, the mutation rate increased significantly in Escherichia coli at the genomic scale (Shibai et al. 2017).

From a biological perspective, effector gene expression can be regulated by internal interactions with other genes in the genome (i.e., epistasis; Phan et al. 2016) or external interactions with metabolites secreted by other species (Bever et al. 2012). For example, as a common regulator of gene expression in eukaryotes (Castel and Martienssen 2013), small RNAs (sRNAs), especially 25–26 nt sRNAs, are associated with silencing of homologous RXLR effector genes in Phytophthora pathogens (Jia et al. 2017; Zhang et al. 2019) and the components in this gene silencing pathway are shared among species (Vetukuri et al. 2011). Co-immunoprecipitation (Co-IP) assays showed that the classes of 20–22 nt and 24–26 nt sRNAs could bind to PiAGO1 and PiAGO4 and modify their translation, respectively, confirming the regulatory function of sRNAs (Åsman et al. 2016). In oomycete, abundant endogenous sRNAs, such as 21 nt and 25 nt sRNAs, were detected in P. infestans, P. sojae, P. parasitica, and P. ramorum (Vetukuri et al. 2012; Qutob et al. 2013; Jia et al. 2017) and silencing of Avr3a (Qutob et al. 2013) and Avr1b (Wang et al. 2019) in P. sojae was confirmed to be associated with homologous sRNA accumulation. Further sequence analysis showed that some flanking regions of effector genes also involved in the epigenetic processes of expression regulation (Wang et al. 2020). Since RNA-seq and sRNA-seq results showed a clear negative correlation between expression of RXLR effector genes and homologous 25–26 nt sRNA accumulation (Jia et al. 2017), the expression levels of effector genes could be directly reflected by their homologous 25–26 nt sRNA accumulation.

Our research on effector genes has mainly focused on functional analysis of their pathogenicity aspects such as virulence/avirulence (Wang et al. 2023) and how they evade host recognition by sequence variation (Avr1, Avr4, Avr3a, Pi20303; Du et al. 2018; Du et al. 2021; Shen et al. 2021; Waheed et al. 2021), but knowledge about how epigenetics and abiotic factors influence evolutionary adaptation of effectors is largely unknown. Especially, we generally lack information on whether the local adaptation of pathogens to ecological niches is associated with regional specific variation in genomic signatures and/or expression of effector genes. Therefore, exploring regional variation in genomic characteristics and expression of effector genes and the contribution of abiotic environments to the variation can provide clues for understanding the role of pathogen effectors in ecological function and resilience. This knowledge is also important in the efficient deployment of resistance genes for sustainable crop production.

In this study, we combined population genomics and potato (Solanum tuberosum)-P. infestans interactions to address these questions. Potato is one of the most important non-grain food crops and plays a crucial role in the world food security system (Guo et al. 2023). However, primary production of potato can be severely threatened by P. infestans, the notorious pathogen that caused the Irish famine in 1840s (Lal et al. 2019). Under the favored temperature and humidity conditions, P. infestans can multiply rapidly and destroy an entire potato crop in a field within days. The pathogen encodes a large number of secreted effector proteins, including the most important RXLR effectors that contain a conserved RXLR (Arg-any amino acid-Leu-Arg) domain (Rehmany et al. 2005). It is estimated that the ~ 240-Mb genome consists of 563 annotated RXLR effector genes, and most of them are located in the rapidly evolving, repeat-rich but gene sparse regions (GSR; Haas et al. 2009), providing unique niches for quick mutation and evolution.

It’s common that the field performance of potato varieties for late blight resistance differs among regions, suggesting regional adaptation of P. infestans populations and the potential role RXLR effector genes may play. To investigate potential involvement of RXLR effector genes in mediating regional adaptation of P. infestans populations, we collected 155 P. infestans isolates from Tianshui, Gansu (TSGS) and Weining, Guizhou (WNGZ), and performed high-throughput sequencing on 24 of these isolates (21 of which were dominant genotypes in both regions). TSGS and WNGZ are among the main potato producing areas in China and differ substantially in ecological environments. TSGS in Northwest China has moderate temperatures and little rainfall, while WNGZ in Southwest China has lower-temperature and more rainfall during the potato planting season. The prevalence of late blight differs significantly between the two regions.

The specific objectives of the study were to (1) investigate the population genetic structure of P. infestans populations in SSR markers and RXLR effector repertoire; (2) identify regional specific RXLR effector genes at genomic and transcriptional levels; and (3) determine ecological and genetic environments contributing the regional differentiation of effector genes.

Results

Shared P. infestans genotypes with variation in virulence in distant regions

A total of 8 genotypes were identified from 155 P. infestans isolates which were collected from potato varieties Longshu 10, Nongtian 1, Tianshu 10, Tianshu 11, Tianshu 12, Tianshu 13, Weiyu 5, Weiyu 3, and Qingshu 9 cultivated in TSGS and WNGZ (Additional file 1: Table S1). The genotypes Gt1 and Gt2 dominated in both regions (Additional file 2: Figure S1a, b). Two (Gt6 and Gt7) genotypes appeared to be regional specific, one in each region. All other six genotypes were detected in both regions, suggesting long-distance transmission of P. infestans.

Virulence assays of 24 P. infestans isolates (12 from each region) on potato cultivar Qingshu 9 showed that 16 of the isolates successfully infected potato leaves (Additional file 1: Table S2, S3). However, symptoms caused by the WNGZ isolates appeared to be more aggressive than those caused by the TSGS isolates, resulting significantly larger lesions (Fig. 1a, b).

Fig. 1

Late blight phenotypes and differential lesion development by P. infestans isolates from TSGS and WNGZ. a Disease phenotype upon infection. A total of 16 isolates (9 from TSGS and 7 from WNGZ) were examined, with lesion diameters induced by the WNGZ isolates significantly larger. The figure shows the disease phenotype of isolate PL1024. b Development of lesion diameters. The lesions were presented as the mean ± SE of at least 3 independent biological replicates. Statistical significance was assessed by Student’s t-test. ns, not significant. ***P < 0.001

Significant variations in P. infestans RXLR effector repertoire

Due to the unsuccessful library construction for 5 P. infestans isolates (PL1007, PL1011, PL1029, PL1038, and PL6080), we performed whole genome re-sequencing for 19 of them, including 8 from TSGS and 11 from WNGZ (Additional file 1: Table S4). RXLR effector repertoire in the isolates ranged from 536 to 548 (Fig. 2a). Among them, 521 annotated RXLR effector genes were mapped to all 19 isolates, while 20 functionally characterized effector genes, including Avr1 and Avr3b, were completely missing from the genomes. Among the 19 genomes, 22 RXLR effector genes showed present or absent variations (PAV; Fig. 2b). However, there is no regional specific PAV in RXLR effector genes between the genomes from two regions.

Fig. 2

RXLR effector repertoire in 19 sequenced P. infestans isolates. a Number of RXLR effector genes in each isolate, including annotated genes and predicted novel genes. b PAV of 22 annotated RXLR effector genes. The colored grids represent the existence of the effector genes in the corresponding isolate. Blue and red represent the isolates from TSGS and WNGZ, respectively. c Conserved amino acid sequences (left) and distribution (right) of novel RXLR effector genes. d Validation of the novel RXLR effector genes. Novel RXLR effector genes were verified by PCR amplification with genomic DNA

We also examined novel RXLR effector genes not annotated in each genome by performing de novo assembly of the unmapped Illumina reads. From the 42.7 Mb contigs that did not map to the reference genome, 28 novel RXLR effector genes were identified (Fig. 2c). Except for Novel_003, most of the novel effector genes also showed PAVs among the sequenced isolates. Eight of the novel effector genes which present in most isolates were confirmed by PCR (Fig. 2d; Additional file 1: Table S5).

RXLR effector genes exhibiting regional differences

Based on the mapping results, a total of 584,869 SNPs were identified among the 19 genomes. Of which 46,405 (7.9%) were enriched in RXLR effector genes (gene body and 1 kb flanking regulatory sequences). Phylogenetic and principal component analysis (PCA) of the SNPs in the whole genomes (Additional file 2: Figure S2a) showed that the TSGS and WNGZ isolates could not be separated according to their geographic origin (Additional file 2: Figure S2b–d), further suggesting close genetic relationship between P. infestans from the two areas. However, phylogenetic analysis of SNPs from RXLR effector genes showed that P. infestans isolates tended to be clustered according to their geographic origins, with some exceptions such as PL1015 (TSGS isolate) and PL6111 (WNGZ isolate; Additional file 2: Figure S2e).

To identify variant sites that represent the differences among 19 P. infestans isolates, the consistent SNPs/Indels among all isolates were filtered out. 644 non-synonymous SNPs and 84 InDels (at least one individual was different with the other 18) were identified in the coding regions of 244 and 63 RXLR effector genes, respectively (Fig. 3a; Additional file 2: Figure S3a). Among them, no variant sites that were present in all isolates of one region but absent in all isolates of another region were detected. After calculating the number of SNPs/Indels in the isolates (Additional file 1: Table S6), we considered a variant that appeared at least three isolates in one region but not in another region to be regional specific. Accordingly, 33 SNPs and 4 InDels showed regional specific (Fig. 3c; Additional file 2: Figure S3c). For instance, the known effector gene Avr-Smira2 exhibited consensus SNPs in six of nine genomes in the TSGS isolates, resulting in a Cly to Arg amino acid change at position 55 (Additional file 2: Figure S3e).

Fig. 3

Sequence variations in P. infestans RXLR effector genes. a, b Number of SNPs (left) and InDels (right) in coding regions and promoter regions of effector genes in TSGS and WNGZ isolates. The number above the black line represents the total SNPs and InDels, which is the union of SNPs and InDels in each isolate. c, d Number and distribution of regional specific SNPs in coding regions and promoter regions. Only part of the whole picture was taken for display. The histogram shows the number of SNPs, and the colored grids represent the existence of SNPs in the corresponding isolates

Analyses of sequence variation in the promoter regions of the RXLR effector genes showed 1650 (0.28%) SNPs and 521 (0.68%) InDels in the regulatory sequences of 385 and 251 RXLR effector genes (Fig. 3b), respectively. The number of differential SNPs and Indels in each genome ranged from 648–860 to 247–310, respectively (Additional file 2: Figure S3b). Among them, we identified 32 and 40 regional specific mutations between the TSGS and WNGZ isolates, respectively (Fig. 3d; Additional file 2: Figure S3d), including the known effector genes Avr-Smira1, Avr-amr1, and Pi22798 (Additional file 2: Figure S3f). Further analyses confirmed that the transcript levels of WNGZ-specific Pi12721, Pi23048, Pi15109, and TSGS-specific Pi19302 exhibited significant regional signatures during plant infection. Pi12721 and Pi15109 were highly expressed in the WNGZ isolates at 12 hpi and 48 hpi, respectively (Additional file 2: Figure S4a, c). Pi19302 was continuously expressed in the WNGZ isolates during infection, and its expression level was significantly higher than that in the TSGS isolates (Additional file 2: Figure S4d). Pi23048 displayed low expression only in the mycelia of the TSGS isolates (Additional file 2: Figure S4b). These results suggest that the contribution of regulatory sequences to the expression variation of regional specific effectors is likely limited.

The Indel mutations identified in RXLR effector genes ranged from 1 to 51 bp in length, with the majority (86.15%) shorter than 15 bp. The 1-bp and 10-bp InDels were more frequent in the flanking regulatory sequences of RXLR effector genes (Additional file 2: Figure S5a). Further analysis showed that the 10-bp InDels were widely distributed (14%) in the regulatory sequences of all genes in P. infestans genome (Additional file 2: Figure S5b), and their frequency was similar to that in the RXLR effector genes (12%). However, there was no regional specific 10-bp InDels in the promoter regions of these effector genes.

Expression variation in RXLR effector genes

We analyzed the accumulation of 25–26 nt sRNAs homologous to RXLR effector genes during vegetative growth in 24 P. infestans isolates by sRNA-seq (Additional file 1: Table S7) and performed transcriptome sequencing of RXLR effector genes (Additional file 1: Table S8). The results showed that the expression of most RXLR effector genes (~ 70%, 390–425 genes) was repressed during the vegetative growth stage (Fragments per Kilobase Million, FPKM < 1) (Additional file 2: Figure S6a). Among the 563 annotated RXLR effector genes, 223 (39.6%)–369 (65.5%) were detected with homologous 25–26 nt sRNAs (RPKM > 1) in each P. infestans isolate (Fig. 4a; Additional file 2: Figure S6b). Among them, the expression of 197–298 (79.8–88.3%) RXLR effector genes was completely silenced (Fig. 4a).

Fig. 4

Regional specific RXLR effector genes enriched with homologous 25–26 nt sRNA. a Number of effector genes enriched with sRNAs. Blue/light blue and red/light red represent the number of silent/expressed effector genes from TSGS and WNGZ isolates, respectively. b Number and distribution of regional specific effector genes. The figure above the bar graph represents the number of genes, and the colored and white grids below represent the present or absent accumulation of sRNAs at these genes in the isolate. c Heatmap visualization of sRNA accumulation levels of regional specific effector genes. d–i Differential accumulation of 25–26 nt sRNAs homologous to effector genes Pi12046, Pi09647, Pi15225, Pi06083, Pi22935, and PexRD2 in TSGS and WNGZ isolates. Statistical significance was assessed by Student’s t-test. *P < 0.05, **P < 0.01

To identify differential expression of RXLR effector genes between the two regions, we examined variations in sRNA accumulation that are detected in at least three isolates from the same region. We identified 277 and 295 sRNA-accumulating RXLR effector genes in TSGS and WNGZ isolates, respectively, with 260 being shared in both regions (Additional file 2: Figure S6c). There were 83 RXLR effector genes that accumulated with the 25–26 nt sRNAs across all isolates (Fig. 4b). Based on their accumulation with homologous 25–26 nt sRNAs in vegetative hyphae, it was notably found that 5 and 10 RXLR effector genes were inferred to be specifically expressed in TSGS and WNGZ isolates, respectively (Fig. 4b). For instance, the most regional specific genes including the TSGS-specific Pi22740 and Pi01905, and the WNGZ-specific Pi07566 showed abundant accumulation of the 25–26 nt sRNAs in isolates from one region but rare from another region. In addition to the qualitative difference in presence or absence, we also found quantitative differences in sRNA accumulation levels between the two regions. In the TSGS isolates, 18 RXLR effector genes accumulated more homologous 25–26 nt sRNAs, whereas 13 accumulated more homologous 25–26 nt sRNAs in the WNGZ isolates (Fig. 4c). For example, Pi12046, Pi09647, Pi15225, and Pi06083 accumulated more 25–26 nt sRNAs in the TSGS isolates, but Pi22935 accumulated more 25–26 nt sRNAs in the WNGZ isolates (Fig. 4d–h). PexRD2 family is a class of RXLR effector genes (Haas et al. 2009), in which five members significantly accumulated more sRNAs in the WNGZ isolates than that in the TSGS isolates (Fig. 4i). Further expression analysis suggested that the high-level accumulation of homologous 25–26 nt sRNAs at Pi09647 and Pi12046 in the TSGS isolates, as well as Pi22935 and Pi18670 in the WNGZ isolates, were negatively associated with the expression of the corresponding RXLR effector genes throughout or at specific stages of infection (Additional file 2: Figure S7a–d).

Enriched accumulation of homologous 25–26 nt sRNA at flanking regulatory sequences of RXLR effector genes

Except for the coding regions, we also found obvious accumulation of 25–26 nt sRNAs at the flanking regions of RXLR effector genes. Overall, there was no significant difference in the number of RXLR effector genes that accumulate homologous 25–26 nt sRNAs in the coding regions, upstream and downstream regulatory sequences (Fig. 5a). However, the abundance of homologous 25–26 nt sRNAs in the upstream and downstream regions was significantly higher than that at coding regions in each isolate (Fig. 5b; Additional file 2: Figure S8a), indicating a higher degree of 25–26 nt sRNA accumulation at the flanking regions of RXLR effector genes.

Fig. 5

Accumulation of homologous 25–26 nt sRNAs in the promoter region of RXLR effector genes. a Number of effector genes with sRNA accumulation in the upstream, downstream, and coding regions. b Average accumulation level of sRNA homologous to effector genes. c Number of effector genes enriched with sRNAs in the promoter regions. d Number and distribution of regional specific effector genes. e Heatmap visualization of sRNA accumulation levels of regional specific effector genes. f–k Differential accumulation of homologous 25–26 nt sRNAs in effector genes Pi15225, Pi12046, Pi06083, Pi09647, Pi10905, and Pi04329 in TSGS and WNGZ isolates. Statistical significance was assessed by Student’s t-test. ns, not significant. *P < 0.05

We further examined the accumulation of 25–26 nt sRNA homologous to the upstream regulatory sequences of RXLR effector genes and compared it with the gene expression level in each isolate at the vegetative growth stage. The results showed that 227–345 RXLR effector genes accumulated homologous 25–26 nt sRNAs in the promoter region and 74.7–83.0% of them were silenced (Fig. 5c). We also found that four and five regional specific RXLR effector genes differentially accumulated homologous 25–26 nt sRNAs in TSGS and WNGZ isolates, respectively (Fig. 5d). It’s also notable that 24 and 3 RXLR effector genes accumulated more homologous 25–26 nt sRNAs than others in TSGS and WNGZ isolates (Fig. 5e), respectively. The TSGS-specific RXLR effector genes Pi15225, Pi12046, Pi06083, Pi04329, Pi01905, and Pi09647 accumulated more homologous 25–26 nt sRNAs in both promoter and coding regions (Fig. 5f–k). Further exploration revealed that there was no correlation between sequence variations in the promoter region and the accumulation of sRNA.

RXLR effector genes differentially accumulated with homologous 25–26 nt sRNAs were virulence factors

To examine whether the sRNA-associated RXLR effector genes contribute to the virulence of P. infestans, the regional specific RXLR effector genes Pi22740 and Pi04329, and region-shared effector genes Pi18487, Pi08399, Pi04350, and Pi17214, all with significant accumulation of homologous 25–26 nt sRNAs, were cloned and transiently expressed in Nicotiana benthamiana by agroinfiltration assays. After 5 days post-inoculation with P. infestans zoospores, N. benthamiana leaves transiently overexpressed with RXLR effector genes all developed significantly larger lesions than the GFP control (Fig. 6a, b). The results showed that the two regional specific RXLR effector genes differentially accumulated homologous 25–26 nt sRNAs were virulence factors, suggesting their contribution to regional specific pathogenicity of P. infestans.

Fig. 6

Transient overexpression of RXLR effector genes enhanced plant susceptibility to P. infestans infection. a Phenotypes of N. benthamiana leaves infected with P. infestans. GFP was transiently expressed on the left of the detached leaves, and the effector was transiently expressed on the right. Photos were taken at ~ 7 days post inoculated with P. infestans zoospores. b Average lesion area of leaves expressing GFP and the effectors. The lesions were presented as the mean ± SE of at least 20 independent biological replicates. Statistical significance was assessed by Student’s t-test. ***P < 0.001

Discussion

Identical genotypes were captured in isolates sampled from two areas approximately 1200 km apart, consistent with the high and long gene flow reported in the pathogen (Cooke et al. 2012). Gene flow homogenizes a species’ genetic variation that would otherwise accumulate due to geographic isolation (MacDonald et al. 2020). In P. infestans, gene flow can occur naturally through the movement of propagules or anthropogenically through the exchange of plant materials (Gao et al. 2020). For example, based on detecting the same SSR genotypes as we used, Tian et al. (2016) documented that gene flow of the pathogen can also occur among distant P. infestans populations.

Conflicting to the population genetic structures in SSR markers, significant spatial population differentiation was detected between the two pathogen samples in RXLR effector genomes, consistent with hypothesis of enhancing evolution in genomes involved in host–pathogen interactions and results from previous studies in P. infestans (Yang et al. 2020; Shen et al. 2021; Waheed et al. 2021) and other pathogens (Mandal et al. 2022). RXLR effectors of P. infestans are mainly located in highly variable regions of the genome surrounded by transposable elements (Haas et al. 2009), allowing the pathogen to generate large amounts of variation including continuous detection of novel effectors and rapidly adapt to the changing biological environments such as resistant genes developed in the plant host (Rietman et al. 2012). Cooke et al. (2012) performed de novo assembly and revealed six candidate RXLR effector genes that were not present in the reference genome. Using the same methods, we identified a total of 28 novel RXLR effector genes in the sequenced isolates (Fig. 2c). Interestingly, 27 novel RXLR effector genes are present in various combinations in all isolates, suggesting that recombination may also contribute to the generation of genetic variation in the effector genes of P. infestans. In terms of disease management in primary food production, this result highlights the importance of exploring the spatial distribution of RXLR effector genes according to ecological context and well explains the lasting efforts in this research field (Gilroy et al. 2011).

We found that 34 and 62 RXLR effector genes contain regional specific variations in the coding or promoter regulatory regions (Fig. 3c, d; Additional file 2: Figure S3c, d). Polarized population differentiation of SSR marker loci and effectoromics as well as the regional specific variation in RXLR effector repertoires also provide strong evidence for local adaptation of the functional genes, contributing to the differences in the frequency and severity of late blight epidemics between WNGZ and TSGS (Todd et al. 2022). Different potato varieties were used in the two regions (Additional file 1: Table S1). These varieties may have different R genes, thereby contributing to the observation of regional differentiation. However, typically only a limited number of effectors are involved in resistance in a given variety (Vleeshouwers and Oliver 2014), and the contribution of this host variation to the observations, if any, is expected to be small. Furthermore, similar spatial differentiation in effector genes was also found when we only used the pathogen genome from potato variety Qingshu 9 grown in both regions.

Although host effects cannot be completely ruled out, we believe that ecological environments, especially climatic factors, also play an important role in the development of population differentiation in pathogen effector genes. Climatic conditions such as temperature, humidity, and ultraviolet radiation vary greatly between the two potato producing areas, which can generate a multifaceted impact on the function and evolutionary development of effector genes, and thus on the adaptation and epidemics of the pathogen (Cheng et al. 2019). Genetically, climatic conditions can directly modulate the generation and preservation of genetic variation in effector genes by altering mutation rates. It can indirectly affect the generation and maintenance of genetic variation through changes in pathogen population size and reproductive strategies and efficiency. For example, temperature can affect the thermal stability of DNA molecules, thereby affecting the mutation rate of the genome including effectoromics during replication (Marguet and Forterre 1994), while ultraviolet radiation is a well-known mutagen that can induce mutations in effectors and other genes (Rastogi et al. 2010), generating spatial population differentiation of effectors along temperature gradients and altitudes in P. infestans (Yang et al. 2022) and many other pathogens (Matsuba et al. 2013). Epigenetically, climatic conditions can also moderate the expression of effectors directly such as by turning on/off promoters (Singh et al. 2014), or indirectly by regulating the production of translational regulators such as homologous sRNAs (Zhang et al. 2016), affecting their functional efficiency during host–pathogen interactions.

Indeed, we found strong evidence that gene expression polymorphisms contribute to the functional diversity of effector genes and potential pathogen adaptation to ecological niches and their expression is associated with accumulation of homologous sRNAs. In P. infestans (Vetukuri et al. 2012), as well as in P. sojae (Qutob et al. 2013), 25–26 nt sRNAs can silence RXLR effector, and we found that most (79.8–88.3%) of the RXLR effector genes that accumulate homologous 25–26 nt sRNAs in our study were unexpressed at mycelium stage (Fig. 4a), which enables the pathogen to produce effector proteins as needed to maximize energy efficiency (Xu et al. 2022). Furthermore, the accumulation of these homologous sRNAs in effector genes was regional specific. Although we identified 21 regional specific sRNA-associated RXLR effector genes in both the TSGS and WNGZ isolates, the total amount of sRNAs homologous to RXLR effector genes in the TSGS isolates was significantly higher than that in the WNGZ isolates (Additional file 2: Figure S8b). We believe that this differential pattern of sRNA accumulation between the pathogen isolates from the two regions may also exist during the infection stage. Due to the negative association between sRNA accumulation and pathogenicity (Qutob et al. 2013; Wang et al. 2019), this result indicates that the pathogen from TSGS causes less disease than that of WNGZ. On the other hand, it is interesting to find that the total expression levels of RXLR effector genes were not significantly different between the two regions (Additional file 2: Figure S8c), indicating that some RXLR effector genes in the TSGS isolates were repressed more heavily by sRNAs, which possibly contribute to less disease in TSGS. Pi21422, a PexRD2 family member, accumulated more sRNAs in WNGZ than TSGS isolates (Fig. 4i). As a virulence factor, Pi21422 interacts with the kinase domain of MAPKKKε to perturb plant defense signaling (King et al. 2014) and higher accumulation of sRNA in Pi21422 may enhance the ability of the pathogen to evade host recognition, thereby reducing host defense in potatoes grown in WNGZ. Taken together, these results may help to explain the reduced severity of late blight in TSGS compared with WNGZ.

Promoters are known to activate gene transcription through the recognition by RNA polymerase. However, sRNAs can mimic the promoter and bind to the σ subunit of RNA polymerase, thereby altering the activity of RNA polymerase and affecting the transcription process (Caswell et al. 2014). In addition, sRNAs can also selectively interfere with gene promoter regions and accidentally activate gene transcription (Li et al. 2005), indicating that sRNAs can mediate gene silencing at the transcriptional level. We detected a large number of homologous sRNAs accumulating in the coding, upstream, and downstream of RXLR effector genes, even more sRNAs accumulation in promoter regions than that in the coding regions (Fig. 5a, b), and the homologous sRNA accumulation levels were also negatively correlated with the expression levels of RXLR effector genes, suggesting that significant sRNA enrichment is characteristic of RXLR effector genes in P. infestans and that epigenetics plays at least a similar role to genetics in regulating effector function.

The accumulation of homologous sRNAs within promoter regions of effector genes was regional specific. In total, 28 and 8 RXLR effector genes with significantly higher sRNAs accumulation were identified in TSGS and WNGZ isolates, respectively, and the expression of more RXLR effector genes may be inhibited in the TSGS isolates can partly explain their weaker virulence. Besides, ~ 50% of effector genes with different sRNA accumulation were consistently found in both gene bodies and regulatory regions (Additional file 2: Figure S8d). Taken these results together, we propose that sRNAs synthesis in the gene bodies and flanking regions of many effectors may be continuous and determined by a consistent epigenetic mechanism in the chromatin regions that regulate the transcription of the sRNA precursor.

Currently, 26 RXLR effector genes with virulence/avirulence have been reported. Among them, only the virulence factor Pi21422 has been identified as a sRNA-associated effector gene. Regardless of whether it was in the gene region or in the promoter regulatory region, the homologous sRNA accumulation levels of these effector genes with known functions were consistently lower than those of sRNA-associated effector genes (Additional file 2: Figure S9). Taking the well-known Avr3a as an example, its homologous sRNA accumulation level was significantly lower than the average accumulation levels of homologous sRNAs in both the sRNA-associated effector genes and entire effector genes (Additional file 2: Figure S9c, f). These results indicate that there may be some new effector genes in the tested isolates, and their differential expression is the cause of virulence variation of P. infestans isolates among regions. The finding of regional specific sRNA-associated effector genes provides both resource of discovering new effectors and clues for understanding the virulence variation, regional adaptation of P. infestans, as well as the effective utilization of resistant potato varieties.

Conclusions

We found spatial variations of RXLR effector genes at the genomic and transcriptional levels in P. infestans from TSGS and WNGZ even though extensive gene flow occurred between the two ecological regions, suggesting that local adaptation in P. infestans is associated with regional variations of RXLR effector genes. However, the underlying functional mechanisms of these eco-specific effector genes remain unknown but has important implications in understanding the adaptive landscape of pathogens such as P. infestans and for the selection and rational deployment of late blight resistance genes for sustainable disease management. This issue deserves further exploration.

Methods

Sampling and purification of P. infestans isolates

P. infestans isolates used in this study were collected from five sampling points in potato fields located at TSGS and WNGZ during the early stages of late blight onset according to the ‘Zigzag’ approach. The isolates were induced from infected potatoes with white mold on the back of leaves on potato tuber as described previously (Tian et al. 2016) and were grown on rye sucrose agar (RSA) medium plates at 16–18°C for 9–12 days. After 2–3 times of purification, the isolates were stored in liquid nitrogen until use.

Plant growth conditions and P. infestans infection assays

S. tuberosum plants were grown in a growth chamber at 23°C with a 16/8 h light/dark cycle. Potato leaves for seven-week-old potato were detached from the plants inoculated on the abaxial surfaces with a 16 μL droplet of spore suspension containing 800 zoospores. The inoculated leaves were incubated at 16–18°C for 3–6 days and lesion diameters were measured.

DNA extraction and SSR analysis

P. infestans isolates grown in RSA broth for 12–14 days at 16–18°C were harvested. Mycelia were quickly frozen with liquid nitrogen and stored at − 80°C. DNA was extracted by Phytophthora DNA extraction solution (7 M Urea, 350 mM NaCl, 50 mM Tris-HCI, 4 mM EDTA-2Na, 1% Na-lauryl sarcosine, 1% Tris-Phenol), and the qualified DNA was stored at − 20°C.

Five polymorphic SSR markers (SSR4, SSR6, SSR8, Pi4B, and Pi02; Lees et al. 2006; Li et al. 2012) were used to genotype the P. infestans isolates according to the protocol described previously (Kiiker et al. 2018). PCR products were separated by electrophoresis in polyacrylamide nondenaturing gel. After silver staining (1% AgNO₃), the gels were developed in 1.5% NaOH solution containing 0.1% formaldehyde until bands were clearly visible and repeatedly washed with distilled water and photographed. Band size was measured by DNA marker.

High-throughput sequencing

The library construction and sequencing of P. infestans genomic DNA were performed using the standard Illumina protocol by Biomarker Technologies (Beijing, CHN). Briefly, the qualified DNA samples were randomly broken into 350 bp fragments by Covaris, and TruSeq Library Construction Kit (Illumina, USA) was used to construct library. DNA fragments were sequenced on Illumina HiSeq 2500 platform.

Total RNA was extracted from the mycelia using TRIzol reagent according to the manufacturer’s protocol (Invitrogen, USA). RNA-seq and sRNA-seq libraries were performed with PE150 on Illumina HiSeq Xten platform by Novogene (Beijing, CHN) with three biological replicates. RNA libraries and sRNA libraries were prepared using 1 µg of total RNA according to the protocols of RNA Library Prep Kit (New England Biolabs, USA) and Small RNA Library Prep Kit (Lexogen, Austria), respectively. For the construction of sRNA library, only 18–45 nt sRNAs were used for sequencing. The raw sequencing data was deposited in the Sequenced Read Archive (BioProject ID: PRJNA1059044).

Bioinformatic analysis

For whole genome resequencing data, quality control and filtering of raw data were performed using FastQC (v 0.11.6) and Trimmomatic-0.36 (-phred33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36), respectively. The clean reads were mapped to the reference genome of P. infestans (T30-4; v 1. 43) by BWA (v 0.7.15) with a seed length of 38 and a maximum of mismatches allowed of 1 (mem -t 4 -M -R). PCR duplications of each sample were removed from alignments using Picard (v 1.119). After indexing with Samtools (v 1.3.1), we used the HaplotypeCaller tool from GATK (v 4.0) software to call variants and obtain gvcf files. The gvcf files from each sample were consolidated by the CombineGVCFs tool, ultimately resulting in a merged gvcf file. SNPs and InDels were extracted using the SelectVariants tool. Variant sites that represent the differences among the 19 strains were obtained by filtering out the SNPs/InDels that were consistent across all isolates based on the merged gvcf file. Picard (v 1.119), Samtools (v 1.3.1), and GATK (v 4.0) were used for PCR deduplication and SNP/InDel indexing and analysis, respectively. The upstream 1 kb sequence of genes was obtained by TBtools and used as the promoter regulatory sequences of the genes. To predict novel RXLR effector genes, unmapped reads of each isolate were extracted from genome alignments in SAM format by Samtools (v 1.3.1) were assembled using Velvet (v 1.2.10) with a Kmer of 30 to obtain contigs (-cov_cutoff 30 -ins_length 300 -ins_length_sd 100 -min_contig_lgth 200), and were then mapped to the reference genome using NUCmer (v 3.23; -c 200 -g 200 -p). After filtering out the contigs that were equivalent to the reference genome, the remaining contigs were retained for the search of open reading frames through an online website http://emboss.bioinformatics.nl/cgi-bin/emboss/getorf. Signal peptides and transmembrane domain of the candidate effectors were predicted using SignalP 4.1 (https://services.healthtech.dtu.dk/services/SignalP-4.1/) and TMHMM 2.0 (https://services.healthtech.dtu.dk/services/TMHMM-2.0), respectively. The detailed protocol for effector searching was as described by Cooke et al. (2012).

For RNA-seq data, raw sequencing reads were available through FastQC (v 0.11.6) and Trimmomatic-0.36. After removing adaptors and low-quality reads, the clean data were mapped to the reference genome of P. infestans using Hisat2 (v 2–2.0.3; -min-intronlen 20 -max-intronlen 3000). Gene expression was calculated using StringTie (v 1.3.3b; -e -B) and the expression level of genes was normalized by FPKM.

Cutadapt (v 2.9) was used to remove adapters from raw data and filter sRNAs with a length distribution between 18 and 45 nt (-m 18 -M 45). After trimming adaptors, the clean data were exactly mapped to P. infestans reference genome using Bowtie (v 1.1.2). No mismatch was allowed and all alignments for each read were reported (-f -v 0 -a -S). The sRNAs derived from ribosomal DNA (rDNA), mitochondrial DNA (mtDNA), and transfer RNA (tRNA) were identified and filtered out by Bowtie (v 1.1.2), allowing for a mismatch (-f -v 1 -a). The retained sRNAs were used for further analyses. Based on the alignment result and the genome annotation information, the sRNA sequences with start positions felling within the gene region of candidate effectors, and their number was counted and normalized with RPKM (Jia et al. 2017).

Transient agroinfiltration assays

The selected P. infestans RXLR effector genes were amplified by primers (Additional file 1: Table S5) linked to the homology arms of the EcoRI and XbaI restriction sites and then ligated into the pART27-NGFP vector using homologous recombinase. For ligation, Agrobacterium tumefaciens strain AGL1 containing plasmid constructs was grown for 24 h in LB medium at 28°C. The A. tumefaciens cells were collected by centrifugation and resuspended in infiltration buffer (10 mM MES, 10 mM MgCl₂, 0.2 mM acetosyringone, pH = 5.6). After adjusting the A. tumefaciens suspensions (OD₆₀₀ = 0.3), Agro-infiltration assays were performed on leaves of 6 to 8-week-old N. benthamiana plants. After 24 h, leaves were detached and 800 zoospores (16 μL) of P. infestans were inoculated at the infiltration sites (Tian et al. 2016). Symptom development was monitored visually from 2 to 8 days after infiltration.

Quantitative RT-PCR analysis

Complementary DNA was synthesized using PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa) according to the manufacturer’s instructions. RT-qPCRs were performed with FastStart Universal SYBR Green Master (Roche). The P. infestans gene PiUBC (Pi08327) was used as the reference gene. Gene expression relative to the reference gene was quantified using the 2^−ΔΔCt method. All specific primers used for RT-qPCR were listed in Additional file 1: Table S5.