Population Structure of German Cockroaches (Blattodea: Ectobiidae) in an Urban Environment Based on Single Nucleotide Polymorphisms

Abstract

German cockroaches (Blattella germanica L.) harbor and disperse medically important pathogens and are a source of allergens that impact human health and wellbeing. Management of this pest requires an understanding of their distribution and dispersal. In this study, we collected German cockroaches from three apartment buildings in New Jersey, USA. We identified single-nucleotide polymorphisms (SNPs) from DNA extractions using next generation sequencing. We analyzed the SNPs and characterized cockroach population genetic structure using F_st, principal component, phylogenetic, and STRUCTURE analyses. We found significant differences in German cockroach population structure among the buildings. Within buildings, we found variable population structure that may be evidence for multiple colonization events. This study shows that SNPs derived from next generation sequencing provide a powerful tool for analyzing the genetic population structure of these medically important pests.

As generalists, cockroaches (Order: Blattodea) are able to survive with limited resources (Bell et al. 2007). Some species of cockroaches, such as Blattella germanica (L.), Periplaneta americana (L.) (Blattodea: Blattidae), Periplaneta fuliginosa (Serville) (Blattodea: Blattidae), Periplaneta brunnea (Burmeister) (Blattodea: Blattidae), and Periplaneta australasiae (Fabricius) (Blattodea: Blattidae) (Schal 2011), are pest species that can live in a wide range of habitats including hospitals, restaurants, and dwellings. In urban areas, cockroaches live on a variety of plant- and animal-based foods that are unintentionally provided by human beings (Bell et al. 2007). These pest cockroach species are often widely distributed throughout the world. They are able to disperse independently within buildings but have low motility among buildings without human assisted dispersal (Roth 1985). Together these characteristics lead to strong genetic structuring within and among cockroach populations (Vargo et al. 2014). The German cockroach (B. germanica) is a globally distributed pest in human habitats and is the most common cockroach species in the U.S (Wang et al. 2008). Cockroaches are known vectors of human pathogens (Rivault et al. 1993, Salehzadeh et al. 2007). German cockroaches are known to harbor several medically important pathogens [such as Enterobacter cloacae (Jordan) (Enterobacterales: Enterobacteriaceae), Klebsiella pneumonia (Schroeter) (Enterobacterales: Enterobacteriaceae), Pseudomonas aeruginosa (Schroeter) (Pseudomonadales: Pseudomonadaceae), and Serratia marcescens (Bizio) (Enterobacterales: Yersiniaceae)] both within their guts (Menasria et al. 2014) and also attached to their surfaces (Rivault et al. 1993, Salehzadeh et al. 2007). German cockroaches also produce and spread allergens that contribute to asthma morbidity and affect human health (Gore and Schal 2007).

Chronic pest infestations in residential buildings are a common phenomenon in urban environments (Runstrom and Bennett 1984, Wang et al. 2008, Saenz et al. 2012, Barbu et al. 2014). In multi-unit dwellings, pest infestations can come from three sources: within an apartment when treatments were not able to eliminate a preexisting population, adjacent households and independent pest movement, and new colonizations from the external environment where movement is assisted by human action (Saenz et al. 2012). Evidence shows that the colonization of pests in urban environments often originates from adjacent locations (Runstrom and Bennett 1984, Cooper et al. 2015), indicating independent pest movement. Furthermore, humans and pets can be the vector for pests among more distant locations (Szalanski et al. 2014). Ineffective pest control can lead to resurgence and recolonization, which in turn results in wasted time, effort, and resources (Myers et al. 2000). Cockroach control products and services totaled US$16 billion in 2019 in the USA (Wang et al. 2021). Despite the economic impact of B. germanica, our understanding of their patterns of colonization, control or eradication, and recolonization remains incomplete. Cockroach population genetic structure may be a good tool to understand these dynamics. A well-developed understanding of source-sink dynamics and the resulting population structures of pest species will allow for more effective pest control (Pulliam 1988).

The distribution of the German cockroach is affected by the structure of the urban landscape. Building separation is one of the factors that can shape cockroach dispersal and distribution. For instance, Crissman and colleagues (2010) found distinct cockroach populations among separate buildings within the same city. Another study showed evidence that cockroach distributions may be affected by human transportation between locations (Booth et al. 2011). Physical structure within buildings can also affect cockroach distributions. However, it is difficult to define the population structure of German cockroaches with high resolution among closely related populations. One mark-recapture study indicated that individual cockroaches have active ranges across multiple apartments (Owens and Bennett 1982). The earliest population genetic study of B. germanica in urban environments utilized allozyme loci and did not reveal significant differentiation among populations between two cities 900 km apart in France (Cloarec et al. 1999). More recently, Crissman and colleagues (2010) used nine microsatellite markers to show low yet significantly nonzero F_st values (range 0.014–0.028) within apartment buildings, indicating a high rate of gene flow (Holsinger and Weir 2009), and yet some cockroach genetic structure. However, neither STRUCTURE nor neighbor joining phylogenetic trees were able to define population genetic relationships within or among apartment buildings, but did differentiate among apartment complexes (Crissman et al. 2010). In contrast, a similar microsatellite approach found strong B. germanica population differentiation among swine farms across a broader region in North Carolina (Booth et al. 2011).

Next generation sequencing methods allow for identification and analysis of thousands of single-nucleotide polymorphism (SNP) loci (Brookes 1999, Russello et al. 2015). The effectiveness of this method has been shown in many studies of medically important arthropod species such as Culex (Yurchenko et al. 2020), Anopheles (Neafsey et al. 2010, Clarkson et al. 2020), and ticks (Lado et al. 2020). Therefore, SNPs may allow for a more detailed understanding of cockroach population genetic structure. In this study, we used SNPs to quantify the population genetic structure of B. germanica within and among three apartment buildings in New Jersey, USA. Our aim was to determine the extent of long-distance gene flow between the three apartment buildings (21 to 91 km apart) and whether the population genetic structure of B. germanica supports singular colonization of apartment buildings and local spread, or multiple colonization events of the same apartment building. We hypothesize that the three buildings support distinct cockroach populations with minimal gene flow between them, and that each building supports a single population of cockroaches. Alternatively, buildings may have been colonized multiple times resulting in more than one distinct population in each building. Our use of more than 2000 unlinked SNP loci allows us to examine B. germanica population structure in more detail than previously possible. A detailed understanding of B. germanica population structure is essential for understanding how cockroaches colonize apartment buildings and recolonize individual apartments after local eradication.

Materials and Methods

Study Sites and Cockroach Sampling

German cockroaches were collected from high-rise multi-unit dwellings in the cities of Paterson, Irvington, and Trenton, New Jersey, USA (Fig. 1). Hereafter we refer to the buildings as PB (Patterson Building), IB (Irvington Building), and TB (Trenton Building). The buildings were built between 1953 and 1965, are 11–15 stories tall, and have 112–246 apartment units each (Table 1, see Supp Figs. 1, 2, and 3 [online only] for building layouts). Each of the three buildings received monthly pest control services hired by property management. The pest control technicians primarily relied on the application of cockroach gel baits to control cockroaches. However, many infested apartments were not treated during monthly pest control service if residents were not home during the visit. Based on our observations, the treatment was very brief and very little amount of bait was applied per apartment. Some residents used insecticide sprays on their own to suppress cockroach infestations. Despite these efforts, 20–30% of the apartments in each building were found to be infested with German cockroaches based on a building-wide cockroach survey concurrent to the present study (Table 1).

Cockroaches were collected in April (PB) and May 2017 (IB) and February 2018 (TB) using glue boards (Trapper Monitor & Insect Trap, Bell Laboratories Inc., Madison, WI) placed in apartments for approximately 14 d. Three traps were placed in the kitchen and one trap was placed in the bathroom in each apartment (Wang et al. 2019). Upon collection, traps with cockroaches were sealed in individual plastic bags and brought to the laboratory and stored at −20°C before analysis. Most of the apartments sampled were one bedroom or studio apartments, while less than 10% were two-bedroom apartments.

Our goal was to sample five adult cockroaches from each of six apartments distributed among the floors, for a total of 30 samples from each of the buildings. However, TB and IB did not have six apartments with five intact adult samples each, so in these buildings our sampling was less even (see results, Tables 1 and 2 for details).

DNA Extraction, SNP Detection, and Pruning Linked Loci

Cockroach genomic DNA was extracted from the cockroach bodies (excluding the abdomen) using DNeasy Blood & Tissue Kits (Qiagen, Valencia, CA) in July 2018. DNA samples were sent to SNPsaurus (Oregon, USA) for sequencing and single nucleotide polymorphism (SNP) analysis. NextRAD genotyping-by-sequencing libraries were built by SNPsaurus using the method described by Russello et al. (2015). Genomic DNA was first fragmented with Nextra DNA Flex reagent (Illumina, Inc, San Diego, CA), which also ligates short adapter sequences to the ends of the fragments. The Nextera reaction was scaled for fragmenting 40 ng of genomic DNA, although 50 ng of genomic DNA was used for input to compensate for degraded DNA in the samples and to increase fragment sizes. Fragmented DNA was then amplified for 27 cycles at 74°C, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence GTGTAGAGCC. Thus, only fragments starting with a sequence that can be hybridized by the selective sequence of the primer were efficiently amplified. The nextRAD libraries were sequenced on a HiSeq 4000 (Illumina Inc., San Diego, CA) with one lane of 150 bp reads. The genotyping analysis was done by using custom scripts from SNPsaurus. The vcf file (Danecek et al. 2011) was filtered to remove alleles with a population frequency of less than 3%. After the SNP analysis, SNPsaurus provided a vcf file that included the SNP data of the whole genome of B. germanica and a phylip data file (Felsenstein 1985) based on SNPs.

We used the ‘gdsfmt’ package in R to manage and organize the SNP data (Zheng et al. 2012). We used the ‘snpgdsLDpruning’ function from the ‘SNPRelate’ package in R (Zheng et al. 2012) to prune any loci with a linkage disequilibrium value (Zaykin et al. 2008) greater than 0.2.

Data Analysis

Phylogenetic, principal component (PCA), F_st, and STRUCTURE analyses were conducted as follows to describe the population structure of the cockroach populations. To investigate the phylogenetic relationships within and among the populations from the three apartment buildings, the edge-linked partition model in IQ-TREE was used to infer the maximum likelihood tree from the full SNPs dataset (Nguyen et al. 2015, Chernomor et al. 2016, Trifinopoulos et al. 2016). We used the IQ-TREE web server to obtain branch support with the ultrafast bootstrap (Minh et al. 2013, Hoang et al. 2018) implemented in the IQ-TREE software (Nguyen et al. 2015) with 10,000 bootstrap alignments and 1,000 iterations. Based on BIC criteria, IQ-TREE chose a symmetric model with unequal rates but equal base frequencies (Zharkikh 1994) and a discrete Gamma model (Yang 1994) with 4 rate categories. The final phylogenetic tree was presented via iTOL (Ciccarelli et al. 2006).

PCA was conducted using the ‘SNPRelate’ package (Zheng et al. 2012) in R (R Core Team 2021). We took a hierarchical approach to the PCAs, first analyzing the three buildings together, and then analyzing each building individually. PCA was conducted on the unlinked loci using the function ‘snpgdsPCA’ from the ‘SNPRelate’ package in R (Zheng et al. 2012).

We estimated F_st values and their 95% confidence intervals in a hierarchical manner, first among the three buildings and then within each building using the ‘betas’ function following the method of Weir and Hill (2002) from the ‘hierfstat’ package in R (Goudet 2005, Weir and Goudet 2017).

We used STRUCTURE version 2.3.4 (Pritchard et al. 2000; Falush et al. 2003, 2007; Hubisz et al. 2009), which uses a Bayesian algorithm to assign individuals with shared ancestry into population clusters based on the proportion of their genome that matches the cluster, to further explore cockroach genetic structure. We implemented STRUCTURE with an admixture model using the SNP dataset pruned of linked loci. The STRUCTURE algorithm was first applied to the full pruned dataset including samples from all three buildings, then it was applied to each of the three buildings separately. We used 50,000 burn-in periods and 200,000 repeats. We ran five replicates for each presumed number of populations (K), from K = 1 to K = 8 (and up to K = 14 for IB).

Because our sampling was unevenly distributed among apartments within buildings, we utilized the MedMeaK, MaxMeaK, MedMedK, and MaxMedK methods proposed by Puechmaille (2016) to identify the optimal number of population clusters, which avoids downwards the bias in the number of K that may result from the commonly used method proposed by Evanno and colleagues (2005) in such situations. The Puechmaille method begins by assigning sampling groups (apartments) to clusters based on either the mean or the median of the membership coefficients of the samples within the sampling group, repeated for each level of K tested. Next, the number of clusters with one or more sampling groups assigned is counted, and either the median or maximum of these counts across replicate runs of a given K is recorded. The maximum of these values across all levels of K tested is the optimal value of K identified by the approach (Gilbert 2016, Puechmaille 2016). We plotted the optimal population cluster composition graphs with Clumpak (Kopelman et al. 2015) within StructureSelector (Li and Liu 2018).

Results

Cockroach infestation rates in TB were much lower than that in IB and PB (Table 1). Among infested apartments, median cockroach densities were 2.5 times higher in IB than in TB and PB, and mean cockroach densities were 3–4 times higher in PB and IB compared to TB. Because there were not sufficient apartments with five intact, adult cockroaches in IB and TB, we were not able to sample five cockroaches from each of six apartments in these two buildings. See Table 2 for details.

A total of 5,839 SNPs were discovered, with 2,419 remaining after linked loci were pruned.

Cockroach Genetic Population Structure Among Buildings

Among the three apartment buildings in New Jersey, the global F_st value (0.116, Table 3) was significantly greater than zero, and yet relatively low, indicating slight but significant genetic differentiation among buildings. PCA supported this conclusion, with clear clustering of samples within buildings (Fig. 2A). Phylogenetic analysis also supported this conclusion with each of the buildings forming a distinct clade with high bootstrap support (PB = 99.6%, IB = 100%, TB = 100%, Fig. 3). Branch lengths representing genetic distance revealed that cockroach populations at PB and TB are more closely related, while the population at IB is more genetically distant (Fig. 3). Structure analysis also supported population differentiation among the three apartment buildings, as the Puechmaille method (2016) indicated that the optimal number of population clusters is three (Table 3, Fig. 4, Supp Fig. 4 [online only]), with all samples from IB and PB showing majority ancestry (>50%) for their respective population clusters, and 27 of 30 samples from TB showing majority ancestry for the third population cluster and plurality ancestry for 2 of the remaining 3 samples.

Cockroach Genetic Population Structure Within Apartment Buildings

F_st values within each building were also significantly greater than zero, and yet relatively low, suggesting substantial gene flow within buildings yet with some population structure (Table 3). Within PB, the F_st value of 0.076 was lower than in the other two buildings, indicating greater gene flow within PB. In TB, the F_st value was higher than global F_st (0.116) (Table 3). Within buildings, the PCA results didn’t show a clustering pattern as obvious as that found between buildings, with some apartments clearly clustered but others mixed with samples from other apartments and even other floors (Fig. 3B–D).

The phylogenetic analysis showed three patterns of relationships among samples from various apartments within buildings. Samples from some apartments, such as TB 807, IB 7U, and IB 9W, were clustered in monophyletic clades, with 100% bootstrap support. Other clades grouped some but not all samples from multiple apartments with 100% support, such as samples from IB 10E and 10F. Samples from some apartments, such as IB 4W and TB 610, were spread across multiple clades, often with low bootstrap support (Fig. 3). The Puechmaille method (2016) identified 4 optimal population clusters within PB and IB, and 7 population clusters within TB, from the STRUCTURE analysis (Table 3, Fig. 4, Supp Fig. 5, 6, and Fig. 7 [online only]). STRUCTURE analysis generally agreed with PCA and phylogenetic analyses, with samples from some individual apartments forming distinct clusters, such as TB 807, IB 7U, and IB 9W, while other clusters suggest samples from disparate apartments are closely related, such as IB 10E and 10F. Samples from some apartments showed mixed ancestry from the optimal clusters, again IB 4W and TB 610 are examples of this pattern (Table 3, Fig. 4).

Discussion

In this study, we used next generation sequencing and single nucleotide polymorphisms to investigate the population genetic structure of German cockroaches distributed across buildings in three New Jersey cities. At the regional scale, phylogenetic, PCA, F_st, and STRUCTURE analyses show a clear genetic differentiation among the three buildings in New Jersey. The F_st value (0.116) of German cockroaches among the three buildings (which are in three different cities) is close to the average F_st value (0.106) for German cockroaches estimated from microsatellite markers from 17 cities in the United States (Vargo et al. 2014). Our F_st value among buildings of 0.116 is significantly higher than 0, indicating slight but significant genetic differentiation at these scales. This F_st value is higher than that found by Crissman and colleagues (2010) (0.048 ± 0.04, 95% CI), while lower than that found by Booth and colleagues (2011) (0.171 ± 0.024, 95% CI). The scale of the various studies may explain this difference. In our study, the buildings ranged from 21 to 91 km apart. Crissman and colleagues (2010), sampled cockroaches from buildings within 10 km from one another, while Booth and colleagues (2011) examined cockroaches collected from hog farms across a 150 km range. If distance among buildings affects the genetic differentiation of B. germanica, this may explain the difference among the F_st values in these three studies, as human dispersal of cockroaches is likely to be reduced with greater distance (Booth et al. 2011).

The PCA, phylogenetic, and STRUCTURE analyses all suggest clear genetic differentiation among the three buildings. The bootstrap values of the building clades are nearly 100% and thus support the hypothesis that the gene flow is lower between buildings than within buildings. From the PCA and STRUCTURE results, we can see the same separation of populations among buildings. This result agrees with the conclusion of the previous studies that there is genetic structuring and limited gene flow within apartment buildings and among localities (Crissman et al. 2010, Booth et al. 2011). Even at these different scales, buildings are likely to be a major barrier to B. germanica movement and to limit gene flow among populations. STRUCTURE and phylogenetic analyses of the samples showed that the samples from TB and PB are genetically more related than the samples from IB (Figs. 3 and 4). This result suggests that in our case cockroach genetic distance is not directly correlated to the geographic distance of populations because IB is geographically in the middle of the location between TB and PB. Other research has shown that the distribution of cockroaches can be more related to the traffic between locations than the simple geographic distance between them (Booth et al. 2011). Trenton and Paterson are connected by trains, which are an important means of public transportation. Therefore, gene flow between TB and PB may be more strongly mediated by human movement even though they are further apart.

Within buildings, our results show a variety of relationships among cockroaches collected from individual apartments and among cockroaches from elsewhere in the buildings. At one end of the spectrum, we see samples from individual apartments, such as TB 406, 608, and 807, and IB 7U and 9W, that are tightly clustered in PCA, phylogenetic, and STRUCTURE analyses (Figs. 2–4). At the other end of the spectrum, we see individuals from multiple apartments that share ancestry from multiple STRUCTURE clusters and are intermixed in PCA and phylogenetic results, including many of the samples from PB and IB. In between these extremes, we find that some apartments appear to house relatively unrelated individuals, such as TB 610 and PB 10E and 10 F. Our within-building F_st values (0.076~0.146) show that gene flow is generally happening within each building, but more granular PCA, phylogenetic, and STRUCTURE analyses suggest that gene flow may vary within buildings.

We hypothesize that these varying patterns within and among the buildings result from different colonization histories. IB has the highest overall cockroach density (Table 1), relatively intermixed PCA results (Fig. 2C), relatively poorly resolved phylogenetic clades (Fig. 3), and mixed ancestry among most individual cockroaches in most apartments (Fig. 4). Together these results suggest an ongoing and persistent cockroach infestation with a substantial movement among apartments. However, the populations in apartments 9W and 7U are more distinct (Figs. 2C, 3, and 4C), suggesting that these populations are either persistent and long standing (note that other roaches in the building share some ancestry with these clusters [Fig. 4C]), or are relatively new colonizations.

Similarly, PB also has a higher cockroach population density. Here we see two smaller and distinct populations (see the populations that dominate apartments 10E and 11E in Figs. 2B, 3, and 4B) and two large panmictic populations that include cockroaches from all six of the apartments (Figs. 2B, 3, and 4B). Note that the apartments sampled in PB are all relatively close to one another, representing lines E and F, which are adjacent, and floors 7, 9, 10, and 11. Thus it is not surprising that cockroaches are moving freely among these apartments. However, the samples from apartments 10E and 11E show very little admixture, leading us to hypothesize that these populations are recent immigrants to the building.

TB is somewhat different, with lower overall cockroach densities and also higher variability than PB (Table 1). TB appears to have controlled cockroach infestations relatively well as most apartments in TB have very few or no cockroaches while some apartments have more, resulting in a patchwork pattern of occurrence and a higher overall F_st. In TB, we see that three putative STRUCTURE populations (Figs. 2D, 3, and 4D) occur in each of three apartments (406, 608, and 807, respectively) with very little admixture among these individuals. Other clusters and apartments in TB have more admixture and variability. We hypothesize that the patterns we have observed result from multiple distinct infestations to 406, 608, and 807.

A number of factors lead us to interpret our results with caution. First, the three buildings were sampled 10 mo apart (see methods), and thus this time lag could cause the populations in the buildings to appear to be more distant from one another. However, given that cockroach lifespans range from 3 to 6 mo, the time lag between samples among buildings is unlikely to have strongly affected our results. Second, we had many apartments in TB with very few or even only one adult cockroach available for sampling. Thus in particular the STRUCTURE results for TB should be interpreted with caution.

Compared with the previous approaches to estimate the population structure of B. germanica that used Random Amplified Polymorphic DNA or Microsatellite DNA (Cloarec et al. 1999, Jobet et al. 2000, Crissman et al. 2010) our SNPs provided more polymorphic loci (~2,500), and hence more information and resolution to the analysis. For instance, Crissman and colleagues (2010) used 8 loci with 7.68 mean alleles per locus, and examined samples from 3 apartments in each of 3 low rise apartment buildings as well as intensive sampling within some apartments and lower-level sampling in additional apartment buildings. Crissman and colleagues (2010) found relatively weak genetic structure among apartments within buildings based on their phylogenetic and STRUCTURE analyses, but in contrast found rather clear FCA clusters within buildings and strong support for apartment level clusters according to pairwise G-tests (Crissman et al. 2010). Crissman and colleagues (2010) concluded that their results suggest common ancestry with subsequent differentiation among apartments. These results are rather different from ours, which show three different situations: isolated apartment population clusters, panmictic clusters spread across multiple apartments, and apartments that house individuals from distinct clusters (Figs. 2–4). One key factor that may explain the difference between Crissman et al. (2010) and our study is that cockroach density appears to have been much higher in the Crissman et al. (2010) study, who left their traps set for 1–2 d (whereas ours were set for 14 d) and apparently had sufficient samples on every trap (whereas we had traps with as few as only one adult). In buildings with very high cockroach density such as those investigated by Crissman et al. (2010), it may be that cockroaches are less likely to migrate successfully when adjacent apartments are already fully occupied, whereas in our case cockroach population levels may have been lower, giving migratory cockroaches more success upon arrival in a new apartment.

Others have shown that German cockroach infestations are spatially correlated within a building based on whether the apartments shared a wall, ceiling/floor, or were across hallways (Zha et al. 2019). Mark-recapture studies show that German cockroaches are able to move readily among adjoining apartments (Owens and Bennett 1982, Runstrom and Bennett 1984, 1990). However, as Crissman and colleagues noted (2010), the low recapture rate in those studies makes it hard to draw solid conclusions about the level of the gene flow between apartments.

While our study appears to be the first to use SNPs to investigate cockroach population genetic structure, SNPs have been used to study a number of other medically important pest species such as Culex (Yurchenko et al. 2020), Anopheles (Neafsey et al. 2010, Clarkson et al. 2020), and ticks (Lado et al. 2020).

From a pest management perspective, it is important to understand the movement of cockroaches at both large and small scales. Future research that samples DNA from individuals from all apartments within a building and explicitly considers connectivity among apartments, elevators, mechanical chases, and garbage chutes, and extermination and control history, may shed further light on mechanisms of cockroach migration within buildings.

Acknowledgments

We thank Jessica L. Ware, Simon J. Garnier, Gareth J. Russell, Caroline M. DeVan, Rebecca Panko, and Donovan Odelugo for their assistance or advice. This research was supported by the US National Science Foundation (DEB-1321265), National Institute of Food and Agriculture, U.S. Department of Agriculture, 1019198 through the NJ Agricultural Experiment Station, the United States Department of Housing and Urban Development Healthy Homes Technical Studies grant program (grant NJHHU0039-17), and Chengdu Research Base of Giant Panda Breeding (2021CPB-B16).

References Cited