Studying genetic population structure to shed light on the demographic explosion of the rare species Barbitistes vicetinus (Orthoptera, Tettigoniidae)

Insect outbreaks usually involve important ecological and economic consequences for agriculture and forestry. The short-winged bush-cricket Barbitistes vicetinus Galvagni & Fontana, 1993 is a recently described species that was considered rare until ten years ago, when unexpected population outbreaks causing severe defoliations across forests and crops were observed in north-eastern Italy. A genetic approach was used to analyse the origin of outbreak populations. The analysis of two mitochondrial regions (Cytochrome Oxidase I and II and 12S rRNA-Control Region) of 130 samples from the two disjunct ranges (Euganean and Berici Hills) showed high values of haplotype diversity and revealed a high geographical structure among populations of the two ranges. The high genetic variability observed supports the native origin of this species. In addition, results suggest that unexpected outbreaks are not a consequence of a single or few pestiferous haplotypes but rather the source of outbreaks are local populations which have experienced an increase in each area. The recent outbreaks have probably appeared independently of the genetic haplotypes whereas environmental conditions could have affected the outbreak populations. These findings contribute to a growing understanding of the status and evolutionary history of the pest that would be useful for developing and implementing biological control strategies for example by maximizing efforts to locate native natural enemies.


Introduction
Herbivorous insect outbreaks usually involve important ecological and economic consequences for agriculture and forestry. In most cases, it concerns the accidental introduction of an exotic species into a new geographical area due to international trade and human movement [1][2][3]. In this scenario, the species is not under biological control by natural enemies, leading to a quick spread of the population. Instead, the insect outbreak of an indigenous iii. are either of the recent outbreaks caused by a single outbreak population followed by a spread pattern or have they derived from multiple local populations?
This information can help understand the factors related to the outbreak events of B. vicetinus and could provide insight to assist management of this outbreak pests.

Study area
The study area was located in north-eastern Italy (Veneto Region) and included the two disjunct ranges where outbreaks have occurred (Euganean Hills and the Berici Hills) and a third range where no outbreaks have been reported to date and the species is rare (Lessini Mountains). The three ranges are separated from each other by cultivated and inhabited alluvial plains devoid of woody vegetation (Fig 1B).
The Euganean Hills cover an elliptical area of approximately 180 km 2 and comprise approximately 100 hills of volcanic origin emerging from the alluvial sedimentary plain, with the highest elevation of approximately 600 m above sea level. They are characterised by numerous narrow and deep valleys, steep hills with different sun exposures and many microclimatic conditions that influence the vegetation. The annual mean temperature is approximately 12˚C and precipitation ranges from 700 to 900 mm [36] even if they are extremely variable due to the inner geomorphological variability.
The Berici Hills are an isolated plateau situated on the southern Vicenza plain, they cover almost 160 km 2 with a maximum elevation of approximately 450 m above sea level [37]. The climate of the area is characterised by an annual rainfall of 958 mm, and average daily temperature of -1˚C in January and 23˚C in July [38].
The Lessini Mountains are a triangular-shaped tableland, which occupies some 800 km 2 , at the transition between the Fore-Alps and the River Po Plain [39]. The mountain group reaches over 2000 m above sea level and is characterised by multiple valleys and long ridges that descend towards the plain.

Insect sampling
Adults of B. vicetinus were collected, using a sweep net, throughout the spring of three successive years (2015, 2016, and 2017). Specimens from the north, centre and south part of both outbreak areas (the Euganean and Berici Hills) as well as from two sites on the Lessini Mountains were collected from bush and tree canopies on marginal side of the forests (Table 1). To avoid sampling relatives, bush-cricket specimens were collected by sweeping an area at least 100 m 2 at each sampling site. After capturing, samples were immediately kept in 95% ethanol and taken to the laboratory where they were morphologically identified and stored in individual vials at −20˚C until DNA extraction.
[33] consisting of 2 min at 92˚C followed by 35 cycles with a denaturation step of 92˚C for 2 min, an annealing step of 52˚C for 30 sec, and extension step of 60˚C for 3 min, with a final extension of 72˚C for 7 min.
PCR products were checked through electrophoresis on 1.0% agarose gels stained with SYBR1 (Invitrogen), purified using Exonuclease and Antarctic Phosphatase (GE Healthcare) and sequenced at the BMR Genomics Service (Padua, Italy). Primers used for amplification were also used for sequencing.

Data analysis
DNA sequence chromatograms were quality checked, manually corrected when necessary, and aligned using MEGA X [42]. Low-quality regions found at the beginning and end of each sequence were trimmed, while low-quality sequences were not included in the analysis.
Haplotype and nucleotide diversity, as well as the pairwise genetic distances between populations, were calculated with Arlequin 3.5 [43] using a Kimura 2-parameters model. The presence of population differentiation was also tested by conducting exact tests of population differentiation with 100,000 steps in Markov chain, with 10,000 dememorization steps. To compare the partition of genetic variability among sampled populations, an analysis of molecular variance (AMOVA) [44] was performed using Arlequin.  A haplotype parsimony network with a probability cut-off of 95%, was reconstructed using the TCS 1.21 software [45] and PopART 1.7 [46] and used for depicting the geographical relationships among haplotypes. Ambiguous connections (loops) were resolved using approaches from coalescent theory based on three criteria: frequency, network location and geography [47,48].
The study of the past demographic history of the species was inferred using Arlequin 3.5 through the Tajima's D and Fu's Fs tests [49,50] and the mismatch distributions of the pairwise genetic differences [51]. Populations at demographic equilibrium or decreasing in size should provide significant positive D and Fs values with a multimodal distribution of pairwise differences, whereas populations that have undergone a sudden demographic expansion usually show significant negative D and Fs values with a unimodal distribution [51,52]. The sudden expansion model was tested through analysis of the sum of square deviations (SSD) and raggedness index (r) representing the modality of distribution, obtaining the corresponding P values with a parametric bootstrap approach (10,000 replicates).
For those populations that did not deviate from sudden expansion (p > 0.05) the time since expansion was calculated considering the relationship τ = 2 � u � t (where τ = age of expansion measured in units of mutational time, u = mutation rate per sequence and per generation, and t = number of generations since the expansion; [53]). The expansion time was then obtained by dividing the estimate of τ by the product of the sequence length in base pairs and the mutation rate per nucleotide (twice the per-lineage substitution rate; u = 2 μk) in percentage per year. Two substitution rates were used: the classical rate of 2.3% divergence per Myr [54] frequently used in Orthoptera [55][56][57][58][59] and the global mtDNA rate of 2.6% divergence per Myr proposed by Papadopoulou et al. [60].

Sampling
A total of 152 individuals were collected from the three natural ranges of B. vicetinus (Euganean Hills, Berici Hills and Lessini Mountains). Unfortunately, only 12 specimens were collected from the Lessini Mountains due to the rarity of the species in this range (Tab. 1).

Data analysis
The two fragments of the mitochondrial genes, COI-tRNALeu-COII and 12S-CR were successfully amplified and sequenced in 130 samples, with an average of 19.6 specimens for each of the three areas within the two main ranges (north, centre, and south of both the Euganean and Berici Hills). After quality assessment and trimming the sequences, high-quality sequences 799 bp long for the COI-tRNALeu-COII and 898 bp long for the 12S-CR fragment were obtained. Sequences of the COI-tRNALeu-COII fragment were translated with Transeq (EMBOSS: http://www.ebi.ac.uk/Tools/emboss/transeq/index.html) to exclude the presence of any nuclear mitochondrial pseudogenes. This bioinformatics tool translates nucleic acid sequences to their corresponding peptide sequences and identifies stop codons. Since these pseudogenes (NUMTs) are characterized by the accumulation of in-frame stop codons and indels [61] the absence of stop codons in the protein sequence can allow to exclude presence of NUMTs.
A total of 13 variable sites, including 8 parsimony informative sites and 5 singleton sites, were identified by the alignment of COI-tRNALeu-COII sequences while the 12S-CR alignment showed 6 parsimony informative sites and 7 singleton sites (13 variable sites in total). Sequences of each haplotype of COI-tRNALeu-COII and 12S-CR obtained in this study are available through GenBank accession numbers MW405351-MW405381.
Partition homogeneity test confirmed that COI-tRNALeu-COII and 12S-CR fragments bear a homogeneous signal (P = 0.28), allowing data to be pooled for further analyses. A final dataset including 130 concatenated sequences of 1697 bp was obtained. Further analyses were conducted considering the combined dataset.

Genetic variability and population structure
The diversity indexes for the concatenated dataset ranged between 0.28 and 0.90 for the haplotype diversity (H) and between 0.02 and 0.24 for the nucleotide diversity (π). Among populations from both the main outbreak areas, those from the north and the centre part of the Berici Hills showed the highest H and π values (H = 0.92, π = 0.13 and H = 0.80, π = 0.10 respectively) while populations from the south part of the Berici Hills showed the lowest variability. The distribution of haplotypes among all populations analysed and other summary statistics are shown in Table 1.
The presence of population differentiation was confirmed by the tests of population differentiation (P>0.001). When B. vicetinus genotypes were grouped based on the six main collection sites (north, centre, and south of both Euganean and Berici Hills), the locus-by-locus AMOVA revealed a significant geographic structure between populations of the two main disjunct ranges (P = 0.0014) (Euganean vs Berici Hills). To avoid bias as a result of the small sampling size of Lessini Mountains, populations from this area were not included in this analysis. The analysis showed that 40% of the variation was explained by differences among groups (Euganean and Berici Hills), whereas about 48% of genetic variation was explained within populations (Table 2).

TCS Network
TCS Network of the combined dataset revealed the presence of 32 haplotypes of which 15 were exclusive to samples from the Euganean Hills, 13 were exclusive to the Berici Hills and 2 regarded samples only from the Lessini Mountains (Fig 1). Only 2 haplotypes were shared by samples from different ranges: A1 included samples from the Euganean and Berici Hills, and J1 samples from both the Berici Hills and Lessini Mountains. The network evinced a geographical separation among haplotypes, with samples from the Euganean Hills showing a star-like pattern and samples from the Berici Hills connected among them (Fig 1A).
Going deeper into the network's characteristics, the most common haplotype, A6, included 34 samples coming from all three parts of the Euganean Hills (north, centre and south). Ten rare haplotypes including samples exclusively from the Euganean Hills were connected, separated by only one mutational step, to A6 in a star-shape. Five of them (A9, A11, B6, G6, and H6) were represented by only one sequence, six (C6, A10, F6, A8) by two sequences, and one (E6) by 4 sequences. In addition, haplotype A1, which included samples from both the Euganean and Berici Hills, was also connected to A6 separated from it by only one mutational step.
The second most common haplotype, I1, was composed of 20 sequences:15 belonged to samples from the three geographical populations of Berici Hills (north, centre and south)  (Fig 1). The parsimony networks of the single markers (COI-tRNALeu-COII and 12S-CR) showed a similar pattern, with most haplotypes including exclusively samples from only one disjunct distribution range (Euganean Hills or Berici Hills or Lessini Mountains) (S2 and S3 Figs).

Past demographic events
Tajima's D and Fu's Fs tests were applied in populations from the Euganean and Berici Hills in order to check for past demographic events. The null hypothesis of neutrality was rejected in populations from the Euganean Hills (D = -2.05; P = 0.02 and Fs = -12.51; P < 0.001), suggesting a past population expansion after a period of low effective sample size (Table 3). Berici populations showed only significant negative Fu's Fs value (Fs = -6.62; P = 0.004), suggesting that B. vicetinus populations of this hilly area did not conform to the theory of neutral evolution ( Table 3).
The mismatch distribution plots of both ranges were smooth and unimodal, revealing that these populations were undergoing population expansion (  Table 3).
Estimation of the expansion time showed that populations from the Euganean Hills lineage started to expand about 10,700 years ago (with a 2.6% substitution rate) and 12,400 years ago (with a 2.3% substitution rate) ( Table 3). Populations from the Berici Hills probably expanded from about 24,000 to 27,600 years ago with 2.6% and 2.3% substitution rates, respectively (Table 3).

Discussion
This study presents the first population genetic analysis of the species B. vicetinus in its outbreak areas and reveals the presence of a high geographical structuring among populations of the two outbreak ranges analysed (Euganean and Berici Hills). Populations from both ranges showed high values of haplotype diversity, a typical characteristic of ancestral populations [63][64][65], supporting that B. vicetinus is a native species. Conversely, populations of invasive alien species are traditionally thought to have reduced genetic variation relative to their source populations because of genetic founder effects linked to small population size during the introduction and establishment phases of an invasion [66]. Furthermore, the genetic variability found clustered according to geographical ranges is in contrast with the possibility that B. vicetinus may be an invasive species accidentally introduced.
Even though outbreaks have been reported on both the Euganean and Berici Hills, different haplotypes have been found in these two distribution ranges. Ninety-four percent of haplotypes were exclusive to a single distribution range (i.e. either Euganean or Berici Hills) and only one haplotype, scarcely represented, was found in samples from both hilly areas. These results suggest that outbreaks are not a consequence of a single or few pestiferous haplotypes but rather that the source of outbreaks is due to local populations which experienced a demographic increase in each area. Thus, it seems that outbreaks have appeared independently from the genetic origin, as also found in some studies which indicate that outbreak events are often more affected by environmental conditions than by genetic characteristics of the local population [4,[67][68][69].
The geographical separation between populations of the Euganean and Berici Hills, observed in the network and confirmed by population differentiation test and AMOVA, indicates a limited gene flow among populations. In phytophagous insects, dispersal capacity, geographical or reproductive barriers, host plant, and habitat fragmentation are reported as the main drivers of genetic structure [20, 70,71]. The limited dispersal ability of this flightless species could have favoured the lack of gene flow among its distribution ranges. B. vicetinus has been reported to show a low dispersal ability [9], as well as other ground-dispersing species that move only relatively short distances, such as 100-200m during their whole life (e.g. Pholidoptera griseoaptera, [72]). In addition, the mostly lowland areas between the Euganean and Berici Hills, with the presence of agricultural fields and the absence of woody vegetation, might be hostile areas for bush-cricket survival and could have acted as a geographical barrier limiting the effective dispersal of the species. Furthermore, spatial configuration of habitats is another factor influencing genetic structure, that affects mainly species with limited dispersal ability [73]. In B. vicetinus outbreak areas, it has been observed that habitat loss and the presence of patchy areas of the non-host alien tree Robinia pseudoacacia play an important role in reducing population density and dispersion [9]. These factors, coupled to the low mobility of the pest, could have a synergic action that might explain the high level of differentiation observed among populations across its distribution ranges.
High haplotypic diversity, low levels of sequence divergence (nucleotide diversity) and a star-like phylogeny were observed in populations of B. vicetinus from the Euganean Hills and to a lesser extent from the Berici Hills (Fig 1; Table 1). This pattern is consistent with what is expected for populations that have experienced past demographic expansions [52,74]. These results were highlighted by the neutrality tests and the unimodal mismatch distribution in both hilly areas supporting that, besides the current and ongoing outbreaks, populations from the Euganean and Berici Hills underwent an expansion after a period of low effective sample size (S4 Fig). The bush-cricket population of the Euganean Hills experienced a postglacial expansion starting approximately 10,700-12,400 years ago after the Last Glacial Maximum (LGM; 21,000 years ago) whereas the expansion process in populations on the Berici Hills could have occurred at the end of the LGM. During this glacial period, climate effects were minimal in the Euganean Hills, with thermophilic vegetation serving as a refuge for several species [75,76]. Accordingly, B. vicetinus populations could have survived during the climatic oscillations, and exploited these hilly areas with potentially suitable environmental conditions throughout the late Pleistocene. Once climatic conditions were favourable, at the end of the last ice age, populations might have experienced an expansion, shaping the present genetic structure of the species. The slow movement of the species and lack of host plants in the lowlands during these periods probably prevented dispersal of the species between the two hilly areas, favouring differentiation of the mitochondrial haplotypes. Thus, both historical (e.g., expansion after LGM) and recent changes (e.g. habitat loss) have contributed to determining the genetic structure and diversity of the bush-cricket in its outbreak range, resulting in genetically structured populations.
Further studies increasing the sampling size, mostly in areas where the species is rare (e.g. Lessini Mountains), combined with analysis of other markers could help obtain a better picture of B. vicetinus population structure.
The current findings contribute to a growing knowledge of the status and evolutionary history of the pest. Here, by removing doubt about the origin of the B. vicetinus, we have achieved an important step in understanding this native species. Shedding light on the origin of a species is also important for its biological control. Given that natural enemies and their host tend to coevolve, biological control programmes often rely on the use of parasitoids present in the native areas. Failure to identify the correct origin of a pest may lead to the use of unsuitable species as biocontrol agents with negative effects on control programmes [20,77]. Maximizing efforts to locate native natural enemies (native parasitoids), as the egg parasitoid reported in Ortis et al.
[13] will be useful for control programmes. Finally, although B. vicetinus is currently in outbreak status, control programme strategies should consider the intrinsic vulnerability of an endemic species. Its low dispersal ability and small and fragmented distribution range could lead this native species to quickly shift from outbreaking to endangered.