OR I G I N A L A R T I C L E Genomic regions underlying the species-specific mating songs of green lacewings Katherine L. Taylor1,2 | Elizabeth J. Wade3 | Marta M. Wells4 | Charles S. Henry1 1Department of Ecology and Evolutionary Biology, University of Connecticut, Storrs, Connecticut, USA 2Department of Entomology, University of Maryland, College Park, Maryland, USA 3Department of Natural Sciences and Mathematics, Curry College, Milton, Massachusetts, USA 4Department of Ecology and Evolutionary Biology, Yale University, New Haven, Connecticut, USA Correspondence Katherine L. Taylor, Department of Entomology, University of Maryland, College Park, Maryland, USA. Email: taylork@umd.edu Funding information University of Connecticut Abstract Rapid species radiations provide insight into the process of speciation and diversification. The radiation of Chrysoperla carnea-group lacewings seems to be driven, at least in part, by their species–specific pre-mating vibrational duets. We associated genetic markers from across the genome with courtship song period in the offspring of a laboratory cross between Chrysoperla plorabunda and Chrysoperla adamsi, two species primarily differ- entiated by their mating songs. Two genomic regions were strongly associated with the song period phenotype. Large regions of chromosomes one and two were associ- ated with song phenotype, as fewer recombination events occurred on these chromo- somes relative to the other autosomes. Candidate genes were identified by functional annotation of proteins from the C. carnea reference genome. The majority of genes that are associated with vibrational courtship signals in other insects were found within QTL for lacewing song phenotype. Together these findings suggest that decreased recombination may be acting to keep together loci important to reproductive isolation between these species. Using wild-caught individuals from both species, we identified signals of genomic divergence across the genome. We identified several candidate genes both in song-associated regions and near divergence outliers including nonA, fruitless, paralytic, period, and doublesex. Together these findings bring us one step closer to identifying the genomic basis of a mating song trait critical to the mainte- nance of species boundaries in green lacewings. K E YWORD S Chrysopidae, courtship behaviour, lacewings, Neuroptera, QTL mapping, recombination, reproductive isolation, speciation INTRODUCTION Species–specific mating signals are thought to facilitate divergence and speciation (Wilkins et al., 2013). Identifying and characterizing the loci underlying these traits important to reproductive isolation between species is fundamental to understanding how speciation pro- ceeds (Coyne, 1992; Templeton, 1981). Though the genetic underpin- nings of traits essential to reproductive isolation have been identified in some organisms, long-standing questions about the number, relative positions, and effect size of barrier loci persist (Nosil & Schluter, 2011; Seehausen et al., 2014; Wu & Ting, 2004). The genetic basis of acoustic signals key to reproductive isolation between species has been studied extensively in Drosophila. Males pro- duce species–specific mating signals with both wing and abdominal vibra- tions (Fabre et al., 2012; Mazzoni et al., 2013). Many genes have been identified that are necessary for the production of typical wing-produced song phenotype or for the perception of that song in Drosophila (Gleason, 2005; Greenspan & Ferveur, 2000). Some of the same genes, Received: 30 May 2022 Accepted: 21 October 2022 DOI: 10.1111/imb.12815 This is an open access article under the terms of the Creative Commons Attribution-NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes. © 2022 The Authors. Insect Molecular Biology published by John Wiley & Sons Ltd on behalf of Royal Entomological Society. Insect Mol Biol. 2023;32:79–85. wileyonlinelibrary.com/journal/imb 79 mailto:taylork@umd.edu http://creativecommons.org/licenses/by-nc/4.0/ http://wileyonlinelibrary.com/journal/imb http://crossmark.crossref.org/dialog/?doi=10.1111%2Fimb.12815&domain=pdf&date_stamp=2022-11-04 including fruitless and doublesex, are also associated with courtship abdominal vibration phenotype (Fabre et al., 2012). However, those genes and mutations known to modify mating signals in the laboratory may not necessarily be the genes that drive species–specific differences in songs. The genetic basis of species–specific acoustic mating signals has also been studied in crickets of the genus Laupala. Males produce signals using spe- cialized stridulatory organs prior to mating (Mullen & Shaw, 2014). Mating signals in Laupala crickets have been associated with quantitative trait loci (QTL) spread across the genome containing multiple promising candidate genes (Blankers et al., 2019; Xu & Shaw, 2021). When multiple genomic loci are responsible for any complex phe- notype, but especially for those traits important to reproduction isola- tion, strong linkage is expected between loci to keep together these co-adapted genes in the presence of gene flow (Schwander et al., 2014). Signal production and signal preference traits are often closely linked, for example in the wing pattern signalling in Heliconius butterflies and in mating song of Laupala crickets (Kronforst et al., 2006; Xu & Shaw, 2021). These so-called ‘supergenes’ might be expected to occur in genomic regions with decreased recombination, such as areas near centromeres, on sex chromosomes, or within chro- mosomal rearrangements (Noor et al., 2001; Schwander et al., 2014). In the rapid radiation of carnea-group green lacewings, species– specific premating duets isolate more than twenty morphologically cryptic species (Henry et al., 2013). Both male and female lacewings in the car- nea-group produce signals that are exchanged in an intricate duet prior to copulation. The signals are sexually monomorphic, with males and females producing nearly identical songs. These signals are thought, at least in part, to have driven the recent and rapid radiation of species in this group (Henry et al., 2013). Phenotypic information from laboratory crosses sug- gests that mating song phenotypes in the carnea-group are genetically controlled by a few loci of large effect (Henry et al., 2002), but to date, no loci underlying mating song phenotypes have been identified. There are multiple pairs of sympatric species that vary primarily in their courtship signals in the Chrysoperla carnea-group. One such species pair is Chrysoperla plorabunda and Chrysoperla adamsi, which perform courtship songs that are similar in structure except for a consistent major difference in one critical song feature, the volley period (Henry et al., 1993, 2002). The signals are composed of multiple bursts of low- frequency substrate-borne sound which are produced by vigorous but highly controlled shaking of the thorax and abdomen. Volley period is the length of a single burst of sound (Figure 1). These two species lack strong morphological or ecological differences and co-occur in western North America. (Henry et al., 1993, 2002). As larvae, both species are generalist predators of invertebrates, while as adults they feed on nectar and honeydew. Chrysoperla plorabunda is largely found on herbaceous plants, while C. adamsi can be found on herbaceous plants, shrubby veg- etation, and fruit trees. These authors have collected individuals of both species from the same fields at the same time and they will hybridize in laboratory conditions, yet field-caught hybrids have not been identified. Herein, we used quantitative trait mapping and genome scans to identify the genomic architecture underlying the volley period phenotype in two sympatric species in the Chrysoperla carnea-group. We identified genomic regions associated with this trait important to reproductive isolation using an F2 laboratory cross. With this approach, we were able to identify regions truly related to the phenotype; however, fine mapping to small genomic windows was limited by number of individuals in the cross, recombination frequency, marker density, and strength of the effect of the loci (Broman, 2001; Mackay, 2001). We narrowed our search based upon signals of genomic divergence between wild mem- bers of the same two species to reveal patterns related to selection or a lack of introgression near these traits of interest. We then used func- tional annotation to identify candidate genes and found patterns sug- gesting that decreased recombination in these song-associated regions may help to keep song-associated loci together. METHODS Collecting, laboratory rearing, phenotype recording Chrysoperla plorabunda (Fitch) and Chrysoperla adamsi Henry et al. cross parents were collected in western Oregon, United States in 2015. The wild-caught Chrysoperla adamsi female was collected from F I G U R E 1 F2 hybrid cross phenotypes. At the top of the figure are 12-second oscillograms of vibrational mating signals of the two parent species, with Chrysoperla plorabunda in yellow and Chrysoperla adamsi in purple. A single volley period for each of the parental species is indicated with a red bar under the oscillogram. Arrows on the oscillograms indicate where a duetting partner would insert a response. Below are the volley period histograms for the Chrysoperla plorabunda (n = 29), Chrysoperla adamsi (n = 20), F1 (n = 35), and F2 (n= 83) individuals from the laboratory cross. The phenotypes of the individual C. plorabunda and C. adamsi grandparents that founded the cross are indicated with dashed vertical lines. 80 TAYLOR ET AL. 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense shrubby vegetation near the top of Mary’s Peak, OR. The wild-caught Chrysoperla plorabunda male was collected from a clover field along a road in Benton, OR. Six additional wild-caught individuals of Chryso- perla adamsi and Chrysoperla plorabunda were collected from North Winters and Gates Canyon California, Guilford Connecticut, Moscow Idaho, and Benton and Marion Oregon between 2008 and 2016. All lacewings were laboratory-reared on a long day light cycle (16 hours light, 8 hours dark) at 25 � 1�C. Adults were kept in clear plastic 8 oz cups with a hole in the bottom for a cotton wick. Rearing cups were stacked inside a second clear plastic cup containing a water reserve and covered by a petri dish lid. A 2:1:1 mixture of Wheast, honey, and water was provided ad libitum for adults. Larvae from the family cross were reared in clear plastic 1 oz cups with snap-on lids. Sterile Ephes- tia kuehniella eggs were provided to the larvae ad libitum. Mating signals for all individuals were recorded in the laboratory at 25 � 1�C using an optical microphone and waveform software (see Henry et al., 2013). The wild-caught individuals were identified to spe- cies based upon mating signals. A single pair of wild-caught parents were induced to mate by confinement together in a rearing container. F1 hybrid offspring were reared to adulthood and brother–sister mated to produce F2 hybrids. The C. adamsi by C. plorabunda cross produced 78 F2 offspring that were successfully sequenced. Family cross-sequencing and bioinformatic analysis DNA was extracted using a Qiagen DNeasy blood and tissue kit from the two parents founding the line, 13 F1 individuals, and 78 F2 individ- uals. A reduced representation sequencing library of the family cross individuals was prepared by Floragenex with the PstI cutting enzyme and sequenced on an Illumina HiSeq at the University of Oregon. Single-end 150 base pair reads were demultiplexed with Stacks (v. 2.2; Rochette et al., 2019) and mapped against the reference assembly for the closely related species Chrysoperla carnea (GCF_905475395.1; Crowley, 2021) with BWA (v. 0.7.17; Li & Durbin, 2009), and genotypes were called with bcftools (v. 1.9; Danecek et al., 2021) mpileup. Geno- types were filtered with bcftools to include only biallelic single nucleo- tide polymorphisms (SNPs) with a quality score >30, a minor allele frequency >0.05, missing in less than 80% of individuals, lacking signifi- cant deviations (p < 0.01) from an X2 normality test (1:2:1), fixed in the cross founding parents, and heterozygous in the female F1s. The association of song period with genomic variants was mea- sured using a linear mixed model (LMM) in Gemma (v. 0.98.4; Zhou & Stephens, 2012). For this analysis, missing genotypes were imputed using LinkImpute (v. 1.1.5; Money et al., 2015). The effect size of the variants was estimated as β from the LMM and significance was assessed by likelihood ratio test (alpha = 0.01) after a Bonferroni cor- rection for multiple tests. A classical approach to quantitative trait mapping was also completed using the scanone function of the pack- age qtl (v. 1.5; Broman et al., 2003) in R software (v. 4.1.1; R Develop- ment Core Team, 2022). The logarithm of odds (LOD) score was used to measure association with a significance threshold of 4.08 deter- mined by a genome wide permutation test (alpha = 0.01). Song phenotype genomic architecture The minimum number of genetic components underlying the volley period phenotype was estimated from the variance in the phenotype data (Lande, 1981). Additionally, hyperparameters related to the geno- mic architecture of the trait were estimated by a Bayesian sparse lin- ear mixed model in Gemma (Zhou & Stephens, 2012). Mean estimates were calculated from five independent runs; each run was five million sampling steps and the first 500,000 were discarded as burn-in. Linkage disequilibrium Linkage between variants was measured by Plink (v. 1.9; Chang et al., 2015). Variants were binned into 10 Mb distance bins and the mean R2 was plotted across each autosome for all bins with at least 5,000 SNP comparisons. Genome scan DNA was extracted from wild-caught C. adamsi and C. plorabunda using a Qiagen DNeasy blood and tissue kit. A reduced representation sequencing library was prepared by the Center for Genome Innova- tion at the University of Connecticut with the SbfI cutting enzyme and sequenced on an Illumina HiSeq at the University of Connecticut. Alignment and variant identification were performed as described above for the F2 cross data set, except without filtering based upon X2 normality test or cross-parent genotypes. Weir and Cockerham’s windowed FST was calculated by VCFtools (v. 0.1.6; Danecek et al., 2011) across 50 kb windows with a 5 kb step. Significance threshold was set as z transformed FST >3. No windows rose above the more conservative threshold of zFST >6. Candidate gene identification Proteins from the reference genome assembly for Chrysoperla carnea (GCF_905475395.1) were functionally annotated by EnTAP (v. 0.10.8; Hart et al., 2020) based upon similarity to proteins in the UniProt data- base and gene family assignment by eggNOG. We identified genes anno- tated with the gene ontology term “male courtship behaviour, veined wing generated song production” (GO:0045433) as candidate genes. RESULTS Song phenotype genomic architecture Volley periods for hybrid F1 and F2 offspring were intermediate to the parental phenotypes, showing greater phenotypic range in the F2 off- spring (Figure 1). The average volley period phenotype for the C. adamsi and C. plorabunda parents that founded the cross were 3.39 GENOMICS OF SPECIES-SPECIFIC LACEWING SONG 81 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense and 1.17 s respectively. The average volley period was 1.76 s (SD = 0.13) for F1 and 1.80 s (SD = 0.29) for F2 offspring. The mini- mum number of genetic factors estimated from variance in the pheno- type data to underlying volley period phenotype was 2.5. Estimation of the genetic architecture by trait hyperparameter estimates from the Gemma Bayesian sparse linear mixed model (BSLMM) showed that our complete set of genetic variants explained about 70% of phe- notypic variance (PVE = 0.699 [0.506, 0.888]), while approximately 66 major effect variants explained about 40% of the variance in phe- notype (N = 66.191 [0, 269], PGE = 0.403 [0.000, 0.968]). The top 1% of SNPs according to sparse effect size from our BSLMM were found across all six chromosomes. Genotype-phenotype association After filtering, described in the methods, 9,733 variant loci remained with an average site sequencing depth of �26�. Using two QTL mapping approaches we determined that genetic markers spanning chromosomes one and two were significantly associated with volley period phenotype. In our traditional QTL mapping approach using R qtl, we found that markers across much of chro- mosomes one and two and a small window of chromosome six rose above the genome-wide significance threshold (Figure 2a). Our LMM-based approach with Gemma yielded very similar results: the markers on chromosomes one and two had the largest estimated effect size and were the only markers significantly associated with song phenotype. The phenotypes of F2 offspring are plotted by genotype at the top SNP by the LOD on chromosomes one and two, showing that for both chromosomes, F2 offspring homozygous for the allele from the C. plorabunda parent had the shortest volley period, heterozygotes had an intermediate volley period, and those homozygous for the allele from the C. adamsi parent had the lon- gest volley period (Figure 3). We focus on the QTL on chromo- somes one and two, as they were identified as the regions of the largest effect and were recovered as significant in both QTL map- ping approaches. F I GU R E 2 The genomic basis of volley period song phenotype. Variants are plotted across the six chromosomes with points in alternating light and dark grey. (a) QTL identified from an F2 cross between Chrysoperla plorabunda and Chrysoperla adamsi by the association of volley period with genomic markers in F2 offspring (n = 78). Association is measured by LOD and a red line indicates the significance threshold determined by a permutation test. (b) Divergence between wild-caught C. plorabunda (n = 6) and C. adamsi (n = 6) measured by Weir and Cockerham’s windowed FST calculated across 50 kb windows with a 5 kb step. The mean FST for all windows with >5 variants is plotted and the zFST significance threshold is indicated by a dashed red line. (c) Positions of genes associated with wing-generated song production. Genes are labelled when identity was determined; unlabelled genes were not identified. Candidate genes within song QTL and in close proximity to QTL peaks are nona (no-on-transient A), para (paralytic), per (period), and dsx (doublesex). F I G U R E 3 Volley period of F2 hybrid songs from a Chrysoperla plorabunda and Chrysoperla adamsi cross. Song phenotypes are plotted for F2 offspring of each possible genotype, AA homozygous for the C. plorabunda grandparent genotype, AB heterozygous, and BB homozygous for the C. adamsi grandparent genotype at the top variant by LOD on (a) chromosome one and (b) chromosome two. Individual phenotypes are plotted as points over boxplots for each group. 82 TAYLOR ET AL. 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense Linkage disequilibrium (LD) Estimates of LD across the genome in the F2 cross progeny revealed fewer recombination events on chromosomes one and two, which contain major QTL for song phenotype (Figure 4). We binned SNPs by physical distance in the reference genome and found that distant SNPs show higher estimated R2 on chromosomes one and two com- pared to the other three autosomes. Genome scan Using a sample of six individuals from across the ranges of each of these two species, significant divergence peaks, identified by FST, were found on chromosomes 1, 2, 3, and 6 that rose above the zFST threshold of three (Figure 2b). We identified the FST peaks found on chromosomes 1 and 2 as being within the large QTL identified for the volley period, suggesting that divergence between species in these regions of the genome may be related to song phenotype. Candidate genes Thirty-two proteins from the Chrysoperla carnea genome were anno- tated with the gene ontology term “male courtship behavior, veined wing generated song production” (GO:0045433). The number of pro- teins annotated with GO:0045433 from chromosome one to six, respectively, were 5, 19, 1, 3, 0, 4. The majority of those proteins were mapped to chromosomes one and two, which contained large QTL for the volley period (Figure 2c). Gene identity was determined by EggNOG and when unidentified was assigned the identity of the top UniProt similarity search hit to an arthropod. Correspondence between QTL, FST peaks, and the positions of candidate genes was assessed. The outliers on chromosome one between 18.6 and 26.0 Mb were more than 1 Mb from all identified candidate gene; the two closest candidates were slowpoke (0.8 Mb) and moesin (39.9 Mb). Several significant FST outlier windows fall between 81 and 94 Mb on chromosome one. This region contains the gene nonA (87 Mb), as well as other non-candidate genes. On chromo- some two, multiple FST outlier windows were found between 11.3 and 40.1 Mb; this region included the candidate genes fruitless (11.6 Mb), paralytic (20.3 Mb), period (23.9 Mb), doublesex (27.9 Mb), one unidentified candidate gene (36.8 Mb), and other non-candidate genes. In all cases, there were non-candidate genes closer to signifi- cant FST peaks than the candidate genes. Both candidate genes and non-candidate genes may be involved in song phenotype. DISCUSSION We explored the genomic basis of reproductive isolation and specia- tion in a rapidly evolving lacewing group. Using association mapping for an F2 cross between Chrysoperla plorabunda and Chrysoperla adamsi, we identified two large genomic regions associated with mat- ing song phenotype on chromosomes one and two. The majority of genes associated with Drosophila vibrational courtship signals are found within the song-associated regions for Chrysoperla, suggesting that the genomic basis of these traits may be shared across multiple insect orders. Though these results align with our expectations, we cannot rule out the role of non-candidate genes in song-associated regions of the genome. Using small samples of wild-caught lacewings, we identified sig- nals of divergence that may be related to selection or a lack of intro- gression near our trait of interest. As genome-wide signals of divergence between species may also be related to selection for other traits or other demographic processes, we identified divergence sig- nals most likely to be related to song phenotype by narrowing our search to only include regions with elevated divergence in our volley period QTL. Using this targeted approach we identified several song- associated regions with elevated FST containing candidate genes. Some of these candidates are the period gene, plus several genes known to modify the abdominal vibration in Drosophila: fruitless, dou- blesex, and nonA, which is responsible for species–specific differences in vibrational signals in Drosophila (Campesan et al., 2001; Fabre et al., 2012; Gleason, 2005). Though these candidate genes may be responsible for differences in volley period between these species, other non-candidate genes in these regions may be what truly under- lies song phenotype. Functional validation will be necessary to con- firm the relationship between genes and song phenotype and to identify how they might modify song. RNAi for these candidates would be a clear next step, as it is already being successfully used on closely related insects in the genus Chrysopa (Zhang et al., 2022). We identified two major effect regions in our association of geno- type and phenotype, but it is likely that other loci of small effect related to volley period are spread throughout the genome and are undetected in this analysis. Our estimation of trait architecture identi- fied that only about 40% of phenotypic variance can be explained by F I GU R E 4 Linkage across the five Chrysoperla autosomes in the F2 cross progeny with line type indicating chromosome. The R2 for SNPs is plotted against the physical distance between those SNPs in the reference genome. Chromosomes 1 and 2, which both contain large QLT for volley period phenotype, are coloured red. GENOMICS OF SPECIES-SPECIFIC LACEWING SONG 83 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense large effect loci, while up to 70% of the variance can be explained by the total panel of genomic variants, though we note that there were wide credible intervals around these estimates. It is unsurprising that more genomic regions are involved in the phenotype than those iden- tified by quantitative trait mapping, as we expect a bias towards iden- tifying the loci of large effect with this type of approach (Broman, 2001). The effect size estimates for the QTL we are able to detect could be inflated by these undetected regions of small effect (Beavis, 1998). Additionally, because we mapped a polygenic trait using a single cross, it is possible that we have failed to capture other song-associated variants present in the population but not fixed between species and the parents founding this cross (Muranty, 1996). Despite these limitations, we can say with confidence that the regions identified on chromosomes one and two likely have the largest effect on volley period differences between C. plorabunda and C. adamsi. Future analysis with replicated crosses, larger numbers of offspring, or a greater number of generations will be necessary to identify the pre- cise genomic basis of mating song phenotypes in these species and other carnea-group species. We also observed lower rates of recombination on chromosomes 1 and 2 compared to the other autosomes in our laboratory cross. One potential explanation for decreased recombination in this context could be chromosomal rearrangements between species. We do not expect that selection on song phenotype could explain decreased recombination, as we did not select for this in our laboratory experi- ment. Potentially, selection for other traits related to survival in our laboratory setting in the same areas of the genome could explain decreased recombination, or other features of the genomic landscape such as of repetitive content or gene density (Stapley et al., 2017). Whatever the source, the decreased rates of recombination observed here could function to keep together multiple genomic loci important to song phenotype or song preference, especially in the presence of gene flow. Overall, we have made significant steps towards characteriz- ing the genetic architecture of song volley period phenotype, a mating song feature critical to reproductive isolation in the Chry- soperla carnea-group. We identified major effect loci on chromo- somes one and two containing candidate genes and found evidence of decreased recombination in song-associated regions of the genome. AUTHOR CONTRIBUTIONS EJW, MMW, and CSH conceived of the project and secured funding, KLT, MMW and CSH conducted the experiments, KLT analysed the data and wrote the manuscript, all authors edited the manuscript. ACKNOWLEDGMENTS We would like to acknowledge the valuable input by Dr. Jill Wegrzyn on experimental design and analysis approaches, and thoughtful feed- back on this project from Drs. Elizabeth Jockusch, Chris Simon, David Wagner, DeWayne Shoemaker and Megan Fritz. Computational resources were provided by the University of Connecticut Computa- tional Biology Core. FUNDING INFORMATION University of Connecticut Research Excellence Grant #4626370 to CS Henry, MM Wells, and EJ Wade. CONFLICT OF INTEREST The authors have no conflicts to disclose. DATA AVAILABILITY STATEMENT The data that support the findings of this study are openly available in NCBI at https://www.ncbi.nlm.nih.gov/bioproject/, reference number PRJNA880641. REFERENCES Beavis, W.D. (1998) QTL analyses: power, precision, and accuracy. In: Molecular dissection of complex traits. CRC Press, Boca Raton, USA. pp. 145-162. Blankers, T., Oh, K.P. & Shaw, K.L. (2019) Parallel genomic architecture underlies repeated sexual signal divergence in Hawaiian Laupala crickets. Proceedings of the Royal Society B: Biological Sciences, 286(1912), 20191479. https://doi.org/10.1098/rspb.2019.1479 Broman, K.W. (2001) Review of statistical methods for QTL mapping in experimental crosses. Lab Animal, 30(7), 44–52. Broman, K.W., Wu, H., Sen, S. & Churchill, G.A. (2003) R/qtl: QTL mapping in experimental crosses. Bioinformatics, 19(7), 889–890. https://doi. org/10.1093/bioinformatics/btg112 Campesan, S., Dubrova, Y., Hall, J.C. & Kyriacou, C.P. (2001) The nonA gene in Drosophila conveys species-specific behavioral characteris- tics. Genetics, 158(4), 1535–1543. Chang, C.C., Chow, C.C., Tellier, L.C., Vattikuti, S., Purcell, S.M. & Lee, J.J. (2015) Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience, 4(1), s13742-015. https://doi.org/ 10.1186/s13742-015-0047-84 Coyne, J.A. (1992) Genetics and speciation. Nature, 355(6360), 511–515. https://doi.org/10.1038/355511a0 Crowley, L.M. (2021). The genome sequence of the common green lace- wing, Chrysoperla carnea (Stephens, 1836). Wellcome Open Research, 6(334), 334. https://doi.org/10.12688/wellcomeopenres.17455.1 Danecek, P., Auton, A., Abecasis, G., Albers, C.A., Banks, E., DePristo, M.A. et al. (2011) The variant call format and VCFtools. Bioinformatics, 27(15), 2156–2158. https://doi.org/10.1093/bioinformatics/btr330 Danecek, P., Bonfield, J.K., Liddle, J., Marshall, J., Ohan, V., Pollard, M.O. et al. (2021) Twelve years of SAMtools and BCFtools. GigaScience, 10(2), giab008. https://doi.org/10.1093/gigascience/giab008 Fabre, C.C.G., Hedwig, B., Conduit, G., Lawrence, P.A., Goodwin, S.F. & Casal, J. (2012) Substrate-borne vibratory communication during courtship in Drosophila melanogaster. Current Biology, 22(22), 2180– 2185. https://doi.org/10.1016/j.cub.2012.09.042 Gleason, J.M. (2005) Mutations and natural genetic variation in the court- ship song of Drosophila. Behavior Genetics, 35(3), 265–277. https:// doi.org/10.1007/s10519-005-3219-y Greenspan, R.J. & Ferveur, J.-F. (2000) Courtship in drosophila. Annual Review of Genetics, 34(1), 205–232. https://doi.org/10.1146/ annurev.genet.34.1.205 Hart, A.J., Ginzburg, S., Xu, M.(.S.)., Fisher, C.R., Rahmatpour, N., Mitton, J.B. et al. (2020) EnTAP: bringing faster and smarter functional annotation to non-model eukaryotic transcriptomes. Molecular Ecology Resources, 20(2), 591–604. https://doi.org/10.1111/1755-0998.13106 Henry, C.S., Brooks, S.J., Duelli, P., Johnson, J.B., Wells, M.M. & Mochizuki, A. (2013) Obligatory duetting behaviour in the Chryso- perla carnea-group of cryptic species (Neuroptera: Chrysopidae): its role in shaping evolutionary history. Biological Reviews, 88(4), 787– 808. https://doi.org/10.1111/brv.12027 84 TAYLOR ET AL. 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://www.ncbi.nlm.nih.gov/bioproject/ https://doi.org/10.1098/rspb.2019.1479 https://doi.org/10.1093/bioinformatics/btg112 https://doi.org/10.1093/bioinformatics/btg112 https://doi.org/10.1186/s13742-015-0047-84 https://doi.org/10.1186/s13742-015-0047-84 https://doi.org/10.1038/355511a0 https://doi.org/10.12688/wellcomeopenres.17455.1 https://doi.org/10.1093/bioinformatics/btr330 https://doi.org/10.1093/gigascience/giab008 https://doi.org/10.1016/j.cub.2012.09.042 https://doi.org/10.1007/s10519-005-3219-y https://doi.org/10.1007/s10519-005-3219-y https://doi.org/10.1146/annurev.genet.34.1.205 https://doi.org/10.1146/annurev.genet.34.1.205 https://doi.org/10.1111/1755-0998.13106 https://doi.org/10.1111/brv.12027 Henry, C.S., Wells, M.L.M. & Holsinger, K.E. (2002) The inheritance of mat- ing songs in two cryptic, sibling lacewing species (Neuroptera: Chry- sopidae: Chrysoperla). In: Genetics of Mate Choice: From Sexual Selection to Sexual Isolation. Springer, Dordrecht, NL, pp. 269–289. Henry, C.S., Wells, M.M. & Pupedis, R.J. (1993) Hidden taxonomic diver- sity within Chrysoperla plorabunda (Neuroptera: Chrysopidae): two new species based on courtship songs. Annals of the Entomological Society of America, 86(1), 1–13. https://doi.org/10.1093/aesa/86.1.1 Kronforst, M.R., Young, L.G., Kapan, D.D., McNeely, C., O’Neill, R.J. & Gilbert, L.E. (2006) Linkage of butterfly mate preference and wing color preference cue at the genomic location of wingless. Proceedings of the National Academy of Sciences, 103(17), 6575–6580. https:// doi.org/10.1073/pnas.0509685103 Lande, R. (1981) The minimum number of genes contributing to quantita- tive variation between and within populations. Genetics, 99(3–4), 541–553. Li, H. & Durbin, R. (2009) Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics, 25(14), 1754–1760. https://doi.org/10.1093/bioinformatics/btp324 Mackay, T.F.C. (2001) The genetic architecture of quantitative traits. Annual Review of Genetics, 35(1), 303–339. https://doi.org/10.1146/ annurev.genet.35.102401.090633 Mazzoni, V., Anfora, G. & Virant-Doberlet, M. (2013) Substrate vibrations during courtship in three Drosophila species. PLoS One, 8(11), e80708. https://doi.org/10.1371/journal.pone.0080708 Money, D., Gardner, K., Migicovsky, Z., Schwaninger, H., Zhong, G.-Y. & Myles, S. (2015) LinkImpute: fast and accurate genotype imputation for nonmodel organisms. G3 GenesjGenomesjGenetics, 5(11), 2383– 2390. https://doi.org/10.1534/g3.115.021667 Mullen, S.P. & Shaw, K.L. (2014) Insect speciation rules: unifying concepts in speciation research. Annual Review of Entomology, 59(1), 339–361. https://doi.org/10.1146/annurev-ento-120710-100621 Muranty, H. (1996) Power of tests for quantitative trait loci detection using full-sib families in different schemes. Heredity, 76(2), 156–165. https://doi.org/10.1038/hdy.1996.23 Noor, M.A.F., Grams, K.L., Bertucci, L.A. & Reiland, J. (2001) Chromosomal inversions and the reproductive isolation of species. Proceedings of the National Academy of Sciences, 98(21), 12084–12088. https://doi. org/10.1073/pnas.221274498 Nosil, P. & Schluter, D. (2011) The genes underlying the process of specia- tion. Trends in Ecology & Evolution, 26(4), 160–167. https://doi.org/ 10.1016/j.tree.2011.01.001 R Development Core Team. (2022) R: a language and environment for statis- tical computing. R Foundation for Statistical Computing, Vienna, AT. http://www.R-project.org Rochette, N.C., Rivera-Col�on, A.G. & Catchen, J.M. (2019) Stacks 2: analyt- ical methods for paired-end sequencing improve RADseq-based pop- ulation genomics. Molecular Ecology, 28(21), 4737–4754. https://doi. org/10.1111/mec.15253 Schwander, T., Libbrecht, R. & Keller, L. (2014) Supergenes and complex phenotypes. Current Biology, 24(7), R288–R294. https://doi.org/10. 1016/j.cub.2014.01.056 Seehausen, O., Butlin, R.K., Keller, I., Wagner, C.E., Boughman, J.W., Hohenlohe, P.A. et al. (2014) Genomics and the origin of species.Nature Reviews Genetics, 15(3), 176–192. https://doi.org/10.1038/nrg3644 Stapley, J., Feulner, P.G.D., Johnston, S.E., Santure, A.W. & Smadja, C.M. (2017) Variation in recombination frequency and distribution across eukaryotes: patterns and processes. Philosophical Transactions of the Royal Society B: Biological Sciences, 372(1736), 20160455. https:// doi.org/10.1098/rstb.2016.0455 Templeton, A.R. (1981) Mechanisms of speciation—a populations genetic approach. Annual Review of Ecology and Systematics, 12(1), 23–48. https://doi.org/10.1146/annurev.es.12.110181.000323 Wilkins, M.R., Seddon, N. & Safran, R.J. (2013) Evolutionary divergence in acoustic signals: causes and consequences. Trends in Ecology & Evolu- tion, 28(3), 156–166. https://doi.org/10.1016/j.tree.2012.10.002 Wu, C.I. & Ting, C.T. (2004) Genes and speciation. Nature Reviews Genetics, 5(2), 114–122. https://doi.org/10.1038/nrg1269 Xu, M. & Shaw, K.L. (2021) Extensive linkage and genetic coupling of song and preference loci underlying rapid speciation in Laupala crickets. Journal of Heredity, 112(2), 204–213. https://doi.org/10.1093/jhered/esab001 Zhang, T., Liu, X., Zhang, L., Wang, M., Li, Y. & Mao, J. (2022) Four insulin- like peptides orchestrate reproductive signaling of the green lace- wing, Chrysopa pallens (Rambur) (Neuroptera: Chrysopidae). Annals of the Entomological Society of America, 115(4), 352–359. https://doi. org/10.1093/aesa/saac007 Zhou, X. & Stephens, M. (2012) Genome-wide efficient mixed-model anal- ysis for association studies. Nature Genetics, 44(7), 821–824. https:// doi.org/10.1038/ng.2310 How to cite this article: Taylor, K.L., Wade, E.J., Wells, M.M. & Henry, C.S. (2023) Genomic regions underlying the species-specific mating songs of green lacewings. Insect Molecular Biology, 32(2), 79–85. Available from: https://doi. org/10.1111/imb.12815 GENOMICS OF SPECIES-SPECIFIC LACEWING SONG 85 13652583, 2023, 2, D ow nloaded from https://resjournals.onlinelibrary.w iley.com /doi/10.1111/im b.12815 by U niversity O f M aryland, W iley O nline L ibrary on [02/10/2023]. See the T erm s and C onditions (https://onlinelibrary.w iley.com /term s-and-conditions) on W iley O nline L ibrary for rules of use; O A articles are governed by the applicable C reative C om m ons L icense https://doi.org/10.1093/aesa/86.1.1 https://doi.org/10.1073/pnas.0509685103 https://doi.org/10.1073/pnas.0509685103 https://doi.org/10.1093/bioinformatics/btp324 https://doi.org/10.1146/annurev.genet.35.102401.090633 https://doi.org/10.1146/annurev.genet.35.102401.090633 https://doi.org/10.1371/journal.pone.0080708 https://doi.org/10.1534/g3.115.021667 https://doi.org/10.1146/annurev-ento-120710-100621 https://doi.org/10.1038/hdy.1996.23 https://doi.org/10.1073/pnas.221274498 https://doi.org/10.1073/pnas.221274498 https://doi.org/10.1016/j.tree.2011.01.001 https://doi.org/10.1016/j.tree.2011.01.001 http://www.r-project.org https://doi.org/10.1111/mec.15253 https://doi.org/10.1111/mec.15253 https://doi.org/10.1016/j.cub.2014.01.056 https://doi.org/10.1016/j.cub.2014.01.056 https://doi.org/10.1038/nrg3644 https://doi.org/10.1098/rstb.2016.0455 https://doi.org/10.1098/rstb.2016.0455 https://doi.org/10.1146/annurev.es.12.110181.000323 https://doi.org/10.1016/j.tree.2012.10.002 https://doi.org/10.1038/nrg1269 https://doi.org/10.1093/jhered/esab001 https://doi.org/10.1093/aesa/saac007 https://doi.org/10.1093/aesa/saac007 https://doi.org/10.1038/ng.2310 https://doi.org/10.1038/ng.2310 https://doi.org/10.1111/imb.12815 https://doi.org/10.1111/imb.12815 Genomic regions underlying the species-specific mating songs of green lacewings INTRODUCTION METHODS Collecting, laboratory rearing, phenotype recording Family cross-sequencing and bioinformatic analysis Song phenotype genomic architecture Linkage disequilibrium Genome scan Candidate gene identification RESULTS Song phenotype genomic architecture Genotype-phenotype association Linkage disequilibrium (LD) Genome scan Candidate genes DISCUSSION AUTHOR CONTRIBUTIONS ACKNOWLEDGMENTS FUNDING INFORMATION CONFLICT OF INTEREST DATA AVAILABILITY STATEMENT REFERENCES