Research Article Journal of Orthoptera Research 2019, 28(1): 11-19 Epigenetic and genetic variation between two behaviorally isolated species of Neoconocephalus (Orthoptera: Tettigonioidea) GIDEON Ney!, JOHANNES SCHUL! 1 Division of Biological Sciences, University of Missouri, 207 Tucker Hall, Columbia, Missouri, USA. Corresponding author: Gideon Ney (gideon.ney@gmail.com) Academic editor: Alina Avanesyan | Received 23 August 2018 | Accepted 31 January 2019 | Published 17 May 2019 http://zoobank.org/AD0139D9-3350-4863-95B5-E19ED5D3F6B2 Citation: Ney G, Schul J (2019) Epigenetic and genetic variation between two behaviorally isolated species of Neoconocephalus (Orthoptera: Tettigonioidea). Journal of Orthoptera Research 28(1): 11-19. https://doi.org/10.3897/jor.28.28888 Abstract Epigenetic variation allows for rapid changes in phenotypes without alterations to nucleotide sequences. These epigenetic signatures may di- verge over time among isolated populations. Epigenetic incompatibility following secondary contact between these populations could result in the evolution of reproductive isolating mechanisms. If epigenetic incompati- bility drove the evolution of species isolating mechanisms, we expect to see significant epigenetic differentiation between these species. Alternatively, epigenetic variation could be the result of predominantly environmental variables and not align along species boundaries. A methylation sensitive amplified fragment length polymorphism analysis was performed on in- dividuals of the closely related katydid species Neoconocephalus robustus and N. bivocatus. We observed significant variation in total methylation levels between species. However, genetic differentiation remained larger than epigenetic differentiation between species groups. We measured a sig- nificant correlation between the epigenetic and genetic distance between individuals. Epigenetic differentiation is therefore likely the result of an interaction between genetic and epigenetic loci and not a mechanism for species differentiation. We therefore did not find evidence to support our hypothesis of an epigenetically mediated mechanism for speciation be- tween N. robustus and N. bivocatus. Key words genetic differentiation, methylation, MSAP, Neoconocephalus, population genetics Introduction Epigenetic variation can lead to changes in phenotypic expres- sion without any change to the nucleic acid sequence (Berger et al. 2009). Unlike genetic variation, epigenetic changes can pro- duce reversible, heritable phenotypic changes within a lineage Jablonka and Lamb 1998). Epigenetically controlled phenotypic plasticity may therefore play an important role in rapid, adaptive changes (Led6on-Rettig et al. 2012, Burggren and Crews 2014). DNA methylation, the most commonly studied form of epige- netic modification, involves the addition of methyl groups, usu- ally to CpG dinucleotides that regulate gene expression (Boyes and Bird 1991). Methylation regulated phenotypic plasticity has been documented across multiple taxonomic groups (Braam and Davis 1990, Rapp and Wendel 2005, Varriale and Bernardi 2006, Elango et al. 2009, Herrera and Bazaga 2010). Until recently, epi- genetically dependent variation has been largely overlooked as an evolutionary mechanism contributing to speciation (Smith et al. 2016). Facultative phenotypes can be maintained in a population as adaptive alternatives to divergent environmental conditions. These alternative phenotypes may ultimately become fixed in a population by way of genetic assimilation (Pal and Miklos 1999, West-Eberhard 2005). Phenotypic variation can also evolve in response to epigenetic incompatibility (= hybrid incompatibility caused by epialleles; Bente and Scheid 2017, Blevins et al. 2017) between groups. Fol- lowing the epigenetic diversification of groups in isolation, inter- mediate forms may show reduced fitness following secondary con- tact Jablonka and Lamb 1999). Aphids moved to a new host plant and allowed to reproduce asexually quickly developed a preference for the new host and produced morphological characteristics simi- lar to conspecifics utilizing this host. Backcrossing to the parental line resulted in the production of non-viable offspring (Shaposh- nikov 1966). Heritable phenotypic variation in aphids has been correlated with changes in methylation (Field et al. 1989, Walsh et al. 2010). Methylation dependent incompatibility could be rein- forced by reproductive isolation (Pal and Miklds 1999). In many animals, behavioral isolation plays a significant role in maintaining species boundaries. Acoustic communication has been studied as a mechanism for reproductive isolation in both vertebrate and invertebrate groups (Gerhardt and Huber 2002). Neoconocephalus, a New World genus of katydids, possesses a di- verse range of habitat preferences (Walker et al. 1973, Frederick 2013) and call phenotypes (Schul and Patterson 2003, Beckers and Schul 2008, Deily and Schul 2009). Females utilize these calls in mate recognition and in phonotaxis, the directional movement to- wards a calling male (Greenfield 1990, Triblehorn and Schul 2009, Beckers and Schul 2010). The acoustic communication system of Neoconocephalus allows for reproductive isolation among the mul- tiple species that may be found in sympatry (Schul et al. 2014). The sibling species (Snyder et al. 2009) N. robustus (Scudder, 1862) and N. bivocatus Walker, Whitesell, & Alexander, 1973 dif- JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) 12 fer both in their habitat preferences and call types. Neoconocepha- lus robustus is a grassland generalist and utilizes grasslands with a wide range of flora and environmental conditions. Neoconoceph- alus bivocatus, on the other hand, is a prairie obligate preferring drier habitats composed primarily of tall prairie grasses (Walker et al. 1973). Both species co-occur in some areas, e.g. at the edges of prairie remnants or along creeks within prairies. This species pair is largely morphologically cryptic but maintains divergent call characters and call preferences (Deily and Schul 2004, 2006). While the two species can be distinguished using genetic mark- ers (e.g. mitochondrial haplotypes or AFLP markers (Snyder et al. 2009, Ney and Schul 2017)), genetic differentiation is weak and evidence for hybridization (e.g. intermediate calls) occur fre- quently. Population-level genetic variation was shown to be very low, with almost no geographically associated genetic variation between sites as distant as 448 km apart (Ney and Schul 2017). This suggests that these two species diverged very recently, most likely during the current glacial cycle (Frederick 2013). A likely scenario of this divergence event includes an ancestral population living in both habitat ranges. In each habitat, different epigenetic patterns would be expressed, which then became fixed within each subpopulation, resulting in decreased hybrid fitness. Genetic differences leading to the differences in male calls and fe- male preferences evolved later, in part driven by reproductive rein- forcement (Deily and Schul 2004). Thus we hypothesize here that epigenetic differentiation, driven by the different environmental conditions experienced in their core habitats, drove divergence of the species. The lack of substantial genetic differentiation between N. robustus and N. bivocatus could point to an underlying epige- netic incompatibility driving species differentiation. This hypoth- esis predicts that the epigenetic differentiation between the species would be larger and more stereotyped than the genetic differences. DNA methylation has been described in many insect groups including multiple orthopteran species (Sarkar et al. 1992, Robin- son et al. 2011). Methylation sensitive amplified fragment length polymorphism (MS-AFLP) analysis is one of the few techniques that allows for the quantification of genome-wide patterns of cy- tosine methylation without any previous knowledge of genome sequences (Xiong et al. 1999). The MS-AFLP technique uses two isoschizomer restriction enzymes (MspI and Hpall) that recognize the same restriction site, cutting in the same location (5’-CACGG-3’), but have different sensitivities to the presence of cytosine methylation. By comparing the presence/absence of fragments produced by both enzymes, the methylation status at each restriction site can be evaluated. In addition, genetic poly- morphisms can be evaluated between individuals that lack the restriction site in both digestions as an indication of a change in nucleotide sequence. To distinguish whether divergence of these two species was initiated by epigenetic, rather than genetic variation, we addressed two questions using a MS-AFLP technique. First, we asked to what extent epigenetic and/or genetic differentiation predicted species assignment. Second, we analyzed whether the patterns of genome- wide DNA methylation were correlated with genetic variation. Materials and methods Specimen collection.—We utilized a total of 94 males collected in the summers of 2006, 2013, and 2014 from eight grassland sites around the state of Missouri (Suppl. material 1: Table S1). We used males’ calls to localize individuals in the field and collect- ed them by hand after sunset. We identified males in the field as G. NEY AND J. SCHUL belonging to members of the target group through their call and morphological features, including cone pigmentation and body size, prior to collection (Walker et al. 1973). As the two species are morphologically cryptic, we identified individuals as either N. bivocatus or N. robustus based on their divergent call characteristics as described below. Individuals were collected from among eight grassland sites within Missouri and from within three different years (Suppl. material 1: Table S1). Thus, a range of epigenetic patterns, expressed within variable environmental conditions, was included in each species sample. Fixed differences observed between the two species are therefore likely to be species specific rather than the result of differences in environmentally mediated methylation patterns. Call recordings. —We recorded male calls in 2006 and 2013 within three days of collection using an Audiotechnica ATR 55 micro- phone and a Marantz PMD-671 solid-state recorder (16 bit, 48 kHz sampling rate). Recordings were made outdoors with males placed in individual mesh cages (approximately 10x20x10 cm) spaced at least 3 m apart. In 2014, we recorded male calls in the field immediately preceding collection using a Tascam DR-40 lin- ear PCM recorder (16 bit, 48 kHz sampling rate). Ambient tem- peratures ranged from 22-28°C. Temporal call analysis. —Call temporal patterns were analyzed and species assignments made as described in Ney and Schul (2017). In short, we marked each sound pulse produced during a single closing movement of the wings (Walker 1975) using custom soft- ware. The recordings were rectified and the envelope extracted with a temporal resolution of 0.125 ms. Pulses were detected automatically using a threshold-based algorithm and manually checked before saving the data to text-files. We analyzed about 2 s of each male’s call containing 150-250 pulses. Further analyses were conducted on the text files in MS Excel. Male calls of N. robustus have a single pulse rate of about 200/s, equivalent to a pulse period of about 5 ms (Walker 1975). Females recognize this pattern by the absence of silent gaps longer than about 2 ms (Deily and Schul 2004). In contrast, male calls of N. bivocatus have two alternating pulse periods of about 4.5 ms and 7 ms (Deily and Schul 2004), resulting in a ‘galloping rhythm’ or ‘pulse pairs’. Female N. bivocatus recognize this pattern by the rate of the pulse pairs: pulse pair rates around 87/s are attractive and attractiveness decreases towards faster and slower rates. For the pulse pairs to be detectable, the two alternating pulse periods must differ sufficiently (J. Schul, unpublished work). We quantified the ratio of the means of the alternating pulse periods (longer pp / shorter pp). For the single pulse calls of N. ro- bustus, this should result in values close to one, while in N. bivoca- tus significantly larger values result (Bush and Schul 2010). Among the calls of 94 males analyzed, there were two distinct groups of values, one between 1 and 1.18 and a second group with ratios >1.38. A natural break occurred in the data between 1.18 and 1.38. We classified ratios <1.18 as ‘N. robustus’, and >1.38 as ‘N. bivocatus’ in call type. Spectral call analysis. —Call spectra were analyzed and species as- signments made as described in Ney and Schul (2017). In short, we analyzed the amplitude spectrum of male calls using a fast Fou- rier transformation (FFT, Hamming window, frame length 256) in Audacity v.1.3 (Audacity Team 2008). Spectra were averaged over 2 s of the call. The main energy in the spectrum of Neoconocephalus calls is concentrated in a narrow band between 7 and 15 kHz with JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) G. NEY AND J. SCHUL the peak frequency differing among species (Schul and Patterson 2003). We measured the center frequency of this low frequency band as the geometric mean of the upper and lower cut-off fre- quencies at -3 dB from the peak amplitude. In N. bivocatus, the center frequencies of this band have a mean of about 10 kHz among individuals, ranging from 7 to 15 kHz. Fe- males have little spectral selectivity in this frequency range (Deily and Schul 2006). The center frequencies of calls with a pulse peri- od ratio of >1.38 that also fell into this range of center frequencies were Classified as N. bivocatus. One individual had a pulse period ratio greater than 1.38 and a center frequency lower than 7 kHz; this male was classified as an intermediate caller, as its call charac- teristics differed from that of both species. In N. robustus, the low frequency band is narrower than in N. bivocatus and is typically limited to 10 kHz and below (Schul and Patterson 2003). Indeed, frequencies well above 10 kHz have an inhibitory effect on female phonotaxis (Deily and Schul 2006). Of the individuals with pulse period ratios <1.18, center frequencies of all but five individuals clustered between 6 and 9 kHz. The five remaining individuals had center frequencies ranging from 9.7 to 11.1 kHz, which would be less attractive for N. robustus females. These individuals were classified as possessing an intermediate call type. All intermediate callers were removed from the analysis of epigenetic and genetic differentiation. Molecular analysis.—We aim here to analyze the differences of trans-generational methylation patterns between two species. Therefore, we selected tissue that is likely to have little tissue-spe- cific methylation differences between the two species. We collected tissue from hind femurs (mostly muscle and cuticle) that should be under similar selection in both species and thus little differen- tial tissue-specific methylation. Using leg tissue also avoids the risk of contamination from GI tract microbiota. We removed the hind femurs of collected males and placed them in 95% EtOH for DNA preservation. We later extracted DNA from the hind femurs using the DNeasy Blood & Tissue Kit (Qia- gen Inc., Valencia, CA, USA). DNA quantification was performed on each sample by spectrophotometry (NanoDrop 1000, Thermo Scientific, Wilmington, DE). Genomic DNA was stored at -80°C prior to molecular analysis. We used a MS-AFLP assay modified from Xu et al. (2000). DNA (55 ng) from each sample was digested in two separate dou- ble digest reactions (EcoRI/HpalI and EcoRI/MsplI). EcoRI selec- tively cuts the sequence 5’-GAATTC-3”. Hpall and MspI are isoschi- zomers, meaning they selectively cut at the sequence (5’-CCGG- 3’), but differ in their sensitivity to cytosine methylation at those sites. Hpall will only cut if the external cytosine is hemimethyl- ated. MspI cuts when the internal cytosine is either hemi or fully methylated. Both enzymes will cut if the target site is completely unmethylated. Using this method, we evaluated the CpG methyla- tion of restriction sites by comparing the fragments produced by both digests. Digestion and ligation were carried out together to prevent regeneration of restriction sites. Synthetic double stranded DNA adaptors (Xu et al. 2000) were ligated to the cleaved ends of re- striction sites. The EcoRI/Hpall and EcoRI/MspI digestion/ligation reaction (11 ul final volume) is comprised of 1.1 pL 10X CutSmart- ™ buffer (New England Biolabs, USA), 0.55 pg/pL of Bovine Se- rum Albumin (BSA) (New England Biolabs), 0.3 pL water, 5 U of EcoRI HE, 1 U of either Hpall or MspI, 1 U T4 ligase (New England Biolabs), 1 pL ATP (10 pM), EcoRI adaptors (5 pM), either Hpall or MspI adaptors (50 pM), and 5.5 pL genomic DNA (10 ng/ pL). 13 Table 1. Loci produced by selective primer combinations used in the MS-AFLP analysis. Shown are the primer pair combinations, the number of scored bands per primer pair, and the number of those bands classified as polymorphic methylation sensitive loci (MSL), and polymorphic genetic loci. EcoRI Msp1/Hpall Bands MSL Genetic loci -AAC -ATC 277 183 265 -AGC -AAT 22 49 102 -AGC -ATC 318 170 301 Total 822 402 668 The reaction was incubated at 37°C for 2 hours. Preselective PCR was conducted with 1:10 dilute digestion/ligation products, the Eco+A (5’-GACTGCGTACCAATTCA-3’), and the HpallI/MspI+A (5’-GATGAGTCTIAGAACGGA-3’) primers using thermocycler set- tings as described in Snyder et al. (2009). We performed selective PCR independently with three primer pairs (Table 1) and 1:100 dilute preselective PCR products from both the Hpall and MspI reactions. Fluorescently labeled Eco primers (Eco+AAC (6GFAM), Eco+AGC (PET)) were used in selective PCR (as described in Sny- der et al. 2009) and the products multiplexed and diluted to pro- duce a 1:10 final dilution of each product. Fragments were separat- ed in an ABI 3730 genetic analyzer at the DNA Core Facility, Uni- versity of Missouri. We called MS-AFLP bands using GeneMarker v.1.6 (Hulce et al. 2011) using an automated peak-calling scheme (as described in Holland et al. 2008) and called alleles between 75-500 bp with a minimum peak intensity of 50. Data analysis. —We obtained presence/absence fragment data for both EcoRI/Hpall and EcoRI/MspI datasets from GeneMarker (Hulce et al. 2011). Presence of a fragment in both the EcoRI/ Hpall and EcoRI/MspI digestions (1/1) was described as unmeth- ylated. The presence of a fragment in one digestion but not the other (0/1 or 1/0) was defined as methylated (either hemimethyl- ated or internal cytosine methylated, respectively). If fragments were absent in both digestions (0/0) the loci were considered un- informative as this state could be due to either full methylation at the target site (methylation of both the inner and outer cytosine) or the absence of the fragment due to variation in the nucleotide sequence between individuals (Schulz et al. 2013, Fulneéek and Kovarik 2014). The methylation sensitivity of each locus was identified us- ing the MSAP Package (Pérez-Figueroa 2013) in R v.3.1.2 (R Core Team 2015). Loci were classified as methylation sensitive loci (MSL) if there was evidence of methylation in at least 5% of the sampled individuals at that locus. Genetic markers were extracted from the MS-AFLP loci. Fragments that were present in one or both MspI and Hpall analyses were scored as present. Fragments that were absent in both reactions were scored as absent. This method of using MS-AFLP loci to estimate genetic parameters has been found to produce similar results to standard AFLP markers (Smith et al. 2016). Only polymorphic MSL and genetic loci that differed among sampled individuals were used in further analy- ses. We evaluated variation for each epigenetic (MSL) and genetic locus individually then calculated the mean diversity of MSL and genetic loci using Shannon’s diversity index (I). We compared the relative frequency of total methylation (internal cytosine meth- ylated and hemimethylated fragments) and unmethylated MSL between species. Significant variation between species was tested using a Mann-Whitney U test. JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) 14 We tested for differences in MSL and genetic differentiation us- ing two-way analyses of molecular variance (AMOVA; Excoffier et al. 1992) that grouped individuals again by species. Significance of the test statistic’s (O,,) deviation from zero was estimated based on 1000 permutations of individuals among groups. Principal co- ordinate analyses (PCoA) of both epigenetic and genetic loci were performed using the R stats package v.3.1.2 (cmdscale; R Core Team 2015) to visualize the Euclidean distance between species. A population structure analysis was implemented for all individuals using the program STRUCTURE v.2.3.3 (Pritchard et al. 2010). The admixture model was used, allele frequencies correlated, with a run length of 100,000 (Burnin = 10,000) for 10 replicates each of K = 1-10 (genetic clusters). The estimate of the most well supported K was calculated as described in Evanno et al. (2005) and imple- mented in Structure Harvester v.0.6.94 (Earl 2012). The program Clumpp v.1.1.2 (Jakobsson and Rosenberg 2007) was used to align the 10 repetitions of the most well supported number of clusters. If similar signals exist for the MSL and genetic structure then epigenetic and genetic distance may be significantly correlated. We estimated the Euclidean distance between individuals for both epigenetic and genetic datasets using the R stats package v.3.1.2 (R Core Team 2015) and compared the distance matrices using a Mantel test (Mantel 1967; ade4 package v.1.7.2; Dray et al. 2007). We utilized 10,000 permutations of the design matrix to deter- mine the significance of the correlation coefficient. Results Genome-wide variation in methylation.—We included 94 males in our analysis. We determined the males’ species assignments based on the temporal and spectral frequency call preferences of N. ro- bustus and N. bivocatus. The call analysis classified the 94 males as 31 N. bivocatus, 57 N. robustus, and 6 having an intermediate call type; five of the intermediates had the N. robustus temporal pat- tern and one the N. bivocatus pattern, with frequency spectra that fell outside of these respective species’ specific pattern (Fig. 1). In- dividuals possessing intermediate call phenotypes were removed from further analyses. The MS-AFLP analysis yielded 227, 277, and 318 fragments in each of three selective PCR reactions (Table 1). Of these fragments, a total of 364 were polymorphic for their meth- ylation status among sampled individuals and were classified as MSL. A loci-calling scheme, utilizing the MS-AFLP dataset, allowed for the inference of the genetic state of fragments among individu- als, similar to a traditional AFLP analysis. In total, 668 polymor- phic genetic loci were produced from 822 MS-AFLP fragments. We compared whether species differed significantly in their proportion of genome-wide methylated sites. The relative frequen- cy of genome-wide methylation (combined hemimethylation and internal cytosine methylation) showed low (Fig. 2; N. robustus = 0.272; N. bivocatus = 0.240) but significant variation between spe- cies (Mann-Whitney U, W = 587, p = 0.001243). Epigenetic and genetic structure.—We examined epigenetic and ge- netic diversity as they relate to species assignment. The within-spe- cies epigenetic Shannon diversity index was 5.0616 + 0.2054 and 5.0094 + 0.2087, within N. robustus and N. bivocatus, respectively. N. robustus and N. bivocatus’ genetic Shannon diversity indexes were 5.2802 + 0.2242 and 5.2179 + 0.1943, respectively. Diver- sity in the genetic loci, as measured in this study, was significantly greater than that of the MSL (epigenetic) diversity measured in both species (Wilcoxon rank sum test; N. robustus, W = 2760, p < 0.0001; N. bivocatus, W = 779, p = 0.0001). G. NEY AND J. SCHUL 14, = 134 ° = x 12) e S11 *e ee e e $10 tM ai a a e e e | e uc 9 : eh - | 2 5 ¥ Cc 5 7 | N. bivocatus e 6- 1 144 #12 13 14 #15 416 «447 ~«218~«19 PP Ratio (PP2/PP1) Fig. 1. Species assignment based on call pulse period ratio and center frequency. Labeled boxes indicate the calls classified as N. robustus and N. bivocatus. Individuals that fall outside of species classifications were removed from further epigenetic and genetic analyses (as described in Ney and Schul 2017). 100 > 90 - @ Methylated O Unmethylated cet 80 - 0 Absence of restriction site Oo L& arin FA een cD) 2 Fa 60 - ** S — ” I l Cc = 504 & ~_ = o 4 2 40 ro) S 30- = ro oO a. 207 10 - Q 7— = N. robustus N. bivocatus Fig. 2. Comparison of genome-wide methylation levels between species. Mann-Whitney U test; p<0.005 (**). Between-species sig- nificant variation is in total methylation (internal cytosine methyl- ated and hemimethylated fragments). Epigenetic and genetic structure was evaluated using a two- level AMOVA. The between-species epigenetic divergence was ®,, = 0.0504 (P < 0.0001). This is slightly less than one-third of the genet- ic divergence observed between species, O,, = 0.1591 (P < 0.0001). Greater genetic variation between species would suggest that genetic mechanisms are underlying species differentiation (Table 2). JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) G. NEY AND J. SCHUL Table 2. Two-level AMOVA of MSL or genetic loci produced from MS-AFLP markers among populations grouped by species as- signment. Included are genetic and epigenetic variance between groups, ®,,, and the corresponding p value. Between group Within group ©, pvalue variance variance MSL (epigenetic) 3.429 (5.04%) 64.54(94.96%) 0.0504 <0.0001 Genetic loci 15.38 (15.91%) 81.28 (84.09%) 0.1591 <0.0001 The principal coordinate analyses (PCoA), calculated using MSL and genetic profiles, showed that genetic variance was smaller than epigenetic variance within species (Fig. 3). The epigenetic PCoA A C1 (11.4%) 15 between species elucidated little meaningful differentiation among groups, with ellipses (95% confidence intervals) overlapping al- most entirely. The PCoA of genetic data between species explained 17.7% of the genetic variance in the first two coordinates (Fig. 3A) and showed moderate variation between individuals based on call phenotype, with species showing divergence in Euclidean space. Three N. robustus individuals clearly fell into the N. bivocatus cluster within the genetic PCoA (Fig. 3A), indicating a possible mismatch between phenotypic and genotypic assignments. The Bayesian analysis of genetic structure revealed the best- supported number of genetic clusters to be K = 2, based on AK values (Fig. 4B; Evanno et al. 2005). These two genetic clusters aligned closely to species boundaries (Fig. 4A). As in the PCoA, wv N CeSNAEA SO Zo OBS a te aoe ~~ _ Ol CN © Ty Ss © +5 0 5 C1 (10.5%) Fig. 3. PCoA of N. robustus and N. bivocatus utilizing genetic (A) and epigenetic (B) data. Plotted are the two most informative principal components calculated for the genetic and epigenetic loci datasets, as derived from the MS-AFLP fragment analysis. A. Genetic Euclid- ean distance with individuals grouped by species assignment. B. Epigenetic Euclidean distance with individuals grouped by species as- signment. Group labels show the centroid of the points for each group. The long axis of the ellipse represents the direction of maximum dispersion and the short axis the direction of minimum dispersion. > Genetic N. bivocatus ra Epigenetic N. bivocatus N. robustus D N. robustus Fig. 4. Consensus shared ancestry population structure for epigenetic and genetic loci. A. and C. Bar plots using MS-AFLP loci to estimate ge- netic (A) and epigenetic (C) structure among N. robustus and N. bivocatus using the software package STRUCTURE. B. and D. Delta K graphs for K = 1-10 genetic clusters showing moderate support for K = 2 genetic clusters (B) and low support for K = 4 epigenetic clusters (D). JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) 16 ice) “ : o ea mee ca ee a BN — 7 wer. nD 2 po Mantel R = 0.8448 > rar p = 0.0001 wT a a 8 10 12 14 16 18 Genetic distance Fig. 5. Scatterplot of between-individual Euclidean genetic and epigenetic distance showing significant positive correlation be- tween genetic and epigenetic differentiation. The correlation was tested using a Mantel test and 10,000 permutations of the design matrix to determine significance. the same three individuals identified as N. robustus were assigned primarily to the N. bivocatus genetic cluster, indicating that they possess the N. robustus call type but a N. bivocatus genotype. As in the PCoA, the epigenetic structure analysis did not show support for significant epigenetic (MSL) structure between species groups. AK values did not rise above 30 for any permutation in the num- ber of clusters (Fig. 4D), an indication of low support for epige- netic structure. In addition, the best supported number of clusters, K = 4, showed no grouping of epigenetic clusters by species assign- ment (Fig. 4C). While significant epigenetic variation was observed between species, significant genetic variation was also detected. We inves- tigated the correlation between inter-individual genetic and epi- genetic Euclidean distance. Between-individual epigenetic and ge- netic distance showed a strong positive correlation (Fig. 5, Mantel R = 0.8448, p = 0.0001). Discussion We found significant differentiation among species in both epigenetic and genetic markers. Genetic differentiation, however, was larger than epigenetic differentiation between species. The AMOVA and PCoA analyses indicated that methylation patterns varied among individuals but showed little differentiation be- tween species. Epigenetic distance among individuals correlated with inter-individual genetic differentiation, suggesting that the low levels of epigenetic differentiation between species have been pulled along by genetic differentiation. Genetic differentiation.—Genetic differentiation between species was low but significant, as would be expected between two closely related taxa (Snyder et al. 2009, Ney and Schul 2017). Genetic structure aligned with species boundaries, with the exception of the three phenotypically N. robustus individuals that showed a G. NEY AND J. SCHUL majority assignment to the primarily N. bivocatus genotypic clus- ter. These mismatch genotype/phenotype assignments were con- firmed for these three individuals in both the genetic structure analysis (Fig. 4A) as well as the genetic PCoA (Fig. 3A). Similar occurrences of mismatched genotype/phenotype individuals have been found between N. robustus and N. bivocatus previously (Fred- erick 2013, Ney and Schul 2017). These mismatched genotype/ phenotype individuals are likely the result of recent hybridiza- tion events. A severe drought during the sampling period (Ney and Schul 2017) led locally to sharp population declines of one or both species. Both reductions in population sizes and environ- mental disturbances can increase rates of hybridization between closely related species (Lamont et al. 2003, Seehausen 2006). Epigenetic variation.—N. robustus showed a significantly higher lev- el of gnome wide methylation than N. bivocatus (Fig. 2). This may be due to variation in the species’ habitat preference. Variable en- vironmental conditions can cause shifts in epigenetics, observed between diverging natural habitats as well as between natural and altered habitat types (Gao et al. 2010). N. robustus, unlike N. bivo- catus, is a grassland generalist and may show epigenetic patterns resulting from exposure to a more variable set of environmental conditions during a lifetime and across generations. Because of this, N. robustus may possess a larger repertoire of adaptive genes for the varied environmental conditions it may utilize. The higher level of methylation found in N. robustus may therefore be the re- sult of higher levels of methylation required to silence the large repertoire of adaptive genes when not in use. While epigenetic differentiation between species was signifi- cant, it remained lower than genetic differentiation. This study did not show support for the epigenetic regulation of species-specific phenotypes; however, this does not eliminate a possible epige- netic mechanism underlying phenotypic differentiation. While the MS-AFLP technique has many benefits, there are also some inherent limitations to its application. For example, MS-AFLPs can underestimate genome-wide levels of methylation (Fulneéek and Kovarik 2014). As a genome-wide scanning method MS-AFLP only detects methylation at 5'-CCGG-3' sites. In addition, it is unable to discriminate between full methylation states (hypermethylation of both cytosines) and sequence variation at or near the restric- tion site (Xiong et al. 1999). Phenotypic variation controlled by a smaller number of methylated sites is unlikely to be detected. A similar investigation into gnome-wide methylation found no sig- nificant differentiation in methylation between the solitarious and gregarious phases of Locusta migratoria (Robinson et al. 2011), de- spite divergent phenotypes and the identification of differentially expressed genes between morphs (Kang et al. 2004). Thus, this technique is most valuable as a first step in investigating genome- wide methylation patterns. Correlation between epigenetic and genetic diversity.—MS-AFLP vari- ation often shows correlations with genetic variation (Liu et al. 2012). In our study, individuals showing greater genetic variation on average showed greater epigenetic variation as well (Fig. 5). Changes in DNA methylation over time are correlated with ge- netic relatedness and suggest that DNA methylation maintenance may be under genetic control (Bjornsson et al. 2008). In addition, genetic variation in retrotransposons (i.e. their presence/absence) could affect the methylation state of retroelements, requiring more methylation if more repetitive elements are present (Michaud et al. 1994). Repetitive elements have been estimated to make up thirty percent of the genome of the orthopteran L. migratoria (Wil- JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) G. NEY AND J. SCHUL more and Brown 1975). In addition, the retrotransposon SINE is differentially methylated between the solitarious and gregarious phases of the species (Guo et al. 2010). Mechanisms of neutral evolution could also account for the correlation between epigenetic and genetic variation. In the event of substantial gene flow between species, strong divergent selec- tion would be needed to maintain divergent epigenetic variation between species. Gene flow between groups would reduce differ- entiation accumulated via drift. Evidence from this study, however, suggests that genetic differentiation between species is significant and gene flow therefore relatively low (Table 2). Stochastic process- es of drift then could allow epigenetic patterns to diverge between species in parallel, resulting in correlations between epigenetic and genetic variation without a causal link (Richards et al. 2010). Conclusions We hypothesized that the phenotypic variation observed be- tween N. robustus and N. bivocatus was the result of an epigenetic- mediated mechanism of species differentiation that led to genetic isolation. While we found clear evidence of genomic methylation in Neoconocephalus, our findings did not support a role for meth- ylation in species isolation. Both the lower level of epigenetic dif- ferentiations and the correlation between inter-individual epige- netic and genetic diversity support the alternative hypothesis, i.e. that differences in methylation patterns between species evolved in response to genetic variation. Epigenetics may still play a key role in phenotypic differentiation within Neoconocephalus katydids through the differential regulation of a relatively small number of genes of large effect; a mechanism not detectable with the meth- ods used here. Further work identifying differentially expressed genes between call types could allow for the targeted analysis of methylation patterns at these sites. Acknowledgements We thank Katy Frederick and Nathan Harness for assistance with specimen collection and recording males. We thank Kim Hunter, Katy Frederick, and Megan Murphy for valuable feedback on the manuscript. References Audacity Team A (2008) Audacity (version 1.3) [computer program]. htt- ps://old.audacityteam.org/ Beckers OM, Schul J (2008) Developmental plasticity of mating calls ena- bles acoustic communication in diverse environments. Proceedings of the Royal Society B: Biological Sciences 275: 1243-1248. https:// doi.org/10.1098/rspb.2007.1765 Beckers OM, Schul J (2010) Female adaptation to developmental plasticity in male calling behavior. Behavioral Ecology and Sociobiology 64: 1279-1290. https://doi.org/10.1007/s00265-010-0942-z Bente H, Scheid OM (2017) Epigenetic contribution to diversification. Proceedings of the National Academy of Sciences 114: 3558-3560. https://doi.org/10.1073/pnas.1702748114 Berger SL, Kouzarides T, Shiekhattar R, Shilatifard A (2009) An operational definition of epigenetics. Genes & Development 23: 781-783. https:// doi.org/10.1101/gad.1787609 Bjornsson HT, Sigurdsson MI, Fallin MD, Irizarry RA, Aspelund T, Cui H, Yu W, Rongione MA, Ekstr6m TY, Harris TB, Launer LJ, Eiriksdottir G, Leppert ME, Sapienza C, Gudnason V, Feinberg AP (2008) Intra-indi- vidual change over time in DNA methylation with familial clustering. JAMA 299: 2877-2883. https://doi.org/10.1001/jama.299.24.2877 17 Blevins T, Wang J, Pflieger D, Pontvianne F, Pikaard CS (2017) Hybrid incom- patibility caused by an epiallele. Proceedings of the National Academy of Sciences 114: 3702-3707. https://doi.org/10.1073/pnas. 1700368114 Boyes J, Bird A (1991) DNA methylation inhibits transcription indirectly via a methyl-CpG binding protein. Cell 64: 1123-1134. https://doi. org/10.1016/0092-8674(91)90267-3 Braam J, Davis RW (1990) Rain-, wind-, and touch-induced expression of calmodulin and calmodulin-related genes in Arabidopsis. Cell 60: 357-364. https://doi.org/10.1016/0092-8674(90)90587-5 Burggren WW, Crews D (2014) Epigenetics in comparative biology: Why we should pay attention. Integrative and Comparative Biology 54: 7-20. https://doi.org/10.1093/icb/icu013 Bush SL, Schul J (2010) Evolution of novel signal traits in the absence of female preferences in Neoconocephalus katydids (Orthoptera, Tet- tigoniidae). PLoS ONE 5: e12457. https://doi.org/10.1371/journal. pone.0012457 Deily JA, Schul J (2004) Recognition of calls with exceptionally fast pulse rates: female phonotaxis in the genus Neoconocephalus (Orthoptera: Tettigoniidae). Journal of Experimental Biology 207: 3523-3529. https://doi.org/10.1242/jeb.01179 Deily JA, Schul J (2006) Spectral selectivity during phonotaxis: a comparative study in Neoconocephalus (Orthoptera: Tettigoniidae). Journal of Experi- mental Biology 209: 1757-1764. https://doi.org/10.1242/jeb.02189 Deily JA, Schul J (2009) Selective phonotaxis in Neoconocephalus nebras- censis (Orthoptera: Tettigoniidae): call recognition at two temporal scales. Journal of Comparative Physiology A 195: 31-37. https://doi. org/10.1007/s00359-008-0379-2 Dray S, Dufour AB, Chessel D (2007) The ade4 package-II: Two-table and K-table methods. R News 7: 47-52. https://pbil.univ-lyon1.fr/ade4/ article/rnews2/rmews2.pdf Earl DA (2012) Structure harvester: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Con- servation Genetics Resources 4: 359-361. https://doi.org/10.1007/ $12686-011-9548-7 Elango N, Hunt BG, Goodisman MA, Soojin VY (2009) DNA methyla- tion is widespread and associated with differential gene expression in castes of the honeybee, Apis mellifera. Proceedings of the National Academy of Sciences 106: 11206-11211. https://doi.org/10.1073/ pnas.0900301106 Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clus- ters of individuals using the software structure: a simulation study. Molecular Ecology 14: 2611-2620. https://doi.org/10.1111/j.1365- 294X.2005.02553.x Excoffier L, Smouse PE, Quattro JM (1992) Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131: 479- 491. http://www.genetics.org/content/131/2/479.short Field LM, Devonshire AL, Ffrench-Constant RH, Forde BG (1989) Changes in DNA methylation are associated with loss of insecticide resistance in the peach-potato aphid Myzus persicae (Sulz.). FEBS Letters 243: 323-327. https://doi.org/10.1016/0014-5793(89)80154-1 Frederick FH (2013) Investigating an adaptive radiation in temperate Neo- conocephalus (Orthoptera: Tettigoniidae). PhD Thesis, Columbia, Mis- souri: University of Missouri. https://mospace.umsystem.edu/xmlui/ handle/10355/44039 Fulneéek J, Kovarik A (2014) How to interpret Methylation Sensitive Am- plified Polymorphism (MSAP) profiles? BMC Genetics 15: 2. https:// doi.org/10.1186/1471-2156-15-2 Gao L, Geng Y, Li BO, Chen J, Yang JI (2010) Genome-wide DNA meth- ylation alterations of Alternanthera philoxeroides in natural and ma- nipulated habitats: implications for epigenetic regulation of rapid responses to environmental fluctuation and phenotypic variation. Plant, Cell & Environment 33: 1820-1827. https://doi.org/10.1111/ j.1365-3040.2010.02186.x Gerhardt HC, Huber F (2002) Acoustic Communication in Insects and Anurans: Common Problems and Diverse Solutions. University of Chicago Press. JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) 18 Greenfield MD (1990) Evolution of acoustic communication in the genus Neoconocephalus. In: Bailey WJ, Rentz DCF (Eds) The Tettigoniidae: Biology, Systematics and Evolution. Srpinger-Verlag, Berlin, 71-97. https://doi.org/10.1007/978-3-662-02592-5_5 Guo W, Wang XH, Zhao DJ, Yang PC, Kang L (2010) Molecular cloning and temporal-spatial expression of I element in gregarious and soli- tary locusts. Journal of Insect Physiology 56: 943-948. https://doi. org/10.1016/j.jinsphys.2010.05.007 Herrera CM, Bazaga P (2010) Epigenetic differentiation and relationship to adaptive genetic divergence in discrete populations of the violet Viola cazorlensis. New Phytologist 187: 867-876. https://doi.org/10.1111/ j.1469-8137.2010.03298.x Holland BR, Clarke AC, Meudt HM (2008) Optimizing automated AFLP scoring parameters to improve phylogenetic resolution. Systematic Biology 57: 347-366. https://doi.org/10.1080/10635150802044037 Hulce D, Li X, Snyder-Leiby T, Johathan Liu CS (2011) GeneMarker geno- typing software: Tools to increase the statistical power of DNA frag- ment analysis. Journal of Biomolecular Techniques 22: S35-S36. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3186482/ Jablonka E, Lamb MJ (1998) Epigenetic inheritance in evolution. Jour- nal of Evolutionary Biology 11: 159. https://doi.org/10.1007/ s000360050073 Jablonka E, Lamb MJ (1999) Epigenetic Inheritance and Evolution: The Lamarckian Dimension. Oxford University Press on Demand. Jakobsson M, Rosenberg NA (2007) CLUMPP: a cluster matching and per- mutation program for dealing with label switching and multimodal- ity in analysis of population structure. Bioinformatics 23: 1801-1806. https://doi.org/10.1093/bioinformatics/btm233 Kang L, Chen X, Zhou Y, Liu B, Zheng W, Li R, Wang J, Yu J (2004) The analysis of large-scale gene expression correlated to the phase changes of the migratory locust. Proceedings of the National Academy of Sciences 101: 17611-17615. https://doi.org/10.1073/ pnas.0407753101 Lamont BB, He T, Enright NJ, Krauss SL, Miller BP (2003) Anthropogenic disturbance promotes hybridization between Banksia species by al- tering their biology. Journal of Evolutionary Biology 16: 551-557. https://doi.org/10.1046/j.1420-9101.2003.00548.x Ledon-Rettig CC, Richards CL, Martin LB (2012) Epigenetics for behavioral ecologists. Behavioral Ecology 24: 311-324. https://doi.org/10.1093/ beheco/ars145 Liu S, Sun K, Jiang T, Ho JP, Liu B, Feng J (2012) Natural epigenetic varia- tion in the female great roundleaf bat (Hipposideros armiger) popula- tions. Molecular Genetics and Genomics 287: 643-650. https://doi. org/10.1007/s00438-012-0704-x Mantel N (1967) The detection of disease clustering and a generalized re- gression approach. Cancer Research 27: 209-220. https://www.ncbi. nlm.nih.gov/pubmed/6018555 Michaud EJ, Van Vugt MJ, Bultman SJ, Sweet HO, Davisson MT, Woychik RP (1994) Differential expression of a new dominant agouti allele (Aiapy) is correlated with methylation state and is influenced by parental lineage. Genes & Development 8: 1463-1472. https://doi. org/10.1101/gad.8.12.1463 Ney G, Schul J (2017) Low genetic differentiation between populations of an endemic prairie katydid despite habitat loss and fragmenta- tion. Conservation Genetics 18: 1389-1401. https://doi.org/10.1007/ $10592-017-0987-x Pal C, Miklos I (1999) Epigenetic inheritance, genetic assimilation and speciation. Journal of Theoretical Biology 200: 19-37. https://doi. org/10.1006/jtbi.1999.0974 Pérez-Figueroa A (2013) MSAP: a tool for the statistical analysis of meth- ylation-sensitive amplified polymorphism data. Molecular Ecology Resources 13: 522-527. https://doi.org/10.1111/1755-0998.12064 Pritchard JK, Wen X, Falush D (2010) Documentation for STRUCTURE software, version 2.3. University of Chicago, Chicago, IL. http://bur- fordreiskind.com/wp-content/uploads/Structure_Manual_doc.pdf R Core Team (2015) R: A Language and Environment for Statistical Com- puting (Version 3.1. 2): R Foundation for Statistical Computing. Vi- enna, Austria. http://www.R-project.org G. NEY AND J. SCHUL Rapp RA, Wendel JF (2005) Epigenetics and plant evolution. New Phytolo- gist 168: 81-91. https://doi.org/10.1111/j.1469-8137.2005.01491.x Richards CL, Bossdorf O, Verhoeven KJ (2010) Understanding natural epigenetic variation. New Phytologist 187: 562-564. https://doi. org/10.1111/j.1469-8137.2010.03369.x Robinson KL, Tohidi-Esfahani D, Lo N, Simpson SJ, Sword GA (2011) Evi- dence for widespread genomic methylation in the migratory locust, Locusta migratoria (Orthoptera: Acrididae). PLoS ONE 6: e28167. https://doi.org/10.1371/journal.pone.0028167 Sarkar S, Rao SR, Gupta VS, Hendre RR (1992) 5-Methylcytosine content in Gryllotalpa fossor (Orthoptera). Genome 35: 163-166. https://doi. org/10.1139/g92-026 Schul J, Bush SL, Frederick KH (2014) Evolution of call patterns and pattern recognition mechanisms in Neoconocephalus katydids. In: Hedwig B (Ed.) Insect Hearing and Acoustic Communication. Springer-Verlag, Berlin, 167-183. https://doi.org/10.1007/978-3- 642-40462-7_10 Schul J, Patterson AC (2003) What determines the tuning of hearing organs and the frequency of calls? A comparative study in the ka- tydid genus Neoconocephalus (Orthoptera, Tettigoniidae). Journal of Experimental Biology 206: 141-152. https://doi.org/10.1242/ jeb.00070 Schulz B, Eckstein RL, Durka W (2013) Scoring and analysis of methyl- ation-sensitive amplification polymorphisms for epigenetic popula- tion studies. Molecular Ecology Resources 13: 642-653. https://doi. org/10.1111/1755-0998.12100 Seehausen O (2006) Conservation: Losing biodiversity by reverse spe- ciation. Current Biology 16: R334-R337. https://doi.org/10.1016/j. cub.2006.03.080 Shaposhnikov GK (1966) Origin and breakdown of reproductive isolation and the criterion of the species. Entomological Review 45: 1-18. Smith TA, Martin MD, Nguyen M, Mendelson TC (2016) Epigenetic diver- gence as a potential first step in darter speciation. Molecular Ecology 25: 1883-1894. https://doi.org/10.1111/mec.13561 Snyder RL, Frederick-Hudson KH, Schul J (2009) Molecular phylogenet- ics of the genus Neoconocephalus (Orthoptera, Tettigoniidae) and the evolution of temperate life histories. PLoS ONE 4: e7203. https://doi. org/10.1371/journal.pone.0007203 Triblehorn JD, Schul J (2009) Sensory-encoding differences contribute to species-specific call recognition mechanisms. Journal of Neurophysi- ology 102: 1348-1357. https://doi.org/10.1152/jn.91276.2008 Varriale A, Bernardi G (2006) DNA methylation in reptiles. Gene 385: 122-127. https://doi.org/10.1016/j.gene.2006.05.034 Walker TJ (1975) Stridulatory movements in eight species of Neoconocepha- lus (Tettigoniidae). Journal of Insect Physiology 21: 595-603. https:// doi.org/10.1016/0022-1910(75)90163-8 Walker TJ, Whitesell JJ, Alexander RD (1973) The robust conehead: two widespread sibling species (Orthoptera: Tettigoniidae: Neoconocepha- lus “robustus”). Ohio Journal of Science 73: 321-30. http://hdl.handle. net/1811/22182 Walsh TK, Brisson JA, Robertson HM, Gordon K, Jaubert-Possamai S, Tagu D, Edwards OR (2010) A functional DNA methylation system in the pea aphid, Acyrthosiphon pisum. Insect Molecular Biology 19: 215- 228. https://doi.org/10.1111/j.1365-2583.2009.00974.x West-Eberhard MJ (2005) Developmental plasticity and the origin of species differences. Proceedings of the National Academy of Sci- ences 102 (Supplement 1): 6543-6549. https://doi.org/10.1073/ pnas.0501844102 Wilmore PJ, Brown AK (1975) Molecular properties of orthopteran DNA. Chromosoma 51: 337-345. https://doi.org/10.1007/BF00326320 Xiong LZ, Xu CG, Saghai Maroof MA, Zhang Q (1999) Patterns of cytosine methylation in an elite rice hybrid and its parental lines, detected by a methylation-sensitive amplification polymorphism technique. Mo- lecular and General Genetics 261: 439-446. https://doi.org/10.1007/ s004380050986 Xu M, Li X, Korban SS (2000) AFLP-based detection of DNA methyla- tion. Plant Molecular Biology Reporter 18: 361-368. https://doi. org/10.1007/BF02825064 JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1) G. NEY AND J. SCHUL 19 Supplementary material 1 Author: Gideon Ney, Johannes Schul Data type: PDF file Explanation note: Table $1: Sample collection localities, locality coordinates, and number of each species sampled in each year (N. robustus / N. bivocatus). Copyright notice: This dataset is made available under the Open Database License _(http://opendatacommons.org/licenses/ odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/jor.28.28888.suppl1 Supplementary material 2 Author: Gideon Ney, Johannes Schul Data type: CSV file Explanation note: Matrix of MS-AFLP called fragments for all in- dividuals. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/ odbl/1.0/). The Open Database License (ODDbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited. Link: https://doi.org/10.3897/jor.28.28888.suppl2 JOURNAL OF ORTHOPTERA RESEARCH 2019, 28(1)