Insights into the origin of the invasive populations of Trioza erytreae in Europe using microsatellite markers and mtDNA barcoding approaches

Genome-wide characterisation of SSRs

We identified and mapped a total of 428,342 microsatellites across the 47,828 scaffolds of the unpublished genome sequence draft of T. erytreae using the GMATA software³⁵. The SSRs frequency was estimated at 765.6 SSRs/Mb, which means 1 SSR for every 1.09 Kb. In silico identified SSRs were distributed among ten types of in tandem repeated motifs (from di- to deca-nucleotides). Analysis of SSR distribution revealed that the di-nucleotide motifs (340,227) were the most abundant SSRs, with a frequency of 79.4%. Both tetra- (20,902) and tri- (61,839) nucleotide repeats comprised about 5–15% (Fig. 1A; Supplementary Data 1). The remaining motifs, from hepta- to deca-nucleotides, comprised less than 1.5% of total SSRs identified in this study (Fig. 1A). Considering the unknown orientation of DNA strands in the Tery6 draft genome sequence of T. erytreae, a further SSRs characterization was carried out grouping the repeat motifs into pairs of complementary sequences. According to this, GA/TC (36.6%) and CT/AG (31.9%) are the most frequent motif pairs, with a total frequency of 68.5% (Fig. 1B). Grouped motif pairs GC/GC (0.05%) and CG/CG (~ 0.02%) were the least abundant di-nucleotide motifs. In decrease order, the most abundant tri-nucleotide motif pairs were ATT/AAT, ATA/TAT, ACA/TGT, TAA/TTA, AAC/GTT, TTG/CAA, and AAG/CTT, which encompassed 9.8% of all identified grouped motif pairs. Occurrence frequency of the remaining grouped motifs, including the rest of tri- and those from tetra- to deca-nucleotides (552 all together), was less than 11% of all motif pairs (Fig. 1B). Our data analysis reveals that SSR markers of 10 bp were most frequent, accounting for about 10% all SSR markers identified in this study. The overall trend of SSR length distribution in the T. erytreae genome is that the frequency of occurrence of SSRs gradually decreases as their length increases (Fig. 1C).

SSR markers development for T. erytreae

Fifteen SSRs chosen from those repeated motifs identified in silico in this study (Table 1) were used as potential markers to investigate the genetic diversity, structure and phylogeography of T. erytreae individuals from populations in mainland Europe and the archipelagos of the Macaronesia. Scaffolds Tery6_s00034 (274,710 bp), Tery6_s02825 (48,689 bp) and Tery6_s07841 (26739 bp) were randomly selected based on their sequence length (long, medium, and short scaffolds, respectively). SSRs were selected on the base of their type of repeat motif (di, tri-, tetra- and penta-nucleotides), nucleotide composition and length (number of in tandem repeated motifs) (Table 1; Supplementary Data 2). For the scaffold Tery6_s00034, 11 SSR loci were chosen from the total of 106 SSRs identified in silico, three for Tery6_s07841 and one for Tery6_s02825. Selected scaffolds were further investigated to know whether SSR loci mapped into coding or non-coding regions (inter-genic or intron sequences). Although gene annotation of the T. erytreae genome draft is not yet completed, it was possible to get this information for most of the selected SSR loci (data not shown). The scaffolds Tery6_s00034, Tery06, − 07, − 13 and − 14 were found in inter-genic regions, while Tery08, − 12 and − 15 were mapped into introns. For Tery05, − 9, − 10 and − 11 was not possible to establish whether they targeted coding or non-coding regions. SSR loci Tery01, − 02 and − 03 were found in intron regions in the scaffold Tery6_s07841, and SSR locus Tery04 in an inter-genic region in the sequence corresponding to the scaffold Tery6_s02825. For amplification of SSR loci, specific PCR primers were designed on the sequence flanking the in tandem repeated motifs. Blast of the different amplicons against the T. erytreae draft genome sequence showed that PCR primers would result in the specific amplification of their specific SSR locus. Experimental validation of PCR primers was carried out on a testing panel of individuals collected in different locations in the Canary Islands and South Africa. Primers pairs for SSR loci Tery04, − 05, − 06, − 08, − 09, − 10, − 11, − 12, − 13 and − 15 yielded DNA fragments of the expected size and were chosen for carry on further population genetic analysis. These loci contain eight di-nucleotides (AC, AG, GA, CA, GT, TC, TA and TG), one tri-nucleotide (TGA), and three tetra-nucleotides (CATA, CTAC and TACC), which arranged in microsatellites of different length (from 5 to 30 in tandem repeated motifs) (Table 1). Five SSR loci (Tery01, − 02, − 03, − 07 and − 14) were not amplified efficiently and the corresponding primer pairs were discarded for further analysis.

The individuals of T. erytreae collected in different geographical locations in the west coast of mainland Spain and Portugal, the Canary Islands and Madeira, as well as in South Africa and Kenya (Table 2), were analysed using the 10 selected SSR markers designed in this study. The scored allelic data for each SSR marker is summarised in the Table 3. The analysis showed that all SSR markers were polymorphic. Seventy alleles were detected over the ten selected SSR loci, and the average number of alleles per locus (N_a) was seven. SSR markers Tery08 and Tery11 had the highest number of alleles (12 and 20 alleles respectively), whereas Tery13 had the lowest (only two alleles). The expected (H_e) and observed (H_o) heterozygosity per locus in the entire population ranged from 0.20 to 0.77 and from 0.03 to 0.84, respectively. SSRs Tery11 and Tery08 displayed the highest diversity (H_e of 0.77 and 0.72, respectively), and Tery09 and Tery13 (H_e of 0.20 and 0.22, respectively) were the least informative markers. Most of the SSR markers used in this work showed H_e values higher than 0.5, apart from Tery05, − 09 and − 13 (with values of 0.39, 0.20 and 0.22, respectively). With the only exception of Tery04 and Tery15, for most of the analysed SSRs H_e was higher than H_o. It can be also observed that the whole population displayed a deficit of average H_o (0.31) compared with the H_e value (0.51) under Hardy–Weinberg equilibrium. This observation agrees with the positive value of the Wright’s fixation index (F_w) estimated for all analysed SSR markers over the whole population (F_w = 0.41). The SSR markers Tery12 and Tery13 showed F_w values close to 1.0 (0.81 and 0.85, respectively), suggesting that their alleles were considerably fixed in the population.

Population structure based on T. erytreae SSR data

To assess the differentiation and genetic diversity among the local populations of T. erytreae sampled in newly invaded areas from Spain and Portugal, including Madeira and the Canary Islands, and those from the previous invaded areas in Africa (South Africa and Kenya), we used a Bayesian clustering method to analyse the SSR multi-locus genotyping data. The STRUCTURE analysis according to the method of ΔK³⁶ showed that the overall genetic profile of all the individuals sampled could be described with two or three different hypothetically original populations corresponding to the highest ΔK values (Fig. 2). It means that the most likely values of genetic clusters (K) are 2 or 3. Nevertheless, Pritchard’s method³⁷ showed a posterior probability of data at K = 7 (Fig. 2). The estimated likelihood distribution increased from K = 1 to K = 7, and then started to decrease. This implied that seven was the smallest value of K, which was the most likely number of inferred populations in our data set. Interestingly, the value of K at which the likelihood distribution reached its maximum coincided with a further peak value of the ΔK statistic at K = 7, suggesting a more complex hierarchical structure of the T. erytreae populations (Fig. 2). In consequence, we plotted the clustering results for K = 2, K = 3 and K = 7 (Fig. 3). Furthermore, we considered an initial structure of two populations (K = 2) as was suggested by the method of ΔK³⁶ whereby most of the analysed individuals were classified with high probability (Q > 0.90) in two clusters (Fig. 3). Cluster 1 (in green) was exclusively formed by individuals from newly invaded areas in Spain and Portugal, including those from the archipelagos of Madeira and the Canary Islands. On the other hand, Cluster 2 (in beige) was mainly comprised of individuals from Africa, but also included individuals from Camacha (Madeira). The exception to this pattern involved three locations in Madeira (Quebradas, Camacha and Moreno), Pretoria (South Africa), and Homa Bay (Kenya), where almost all individuals consistently had significant membership in both clusters. Looking at K = 3 plot, the Bayesian clustering analysis resolved Cluster 1 into two by reassigning some individuals to Cluster 3 (in purple). Almost of all individuals from Moreno, Poiso, and Farrobo (in Madeira and Porto Santo, respectively) were entirely reassigned to Cluster 3 along with several individuals from the Canary Islands and Galicia (Spain). In addition, individuals from Vairão (Porto) and São Vicente de Pereira Jusã (Aveiro) (both in the northwest coast of Portugal) were also assigned to Cluster 3, while those individuals sampled from southern locations up to Sobreda (Setúbal) were assigned to Cluster 1. The exceptions to this pattern were the individuals from Ribamar (Ericeira), which were assigned to Cluster 3. Most notably, samples from Kenya were genetically different from those of South Africa and grouped in Cluster 1. At K = 7 the population structure scenario was more hierarchical, but 73% of all individuals (108 out from 147) could be assigned to one of the seven clusters with more than 90% probability (Q > 0.9). The assignment of half of the remaining individuals (21 out of 39) could be done with more than 70% probability (Q > 0.7). Among the different groups, Cluster 1 (in green) and 2 (in beige) are restricted to the populations of South Africa and Kenya, respectively, with almost no presence of individuals from any of the newly invaded areas. Clusters 3 (in purple) and 4 (in pink) are mostly exclusive to the individuals from Madeira and Portugal mainland, although with some membership in the Canary Islands and Galicia. Cluster 5 (in light blue) and Cluster 6 (in orange) are represented by individuals from Madeira, the Canary Islands and Galicia, while the individuals from Camacha (Madeira) –the only ones that were collected from Casimiroa edulis La Llave & Lex. (Rutacea: Toddalioideae)—form exclusively Cluster 7 (in dark blue). Remarkably, Q fractions corresponding to Cluster 7 are present in the individuals from Nelspruit, Tzaneen, and some in Pretoria.

Genetic diversity analysis using T. erytreae SSR allelic data

The genetic diversity of T. erytreae populations was also assessed by means of a distance-based clustering method. The scored SSR allelic data obtained from the ten SSR loci developed in this study were used to calculate a genetic dissimilarity matrix and to compute a Neighbor Joining (NJ) tree. A preliminary dendogram constructed using only the African populations of T. erytreae showed that the individuals from South Africa grouped together into a single cluster clearly separated from the Kenyan population. The robustness of the tree clustering was supported by the high bootstrap values obtained for nearly all branches (Fig. 4). To confirm the results obtained from the structure analysis a NJ tree under topological constraints was inferred using as initial tree the population structure of individuals from all the sampled areas with Q > 0.7. The remaining individuals were positioned (constraint) on that previous topology. Inspection of the constrained tree topology revealed seven clusters that were in congruence with the structural population at K = 7 suggested by the STRUCTURE analysis (Fig. 5). It is noteworthy that Cluster 7 emerged as a paraphyletic group in the base of African Cluster 2. The cluster assignments of individuals with low membership coefficients (Q < 0.7) performed well in our distance-based clustering analysis.

Phylogenetic analysis using mtDNA-based barcoding

The maternal phylogenetic relationship among the T. erytreae individuals collected in this study was assessed using mitochondrial DNA (mtDNA) barcoding, the sequence analysis of nucleotide variations in the 5’ region of the Cytochrome C Oxidase I gene (COI)³⁸. Comparison of the nucleotide sequence of COI barcode fragments from this study with other T. erytreae GenBank accessions demonstrated that all sequences are highly conserved (Supplementary Data 3). With the only exception of the fragments amplified from Kenyan individuals, all remaining sequences showed absolute identity (100%) to the GenBank accession numbers that we previously deposited in the GenBank³⁹, corresponding to COI barcode sequences from T. erytreae individuals collected in the Canary Islands (MK285551-MK285553), Galicia (MK285548-MK285550), Madeira (MK285558) and South Africa (MK285554-MK285557, MK285559, MK285560) (Supplementary Data 4). Alignment of sequences amplified in this study from individuals collected in Homa Bay (Kenya) shared absolute identity with those extracted from the entire mitochondrion genome sequences from Eritrea and Uganda²⁷, and 97–98% identity with those COI barcode fragment sequences from individuals from other locations in Kenya³⁹.

The nucleotide sequences of the COI barcode fragments generated in this study (n = 39) and some previously deposited in the GenBank (n = 37) were used to analyse the maternal phylogenetic relationship of T. erytreae populations that have invaded Spain and Portugal, with those from South Africa and Kenya (Fig. 6). From all COI sequences used in this study, 38.2% were obtained from Spain and Portugal (including Madeira and the Canary Islands), 36.8% from South Africa, 19.7% from Kenya, and 5.3% from other African countries (Cameroon, Ethiopia, Tanzania and Uganda). In accordance with the high level of identity of their COI barcode nucleotide sequences, the NJ tree generated from these sequences showed that the individuals from Spain and Portugal, including those from Madeira and the Canary Islands, formed a monophyletic group with the individuals from South Africa (Pretoria, Nelspruit, and Tzaneen). Our phylogenetic analysis reveals a clear differentiation between this monophyletic group and the individuals from Homa Bay (Kenya), as well as from those individuals previously reported in other locations in Kenya³⁹. It was observed that the individuals from Spain and Portugal formed a paraphyletic group with those from Pretoria (Fig. 6), as the remaining South African individuals from Nelspruit and Tzaneen formed a separated clade. Furthermore, our analysis demonstrated the presence of two different T. erytreae lineages in Tzaneen, as most of their individuals formed a paraphyletic group with those from Nelspruit, while the remaining formed a clade with four individuals from West Acres (South Africa)³⁹. The few exceptions to this observation were three South African individuals, one from Pretoria (Pretoria-100), and two from West Acres (TeSA1 and TeSA7), which may correspond to migrants from Nelspruit or Tzaneen, and Pretoria, respectively. In a sister clade position to the South African clade, the GenBank accessions of Kenyan and Tanzania COI sequences³⁹ included in this phylogenetic analysis formed a monophyletic group. The COI barcode sequences from Homa Bay (Kenya) clustered separately as an outgroup in a different clade with the corresponding fragment extracted from the mitochondrion genome sequences from Ethiopia and Uganda (MT416551 and MT416549, respectively)²⁷, and Cameroon (MG989238)⁴⁰ present in GenBank.

Source: Ecology - nature.com