Evidence for intragenic recombination and selective sweep in an effector gene of Phytophthora infestans.

Effectors, a group of small proteins secreted by pathogens, play a critical role in the antagonistic interaction between plant hosts and pathogens through their dual functions in regulating host immune systems and pathogen infection capability. In this study, evolution in effector genes was investigated through population genetic analysis of Avr3a sequences generated from 96 Phytophthora infestans isolates collected from six locations representing a range of thermal variation and cropping systems in China. We found high genetic variation in the Avr3a gene resulting from diverse mechanisms extending beyond point mutations, frameshift, and defeated start and stop codons to intragenic recombination. A total of 51 nucleotide haplotypes encoding 38 amino acid isoforms were detected in the 96 full sequences with nucleotide diversity in the pathogen populations ranging from 0.007 to 0.023 (mean = 0.017). Although haplotype and nucleotide diversity were high, the effector gene was dominated by only three haplotypes. Evidence for a selective sweep was provided by (i) the population genetic differentiation (GST) of haplotypes being lower than the population differentiation (FST) of SSR marker loci; and (ii) negative values of Tajima's D and Fu's FS. Annual mean temperature in the collection sites was negatively correlated with the frequency of the virulent form (Avr3aEM), indicating Avr3a may be regulated by temperature. These results suggest that elevated air temperature due to global warming may hamper the development of pathogenicity traits in P. infestans and further study under confined thermal regimes may be required to confirm the hypothesis.

INTRODUCTION

Pathogens and their host plants are engaged in a never‐ending battle, in which pathogens are continuously evolving invasive mechanisms to break plant defense systems and host plants responding with the constant development of new protection to prevent or mitigate damage (Zhan, Thrall, & Burdon, 2014; Zhan, Thrall, Papaïx, Xie, & Burdon, 2015). Many plant pathogens interact with their hosts following the widely accepted gene‐for‐gene model. This interaction model hypothesizes that plant immune reactions are triggered by effectors, a group of small proteins secreted by pathogens during infection, when they are recognized by receptor proteins produced by corresponding host resistance genes (Hein, Gilroy, Armstrong, & Birch, 2009; Jones & Dangl, 2006). This theory suggests that effector genes may reduce pathogen fitness and therefore are selected against during pathogen evolution. However, recent molecular and functional analyses indicate that many effector proteins are essential components of infection and may manipulate plant cellular processes to increase host susceptibility and suppress host immune systems (Boevink et al., 2016; Wang et al., 2015; Yang, McLellan, et al., 2016). Consequently, effector genes are expected to evolve differently to other parts of pathogen genomes in reflecting this antagonistic coevolution between the pathogen and host and the trade‐off of the advantages and disadvantages any effector change may cause to the pathogen population. On the one hand, effector genes might be expected to evolve at an escalated rate in order to effectively and quickly escape host recognition systems and increase the invasive, survival, and reproductive chances of the pathogen. However, at the same time they are potentially constrained by the need to retain features critical to the biological and ecological adaptation of pathogens.

Effector genes commonly lie in gene‐sparse, transposon‐rich regions of the pathogen genome (Haas et al., 2009; Raffaele & Kamoun, 2012). This physical location provides effector genes a unique opportunity to generate sequence duplication. Functional redundancy encoded by multiple copies of effector genes relaxes selective pressure on one or more of the gene copies, which in turn allows more frequent and abrupt mutations to occur without severe impact on the fitness of the pathogen. Indeed, multiple mutation mechanisms including base substitution, deletion, pseudogenization, and transcriptional silencing have been documented in effector genes (Cooke et al., 2012; Raffaele & Kamoun, 2012). Interaction of these mutation events results in higher genetic variation in effector genes than in the rest of pathogen genome (Karasov, Horton, & Bergelson, 2014; Raffaele & Kamoun, 2012) and enhances response to selection driven by the change in host defense systems. Consequently, virulent types of effector genes can emerge quickly from avirulent types in the plant pathogen populations, leading to the breakdown of many major gene‐mediated plant resistances within a few years of their commercialization (Cooke et al., 2012; Pilet, Pellé, Ellisseche, & Andrivon, 2005).

In addition to their genetic characters, other biological and ecological processes may also influence genetic variation and the evolutionary trajectory of effector genes in nature. In a biological context, genetic variation of effector genes can be affected by recombination occurring either within a single gene (intragenic recombination) or between genes (intergenic recombination). Intragenic recombination can increase genetic variation of effector genes directly through generating new sequences or indirectly by reducing allelic loss associated with random genetic drift and hitchhiking selection. In an ecological context, selection driven by host defenses may be confounded by other biological and environmental factors such as temperature in determining the genetic variation and potential evolutionary landscapes of effector genes. Temperature can exert a critical influence on all aspects of chemical, biological, ecological, and evolutionary processes. It can affect survival, reproduction, and transmission of plants and pathogens both individually and interactively. Indeed, temperature‐mediated rate of spontaneous mutations has been documented in some species (Matsuba, Ostrow, Salomon, Tolani, & Baer, 2013). In host–pathogen interaction, temperature can regulate plant resistance responses to pathogens as well as the expression and competition of effector proteins (Banta, Bohne, Lovejoy, & Dostal, 1998).

Phytophthora infestans (Mont) de Bary is a pathogenic oomycete causing late blight diseases of potatoes and tomatoes. It distributes globally and can cost billions of dollars economic losses annually in the potato industry alone (Fry, 2008). In the potato–P. infestans interaction, resistance in the host is triggered by corresponding effector proteins released by the pathogen. Molecular genomic analyses confirm that P. infestans secretes a wide range of effector proteins into extracellular space (apoplastic effectors) or cell cytoplasm (cytoplasmic effectors). The pathogen genome contains at least 550 effectors all earmarked with a conserved N‐terminal motif of Arginine–X–Leucine–Arginine (RXLR), where X can be one of any 20 amino acids (Haas et al., 2009). However, this RXLR effector repertoire is likely expanding, contributing to the enormous size of the P. infestans genome (Haas et al., 2009; Raffaele & Kamoun, 2012). Biochemical and virulent functions of some RXLR effectors in P. infestans have been characterized recently (Bozkurt et al., 2011; King et al., 2014; McLellan et al., 2013; Wang et al., 2015; Yang, McLellan, et al., 2016).

Avr3a, recognized by the corresponding resistance gene R3a in the host plant, was the first effector gene characterized in P. infestans (Armstrong et al., 2005). It is a RXLR effector containing conserved W, Y, and L motifs in C‐terminal domains. These conserved motifs are believed to be generated by purifying selection imposed by the plant host over the course of the co‐evolutionary history of host and pathogen (Dou et al., 2008; Haas et al., 2009; Jiang & Tyler, 2008). The effector gene is also essential to pathogen infection. Silencing Avr3a compromised the pathogenicity of P. infestans in a susceptible potato cultivar and Nicotiana benthamiana (Bos et al., 2010; Vetukuri et al., 2011). Although both virulent and avirulent types of Avr3a existed in P. infestans, the former are more frequent, accounting for 60%–100% of Avr3a in natural populations (Armstrong et al., 2005; Cárdenas et al., 2011; Wu et al., 2016). Recently, several main groups of isoforms have been functional analyzed. Avr3aEM, containing amino acid residues S19, E80, M103, can evade recognition by plant hosts carrying R3a resistance gene (Bos et al., 2006) and is the most dominant virulent type. Avr3aKI, containing amino acid residues C19, K80, and I103, is an avirulent type. This isoform can be recognized by potato hosts with R3a resistance; it also promotes pathogen pathogenicity by suppressing the programmed cell death of resistant host plants. Avr3aKI/147del and Avr3aKI/Y147S all demonstrate the ability to activate the hypersensitive response in plants carrying R3a resistance but do not enhance the ability of the pathogen to infect (Bos, Chaparro‐Garcia, Quesada‐Ocampo, Gardener, & Kamoun, 2009).

Research on effector genes mushroomed in the last decade. The majority of these studies have focused on the molecular and functional dissection of effector genes. Population genetic analyses of effector genes and their interactions with environmental factors such as temperature are limited but urgently required to properly understand the evolutionary and epidemic behaviors of plant pathogens under a global warming scenario. In this study, we used Avr3a in P. infestans as a model to test the hypothesis that intragenic recombination and thermal conditions play an important role in the evolution of effector genes in plant pathogens. The specific objectives of this study were as follows: (i) to determine evolutionary mechanisms generating genetic variation of Avr3a; (ii) to infer the contribution of natural selection to the spatial population genetic dynamics of Avr3a; and (iii) to evaluate the potential impact of temperature on the spatial distribution and evolution of Avr3a.

MATERIALS AND METHODS

Pathogen isolates

Ninety‐six isolates (Table 1) each with a distinct genotype determined previously by molecular and phenotypic markers (Zhu et al., 2015; Yang, Zhu, et al., 2016) were included in this study. The isolates were collected from six commercial fields, one each from six regions located in Anshun (Guizhou), Fuzhou (Fujian), Guyuan (Ningxia), Kunming (Yunnan), Nanning (Guangxi), and Yangshuo (Ningxia) in the early stage of epidemics between 2010 and 2012. Gansu and Ningxia populations were kindly provided by Professor Weixing Shan at the Northwest A&F University. Detailed information on pathogen collection and isolation can be found in previous publications (Yang, Zhu, et al., 2016; Zhu et al., 2015). Briefly, infected leaf samples (one per plant) were collected randomly from potato plants at a spacing of at least 100 cm. Infected leaves were kept separate to prevent cross‐contamination and transferred to the laboratory on ice within 24 hr. In the laboratory, infected leaves were rinsed with sterilized water and single‐mycelium isolates secured from sporulating lesions using an inoculating needle. The isolates were purified by repeatedly transferring a single piece of mycelium to fresh medium three times. Only one single‐mycelium isolate was kept from each infected leaf.

Table 1

Geographic coordinates, annual mean temperatures, and summary of statistics for Avr3a in six Phytophthora infestans populations

Population	Location	AMTa	No. of sequences	S b	No. of haplotypes	Haplotype diversity	Nucleotide diversity
Ningxia	Guyuan	7.0	19	25	8	0.719	0.007
Gansu	Tianshui	11.7	17	46	11	0.912	0.014
Guizhou	Anshun	14.7	17	46	11	0.882	0.021
Yunnan	Kunming	15.6	17	55	13	0.926	0.018
Fujian	Fuzhou	20.5	11	45	10	0.982	0.019
Guangxi	Nanning	22.6	15	43	10	0.857	0.023
Total			96	117	51	0.884	0.017

AMT = Annual mean temperature in collection sites.

S = Number of variable sites.

DNA extraction and Avr3a sequencing

Phytophthora infestans isolates were retrieved from long‐term storage and cultured on rye B agars at 18°C in the dark. Mycelia (~100 mg) were harvested 15 days after inoculation, transferred into sterile, 2‐mL centrifuge tubes, and lyophilized with a vacuum freeze dryer (Alpha1‐2, Christ, Germany). The lyophilized mycelia were ground to powder with a mixer mill (MM400, Retsch, Germany). Total DNA was extracted using a Plant gDNA Miniprep Kit (GD 2611, Biomiga, China) according to the manufacturer's instructions. The genomic DNA was suspended in 200 μL of ultrapure water and stored at − 20°C until use.

Genomic DNA was amplified with primers (F: 5′‐CCATGCGTCTGGCAATTATGCT‐3′, R: 5′‐CTGAAAACTAATATCCAGTGA‐3′, Armstrong et al., 2005). These primers were originally designed for Pex147 (Armstrong et al., 2005) but were successfully used to amplify Avr3a in a previous study (Cárdenas et al., 2011). PCRs were carried out in a 25 μl reaction volume using Gene Cycler™ (Bio‐Rad, Shanghai). Each reaction contained 1× PCR buffer, 100 μM dNTPs, 1 unit of Taq polymerase (TaKaRa Ex Taq®), 0.32 μM of primers and 20 ng of template DNA. PCR amplification was started with an initial denaturation step of 94°C for 4 min, followed by 35 cycles of 94°C for 40 s, 58°C for 1 min, 72°C for 1.5 min, and ended by an extension cycle of 72°C for 5 min. PCR products were isolated by electrophoresis and purified for single direction sequencing according to the manufacturer's instructions (QIAquick® Gel Extraction Kit), ligated into T5 zero cloning vector and transformed into Trans1‐T1 competent cells by heat‐shock at 42°C for 30 s (pEASY®‐T5 Zero Cloning Kit). Three colonies were picked from each transformation and incubated in liquid LB media at 37°C with shaking. One colony was randomly picked and sequenced by GenScript Biological Technology Co., Ltd. (GenScript, Nanjing, China) using an ABI3730 automated DNA sequencer (Applied Biosystems, USA).

Population genetic analyses and gene networks

Avr3a isoforms were deduced from nucleotide sequences, and multiple sequence alignments were performed using the MUSCLE algorithm (Edgar, 2004) implemented in MEGA5 (Tamura et al., 2011). Haplotypes were reconstructed by PHASRE algorithm implemented in DnaSP 5.10 (Librado & Rozas, 2009). The DnaSP 5.10 program was also used to estimate haplotype diversity, nucleotide diversity, and overall population differentiation (G ST) in Avr3a. Haplotype and nucleotide diversities were estimated for each of the six populations as well as the combined population by pooling the sequences for individual populations. SSR data of the isolates were taken from previous publications (Wu et al., 2016; Yang, Zhu, et al., 2016), and overall population differentiation in the SSR marker loci was estimated by fixation index (F ST) using POPGENE 1.32 (http://www.ualberta.ca/~fyeh/popgene_download.html). The standardized deviation of K ST was generated by bootstrapping with 100 replicates and used to compare G ST in the Avr3a gene and F ST in SSR marker loci by a t test (Gao, Zou, Xie, & Zhan, 2017). A median‐joining (MJ) network illustrating genealogical relationships among haplotypes was generated using Network 5.0 (Bandelt, Forster, & Rohl, 1999).

Intragenic recombination, natural selection, and phylogenetic relationship

Putative intragenic recombination events and parental sequences were identified with seven algorithms (RDP, GENECONV, Bootscan, MaxChi, Chimaera, SiScan, and 3Seq) using the RDP4 suite (Martin, Murrell, Golden, Khoosal, & Muhire, 2015). The likelihood of putative recombination and parental sequence detections was corrected by a Bonferroni procedure with a cut‐off of p < .01, and only detections supported by at least four of the seven algorithms were retained. Recombination events were confirmed and displayed by similarity plots implemented in SimPlot 3.5.1 (Lole et al., 1999), a window size of 20 nucleotides and a step size of two nucleotides. The recombinant sequences were excluded from subsequent analysis of phylogenetic relationships.

Selective neutrality in Avr3a was evaluated by Tajima's D (Tajima 1989) and Fu's F S (Fu, 1997) simultaneously using Arlequin 3.5 (Excoffier & Lischer, 2010). Recent selective sweep or purifying selection can cause an excess of rare polymorphisms to neutral expectation in sequences (Tajima 1989, Fu, 1997), leading to negative estimates of the two indexes in Avr3a. On the other hand, balancing or diversifying selection may generate positive estimates of the indexes as a result of deficiency of polymorphisms with low and high frequencies (Tajima 1989, Fu, 1997). The hypothesis of selective neutrality in Avr3a is retained if the estimated indexes did not differed significantly to the theoretical expectation of zeroes. Selective neutrality in Avr3a was also evaluated by the ratio of nonsynonymous (dN) to synonymous (dS) substitution (ω=dN/dS using four codon‐based algorithms implemented in HyPhy 2.10b (Pond et al. 2005) including FEL (fixed‐effects likelihood), IFEL (internal branches fixed‐effects likelihood), FUBAR (Fast Unconstrained Bayesian AppRoximation), and MEME (mixed effects model of evolution) (Yang 2007). Only sites simultaneously identified by FEL, IFEL, and MEME with p < .05 and >.95 posterior probability identified by FUBAR were considered to be under selection.

A Bayesian phylogenetic tree of Avr3a sequences was reconstructed using MrBayes 3.26 (Ronquist et al., 2012). The Bayesian analysis was carried out under the HKY substitution model, which was selected by MrModelTest (Nylander, 2008). Markov chain Monte Carlo chains (MCMC) were run for 106 generations and sampled every 100 generations. Chain stationary and convergence of running parameters were checked using TRACER 1.6, and node support was evaluated with posterior probability generated from 105 trees (burn‐in of 2500 trees).

Association between haplotype frequency and temperature

Thermal data for each collection site were downloaded from Weather Network (http://www.tianqi.com/). Annual mean temperature at collection sites was estimated using the air temperature across last 10 years as described previously (Yang, Zhu, et al., 2016). An association between isoform frequency (Avr3aEM) and the annual mean temperature at the collection sites was evaluated by Pearson's correlation (Lawrence & Lin, 1989).

RESULTS

Sequence variation in Avr3a

A total of 117 variable sites were detected in the 96 full nucleotide sequences, representing 11–19 from each of the six populations collected across China (Table 1, Figure 1). These variable sites formed 51 nucleotide haplotypes (Table 1) encoding 38 amino acid isoforms. The majority of sequence variations were generated by point mutations. Nearly 60% of the deduced amino acid isoforms differed only in one amino acid residue by a nonsynonymous mutation. Approximately 75% of the nonsynonymous mutations occurred in the effector domain and ~50% of those mutations were observed in the W and Y motifs. One isolate from Yunnan (YN2) had a mutation in the 397th nucleotide (C397T), generating an early termination stop codon causing a 12 C‐terminal amino acids truncation (Figure 2a). Four isolates from Fujian (FJ65 and FJ9), Guizhou (GZ116), and Guangxi (GZ116) populations had a single base deletion in the 10th, 93rd, 257th, and 273rd nucleotide, respectively (Figure 2a). Mutations in both start and stop codons were found in NX21113, a sequence from Ningxia (Figure 2a).

Figure 1

Frequency and spatial distribution of Avr3a haplotypes in the Phytophthora infestans isolates sampled from six geographical locations in China

Figure 2

Mechanisms of generating genetic variation in Avr3a gene of Phytophthora infestans: (a) DNA sequences showing point mutations, early termination, single base deletion, and defeated start and stop codons in Avr3a. Point mutations in start and stop codons are highlighted with bold letters, and the single base deletions are highlighted with bold dash lines; (b) intragenic recombinant generated by YN3, an Avr3aEM form, and YN8, an isoform in Clade 3

H‐1, H‐3, and H‐10 were the most common haplotypes, accounting for 16.7%, 30.2%, and 3.1% of the combined population, respectively (Figure 1). H‐1 and H‐3 were detected in all six populations, but H‐10 was only found in the pathogen populations from Southern China—in Fujian, Guangxi, Guizhou, and Yunnan (Figure 1). All other haplotypes were detected only once and were hence private to each population. The deduced amino acid isoform of H‐3 was identical to virulent Avr3aEM, while the amino acid sequence of H‐1 was one AA residue (R124G) different from virulent Avr3aEM. The deduced amino acid isoform of H‐10 was not found previously in the gene.

The haplotype diversity of nucleotide sequences in the six populations ranged from 0.719 to 0.982 with a grand mean of 0.884, when the 96 sequences from the six populations were combined. Nucleotide diversity of the populations ranged from 0.007 to 0.023 with a grand mean of 0.017 (Table 1). The Fujian population displayed the highest haplotype diversity, and the Guangxi population had the highest nucleotide diversity. The lowest haplotype and nucleotide diversities were found in the Ningxia population. The overall genetic differentiation across the six populations in haplotype (G ST) and eight SSR marker loci were 0.017 and 0.296, respectively. G ST was significantly smaller than F ST (p < .01) by a two‐tailed t test.

Haplotype network and geographical distribution

The 51 haplotypes were grouped into three clades by Bayesian inference (Figure 3a). Clade 1 was composed of all haplotypes translated to the Avr3aEM form, and Clade 2 was composed of haplotypes translated to isoforms similar to Avr3aEM with a single AA residue change. Clade 3 consisted of four haplotypes, all except one (H‐46) translated to the same isoform. Nucleotide sequences from Gansu, Guizhou, Guangxi, and Fujian populations were found in all three clades with more than half being grouped into Clade 1. Sequences from the Yunnan population were also found in all three clades, but more than half were grouped in Clade 2. Sequences from the Ningxia population were only found in Clade 1 and Clade 2 (Figure 3). Nucleotide haplotypes within a clade were connected to each other by less than two base substitutions, and haplotypes among the clades were connected to each other by two to four steps of base substitution through a common progenitor sequence (m/1), which was not detected in the samples (Figure 1b). The combined frequency of Avr3aEM (Clade 1) and its AA residue change alike (Clade 2) in the six populations ranged from 86.7% to 100% with an average of 91.4% (Figure 3b), while that of Clade 3 ranged from 0.0% to 13.3% with an average of 8.6% (Figure 3b). The frequency of the Avr3aEM form was negatively correlated with the annual mean temperature of the collection sites (r = −.826, DF = 4, p = .041, Figure 4).

Figure 3

Phylogenetic relationship and haplotype network in the Avr3a sequences generated from 96 Phytophthora infestans isolates collected from six geographical locations in China; (a) phylogenetic relationships among the 51 Avr3a nucleotide haplotypes of P. infestans reconstructed by a Bayesian inference. Posterior probability of topology was generated by 1,000,000 bootstrapping and is shown on the three main clusters. The most common (H‐1, H‐3, and H‐10) and the recombinant (H46) sequences were highlighted with bold letters; (b) haplotype network of 51 Avr3a nucleotide sequences and circle sizes of the circles represent haplotype frequencies in populations

Figure 4

Correlation between annual mean temperature in collection sites and the frequency of Avr3aEM

Intragenic recombination and neutrality test in Avr3a

Intragenic recombination was detected in one (YN17) of the sequences from the Yunnan population by four of the seven models with high confidence (GENECONV, p = 1.5−4; MAXCHI, p = 7.22−9; CHIMAERA, p = 5.66−9; and 3SEQ, p = 1.08−13). The recombinant was probably generated through interaction between YN3, an Avr3aEM form sequence, and YN8, a new isoform in Clade 3 (Figure 2b). As with the recombinant sequence, both parental ones were detected in the Yunnan population. Tajima's D and Fu's F S tests in the combined population from the six locations were −2.214 (p = .003) and −24.967 (p = .000), respectively. Both of these values were significantly lower than that expected under selective neutrality (zero), suggesting the existence of purifying selection acting on Avr3a. The ratio of nonsynonymous to synonymous substitution (dN/dS) was less than one in all variable sites but significance was found only in the 61st (p = .04) and 126th (p = .03) codons.

DISCUSSION

Genetic diversity of Avr3a gene was found to be high in the current study, and nucleotide diversity (Table 1) was positively correlated (r = .72, p = .10, data not shown) to SSR diversity (Qin et al., 2016). More than 50 nucleotide haplotypes were detected in 96 sequences, and the combined nucleotide diversity was 0.017 (Table 1). Effector genes are critical for the invasion, colonization, and reproduction of pathogens on hosts, and the finding of high genetic diversity in the effector gene Avr3a is consistent with the evolutionary hypothesis postulating that genes involving in antagonistic host‐pathogen coevolution have higher evolutionary rates compared to other genes as documented both in P. infestans (Cárdenas et al., 2011; de Vries et al., 2017) and many other species (Allen et al., 2004, 2008; Raffaele, Win, Cano, & Kamoun, 2010). This high evolutionary rate may well contribute to the rapid “breakdown” of many host resistances mediated by major genes (Cooke et al., 2012; Pilet et al., 2005). In previous studies, only two haplotypes formed by three SNPs were detected in 55 isolates of P. infestans (Armstrong et al., 2005) and six haplotypes generated by 12 SNPs were found in a different set of 88 sequences (Cárdenas et al., 2011). High haplotype diversity in the current study is unlikely due to artifacts generated by PCR amplification. The error rate for the standard Tag PCR‐based sequencing is 2.0 × 10−5 (Lundberg et al., 1991), and the error rate for the Ex Tag PCR‐based sequencing is approximately 4.5 times lower than that generated by standard Tag PCR based (http://www.clontech.com/MR/Products/PCR/High_Yield_PCR/Ex_Taq_DNA_Polymerase). Avr3a is 1,448 bp in length, and 96 haplotypes were sequenced in the study. Using this information, we estimated that the potential number of artifact sequences in our study was less than one (2 × 10−5 × 1,448 × 96/4.5 = 0.62). Furthermore, when isolates with “minor” haplotypes were resequenced with five independent clones, all but one “minor” haplotypes were recaptured in the new amplification.

Two factors may contribute to the differences in diversity found between the current and previous studies. Genetic variation, and therefore evolutionary potential, of populations can be underestimated when samples are over‐represented by identical genotypes (clones) generated by asexual reproduction of parents. This sample error is usually corrected using only one representative of the same genotype in the estimation of genetic variation (Zhan, Pettway, & McDonald, 2003). Unlike previous studies, isolates in the current study were prescreened molecularly and morphologically and only isolates with distinct genotypes were selected for sequence analysis (Wu et al., 2016). China also plants ~30% of the global potato acreage (http://faostat.fao.org/), thereby possibly hosting the largest P. infestans population in the world. Moreover, epidemiological analysis also suggests that the overall Chinese population of P. infestans has continued to expand in recent decades (unpublished data). Large population size coupled with its continuing expansion creates optimal conditions for the generation and maintenance of genetic variation in the Chinese pathogen population. Interestingly, no shared AA sequences were found between the current and previous studies except the Avr3aEM virulent form (Armstrong et al., 2005; Cárdenas et al., 2011). In 2005, the frequency of Avr3aKI accounted for ~15% of the population (55 isolates, Armstrong et al., 2005) but its frequency reduced to ~7% in 2011 (88 isolates, Cárdenas et al., 2011). In this study, no Avr3aKI sequences were detected in 96 isolates, suggesting either this isoform has gradually disappeared from the pathogen populations or is patchily distributed in space. Instead, Avr3aKI was replaced by a new isoforms deduced from H‐10, H‐15, H‐20, H‐23, and H‐27 with a sequence identical to Pex‐147‐2 (Armstrong et al., 2005).

Sequence analysis indicates that, as for many other effector genes (Cooke et al., 2012; Raffaele & Kamoun, 2012), various types of mutations are responsible for the level of genetic variation detected in Avr3a including point mutation, early termination, frameshift, and defeated start and stop codons (Figure 2a). Although nucleotide substitution was the main mechanism, other types of mutations account for ~15% of haplotype variation observed. Avr3a proteins secreted from P. infestans during infection trigger defense responses in hosts carrying corresponding resistance genes (Lo Presti et al., 2015) and are also important for the pathogenicity of plant pathogens by increasing host susceptibility or suppressing host immunity (Bos et al., 2010; Haas et al., 2009; Vetukuri et al., 2011). Although increasing the invasive opportunity of P. infestans due to the prevention of its recognition by potato hosts carrying with R3a gene, high mutation rates in Avr3a particular through pseudogenization also abruptly disrupt the normal biochemical functions of the gene. This result is consistent with previous hypothesis that gene duplication may occur frequently in Avr3a of P. infestans (Armstrong et al., 2005). Alternatively, it indicates secondary or tertiary mutations may occur in the genome to compensate for the fitness penalty association with loss of the Avr3a functions required by the pathogen. These compensatory mutations have been widely documented (Jochumsen et al., 2016; Sato et al., 2016). Continuing genome expansion (Raffaele & Kamoun, 2012) may provide a unique environment for the development of compensatory mutations leading to functional redundancy or epitasis (Birch et al., 2008; Oh et al., 2009) in the effector genes of P. infestans.

In addition to mutation, our results confirm that intragenic recombination also contributes to the generation of high genetic variation in Avr3a. Unlike intergenic recombination that generates new genotypes through random reshuffling of unaltered existing genes, intragenic recombination occurs among varying nucleotide sequences within a single‐gene locus, generating new allelic variation with potentially novel properties (Watt, 1972). When one considers that chromosomes are essentially simple strings of nucleotides separated into many genes of varying biological and ecological functions, recombination can potentially occur at any point regardless of whether it is within a particular gene or among genes. Hence, it is logical to expect that the number of intragenic recombination events may be similar to the number of intergenic ones although an unknown proportion of intragenic recombination may be deleterious or lethal. Despite this possibility, its role in the evolution of species, largely underestimated until recent advances in molecular and computation technologies, has been shown to be of great importance occurring in almost all species and kingdoms ranging from virus (He & Ding, 2012; He et al., 2010; Phan, Okitsu, Maneekarn, & Ushijima, 2007) to protists (Ferreira & Briones, 2012), fungi (Stergiopoulos et al., 2014), plants (Kelly et al., 2010; Ortega, Bošković, Sargent, & Tobutt, 2006; Städler & Delph, 2002), and animals (Godinho, Mendonça, Crespo, & Ferrand, 2006; Marthaler et al., 2014). In the current study, we have two lines of evidence to support the hypothesis that intragenic recombination also contributes to the high genetic variation and evolutionary potential of Avr3a. A chimeric sequence and its putative parents were detected in the isolates from the same population. This was supported by several algorithms with a high level of confidence (Figure 2b). In addition, a reticulation pattern in the haplotype network of some sequences (Figure 3b) also points to the possibility of intragenic recombination in Avr3a. Sequences chimeras could also be generated by joining two biological sequences together during PCR amplification, and this type of sequence chimeras usually disappear after one repeat of PCR re‐amplification (Smyth et al., 2010). However, the chimeric feature of YN17 sequence persists when the isolate was reamplified and resequenced, suggesting it is unlikely to be the artifact of PCR amplification.

Although large numbers of nucleotide haplotypes (51) and amino acid isoforms (38) were detected in the 96 sequences, only a single amino acid form (AVR3aEM) representing two nucleotide haplotypes (H‐1 and H‐3) dominated the P. infestans populations. All other isoforms including the new isoforms deduced from H‐10, H‐15, H‐20, H‐23, and H‐27 had a frequency of less than 5%. These results suggest that Avr3a gene might have experienced a selective sweep during its evolutionary history (Armstrong et al., 2005). Selective sweeps, usually associated with purifying selection, occur when environments select for the same traits/genotypes among populations inhabiting different ecological niches, thereby reducing genetic differentiation in related traits/genotypes among geographic populations (Zhan et al. 2004; Karasov et al., 2014). The hypothesis of a selective sweep occurring in the effector gene was supported by a comparative population genetic analysis of Avr3a and SSR marker loci which showed that genetic differentiation in Avr3a (assessed through the G ST statistic) was significantly lower than genetic differentiation in the SSR marker loci (assessed through the F ST statistic, Zhan et al., 2005). Furthermore, neutrality tests showed both Tajima's D and Fu's FS to be significantly less than zero. Lower genetic differentiation in Avr3a than the SSR marker loci was also detected previously by phenotypic (virulence) analysis of P. infestans isolates on potato differentials (Wu et al., 2016). It is likely that H‐1 (Avr3aEM) survived the selective sweep and other haplotypes (isoforms) were derived from the progenitor through mutation and intragenic recombination. The observation of selective sweep indicates Avr3a is an important gene contributing to survival, reproduction, and transmission of P. infestans, consisting with functional analysis of the gene (Bos et al., 2010). Indeed, six isolates with nonfunctional Avr3a sequences (Figure 2a) produced ~60% less disease on a susceptible cultivar (data not shown).

Host resistance is expected to be the primary force responsible for the observed pattern of natural selection in Avr3a. It is likely that potato cultivars with R3a gene have been intentionally or unintentionally used across many regions of the world, selecting for Avr3a isoforms able to circumvent the resistance. Indeed, the resistance gene has been detected in many potato cultivars recently (Nowicki, Foolad, Nowakowska, & Kozik, 2012; Zhu, Li, Vossen, Visser, & Jacobsen, 2012). In addition, environmental factors such as temperature may interact with host resistance influencing the fitness of Avr3a in P. infestans. Temperature is one of the most important environmental parameters crucially impacting many aspects of host–pathogen interactions including host susceptibility (Menna, Nguyen, Guttman, & Desveaux, 2015), pathogen density (Mikkelsen, Jørgensen, & Lyngkjær, 2015), and effector gene expression (Banta et al., 1998). The hypothesis of temperature‐mediated Avr3a fitness is supported by significant associations of isoform frequency with the mean annual temperature in the sample sites (Figure 4). Like other studies (Stefansson, Willi, Croll, & McDonald, 2014; Yang, Zhu, et al., 2016), we correlated Avr3a frequency to annual mean temperature rather than mean temperature in growing season for two reasons: (i) Pathogen adaptation to temperature occurs both in the epidemic and saprotrophic phase; and (ii) the length of the growing season is difficult to determine for many pathogens with multiple hosts due to cross‐species transmission. Regardless, a similar pattern of association was found when the mean temperature of potato growing season was used but the level of significance was reduced (data not shown). This result suggests that the pathogen with the Avr3aEM form may adapt better to lower than higher temperatures. Avr3aEM cannot be recognized by R3a (Armstrong et al., 2005) and therefore is a virulent type of Avr3a. Indeed, most isolates with Avr3aEM sequences induced late blight symptoms on R3a differential plants in an assay (Wu et al., 2016). For example, Avr3aEM (Clade 1 and 2) in isolates from Ningxia and Gansu accounted 100% and 94% of the populations, respectively, which correlated well with the frequency of virulent phenotypes. Isolates from Fujian matched least between genetic and phenotypic data. In this population, Avr3aEM and virulent phenotype accounted for ~90% and ~70%, respectively. The lower phenotypic than genotypic frequency in Fujian and other three populations likely results from recessive nature of Avr3aEM or gain of function mutants (Armstrong et al., 2005; Bos et al., 2009). This result indicates that increases in air temperature during global warming may slow the appearance of Avr3a and other virulence factors in P. infestans (Wu et al., 2016).

Directional selection on plant pathogens, which can generate selective sweep under extreme conditions, is the main mechanism driving the rapid breakdown of host resistance observed in agriculture. Unlike intragenic recombination, the direction and magnitude of selection on plant pathogens can be manipulated through the change in agricultural practices. To slow down the evolution of plant pathogens and achieve durable resistance in hosts in modern agriculture, it is important to create field conditions favoring disruptive selection through the evolutionary deployment of host resistance such as a cultivar mixture (Zhan et al., 2014, 2015). In addition, elevated air temperature associated with anthropogenic activities may greatly threaten food security due to extended seasonality supporting more plant diseases or enhanced pathogen growth thereby driving more intense disease epidemics. The negative association between virulence frequency of plant pathogens and local temperature found in the current and previous studies (Wu et al., 2016) raises the possibility that increasing air temperature during global warming may negatively affect the generation and maintenance of novel virulence genes in P. infestans, extending the life span of resistance genes. However, further investigations by functional, competition, and experimental evolution analyses of effector genes under controlled temperatures are required to confirm the result.

CONFLICT OF INTEREST

None declared.

DATA ARCHIVING STATEMENT

Haplotypes generated in this study are deposited in GenBank under accession numbers: MH043150 (FJ‐F002, hap‐1), MH043151 (FJ‐F041, hap‐3), MH043152 (FJ‐F050, hap‐4), MH043153 (FJ‐F052, hap‐5), MH043154 (FJ‐F057, hap‐6), MH043155 (FJ‐F064, hap‐7), MH043156 (FJ‐F103, hap‐9), MH043157 (FJ‐F105, hap‐10), MH043158 (GZ‐GA005, hap‐11), MH043159 (GZ‐GA007, hap‐12), MH043160 (GZ‐GA009, hap‐13), MH043161 (GZ‐GA014, hap‐14), MH043162 (GZ‐GA029, hap‐15), MH043163 (GZ‐GA038, 16), MH043164 (GZ‐GA041, hap‐17), MH043165 (GN‐GN003, hap‐20), MH043166 (GN‐GN008, hap‐21), MH043167 (GN‐GN010, hap‐22), MH043168 (GN‐GN013, hap‐23), MH043169 (GN‐GN016, 24), MH043170 (GN‐GN033, hap‐25), MH043171 (GN‐GN057, hap‐26), MH043172 (GS‐Pd11201, hap‐27), MH043173 (GS‐Pd11220, hap‐28), MH043174 (GS‐Pd11239, hap‐29), MH043175 (GS‐Pd11248, hap‐30), MH043176 (GS‐Pd11251, hap‐31), MH043177 (GS‐Pd11322, hap‐32), MH043178 GS (Pd11330, hap‐33), MH043179 (GS‐Pd13213, hap‐34), MH043180 (GS‐Pd13220, hap‐35), MH043181 (NX‐Pd21203, hap‐37), MH043182 (NX‐Pd213100, hap‐38), MH043183 (NX‐Pd213183, hap‐39), MH043184 (NX‐Pd21366, hap‐40), MH043185 (NX‐Pd21418, hap‐41), MH043186 (YN‐YN003, hap‐43), MH043187 (YN‐YN009, hap‐44), MH043188 (YN‐YN015, hap‐45), MH043189 (YN‐YN017, hap‐46), MH043190 (YN‐YN021, hap‐47), MH043191 (YN‐YN038, hap‐48), MH043192 (YN‐YN055, hap‐49), MH043193 (YN‐YN068, hap‐50), MH043194 (YN‐YN071, hap‐51).