Urbanization reduces gut bacterial microbiome diversity in a specialist ground beetle, Carabus convexus

Urbanization is rapidly shaping and transforming natural environments, creating networks of modified land types. These urbanization‐driven modifications lead to local extinctions of several species, but the surviving ones also face numerous novel selection pressures, including exposure to pollutants, habitat alteration, and shifts in food availability and diversity. Based on the assumption that the environmental pool of microorganisms is reduced in urban habitats due to habitat alteration, biodiversity loss, and pollution, we hypothesized that the diversity of bacterial microbiome in digestive tracts of arthropods would be lower in urban than rural habitats. Investigating the gut bacterial communities of a specialist ground beetle, Carabus convexus, in forested rural versus urban habitats by next generation high‐throughput sequencing of the bacterial 16S rRNA gene, we identified 3839 bacterial amplicon sequence variants. The composition of gut bacterial samples did not significantly differ by habitat (rural vs. urban), sex (female vs. male), sampling date (early vs. late spring), or their interaction. The microbiome diversity (evaluated by the Rényi diversity function), however, was higher in rural than urban adults. Our findings demonstrate that urbanization significantly reduced the diversity of the gut bacterial microbiome in C. convexus.

recalcitrant plant polymers, the production of essential vitamins and nutrients (Engel & Moran, 2013), the degradation of toxins absorbed from the environment or ingested with food items (Kikuchi et al., 2012), regulating the immune system and cellular homeostasis of the host (Bai et al., 2021;Douglas, 2015), as well as enhancing the host resistance to pathogens and/or parasites (Douglas, 2015;Engel & Moran, 2013).
Gut bacterial microbiome diversity in wild insects can be linked to the biodiversity of the habitats where the host organisms live, since the gut microbiomes are likely acquired from food items and other sources in the habitat (Borer et al., 2013;Engel & Moran, 2013).Thus, changes in habitat parameters, especially the impoverishment of food supply and biodiversity would influence the composition of the gut microbial communities (Wang et al., 2011).
The exposure to various xenobiotics can also reduce the diversity of gut microbiome (Motta et al., 2018;Yang et al., 2018).While effects of human disturbances associated with agriculture on insect gut microbiome have been studied (Giglio et al., 2021;Magagnoli et al., 2022;Motta et al., 2018), the impact of urbanization on the same was addressed only on flying insects (Bosmans et al., 2018;Nobles & Jackson, 2020).Studies on dispersal limited, terrestrial insects are completely absent.To fill this gap, we studied the impact of urbanization-driven environmental disturbances on the gut microbiome diversity of Carabus convexus Fabricius, 1775 (Coleoptera: Carabidae), a flightless habitat specialist ground beetle.We selected a ground beetle because it is a representative of a speciose insect family whose taxonomy, behaviour and ecology are well studied (Lövei & Sunderland, 1996), including their responses to urbanization (Magura & Lövei, 2021).
Specifically, we hypothesized that (1) the diversity of gut bacterial microbiome would be lower in urban beetles compared to their rural conspecifics, due to the reduced food spectrum, the general impoverishment of biodiversity, and the presence of xenobiotics in urban habitats.Further, (2) we hypothesized sex-specific differences in gut bacterial communities, as males of ground beetles are more mobile than females, especially during the reproductive period, therefore may acquire a more diverse set of facultative symbionts from their environment.
Our results supported the first but not the second hypothesis.
Urban beetles had a less diverse gut microbial community than their rural conspecifics, but we found no sex-or season-specific differences.

| Study area and sampling design
The study area was on the eastern part of the Hungarian Plain, in and near the city of Debrecen (47°32′ N; 21°38′ E, 202,402 inhabitants and an area of 461.7 km 2 ).In the lowland forest adjacent to the city, four rural stands were selected.All forest stands were part of a protected area (the Great Forest of Debrecen), and part of the Natura 2000 network (site code: HUHN20033).This lowland forest is recently embedded within a matrix composed of undisturbed or moderately disturbed meadows, grasslands, pastures, and farmland.All selected sites (3.71, 3.75, 3.94, and 3.04 ha) were mature lowland forest stands (>120 years) dominated by English oak (Quercus robur).The formerly continuous lowland forest adjacent to the city historically extended into the northern part of the city, but today only fragmented, isolated patches remain.From these patches (also dominated by English oak) four fragments (3.21, 3.86, 3.98, and 3.33 ha) were selected.Rural and urban sites were at least 250 m from each other (mean distances: 396.5 m between rural and 702.2 m between urban sites).The classification of rural and urban sites was based on the proportion of built-up areas within a 1000 m radius circle around the given site.In the rural forest, there was no built-up area, while in the urban area >60% of the surface was built-up or was drastically different from the rural habitats.In the rural forest stands, there was no regular forestry intervention, with few visitors, and minimal trampling.
Contrary to this, in the urban forest fragments the maintenance operations were regular (smaller fallen branches and trunks were removed, although larger ones were cut into smaller pieces and left on the ground), and the shrub layer was strongly thinned.
These fragments were also isolated by asphalt-or gravel-covered paths, and had numerous visitors, and substantial trampling.
Furthermore, urban sites were characterized by environmental changes associated with urbanization, expressed in higher soil and soil surface temperature, higher soil pH values, lower amount of decaying wood and higher concentration of Ca and Zn in the soil (Bogyó et al., 2015).The minimum distance between the rural and urban areas was 3 km.

| Model organism and sampling
The chosen ground beetle, Carabus convexus is a widespread Eurasian species (Turin et al., 2003).In Central Europe, this nocturnal, predatory species with extraintestinal digestion reproduces in early spring, and becomes active in late March -early April.Teneral individuals appear in late July and adults go to overwinter in November.
This medium large (14-23 mm), wingless (brachypterous) species has limited dispersal power (Turin et al., 2003).In the studied region (Great Hungarian Plain) C. convexus is a forest specialist (Magura, Tóthmérész, et al., 2008); although also occurs in urban forest fragments, where its abundance is significantly lower than in rural forest stands (Magura, Lövei, et al., 2008) due to being very sensitive to urbanization-driven disturbance (Martinson & Raupp, 2013;Niemelä et al., 2002).Adults of C. convexus were collected during the spring breeding period, from the end of March to mid-May in 2021, using 15 live, unbaited pitfall traps at each site (2 areas × 4 sites × 15 traps = 120 traps in total).Traps were placed randomly, at least 10 m apart from each other, and at least 50 m from the nearest forest edge in order to avoid edge effects (Molnár et al., 2001).Traps were square plastic containers (170 mm long × 110 mm wide × 105 mm deep) with shredded leaves to provide shelter, thus reducing both intra-and interspecific predation.Traps were covered with fibreboard roofs (200 mm × 200 mm) to protect them from litter and rain, and to prevent bycatch.Traps were checked twice per week; trapped beetles were transported to the laboratory in cool boxes, identified to species level, sexed, and put individually in micro centrifuge tubes (2 mL) in 96% ethanol and stored at +5°C until dissection.

| Gut dissection, DNA extraction, PCR amplification and Illumina MiSeq analysis
Beetles (15 rural females, 16 rural males, 8 urban females, and 9 urban males) were dissected within 3 days after sampling to minimize changes in the gut microbial community (Kudo et al., 2019).Beetles were surface sterilized by immersing them twice in 96% ethanol for 1 min each, followed by washing in sterile Ringer's solution.After removing the legs and elytra, individuals were pinned and rewashed with sterile Ringer's solution.The tergum plates (tergites) were cut lengthwise on both sides and the cuticle was opened, taking care not to damage the intestinal tract.The foregut and midgut were extracted, carefully washed with ethanol (96%), placed in sterile 1.5 mL micro centrifuge tubes filled with 96% ethanol and stored at − 80°C.
The dissections were performed under a stereo microscope (8-80× magnification, Delta IPOS-810, Poland).After every dissection, the immersing solution (96% ethanol) was replaced and the plate of the microscope, as well as the equipment used for dissection (scalpel, forceps, dissecting needle) were cleaned and disinfected with 96% ethanol.
Microbial genomic DNA was extracted from the gut samples using the QIAamp PowerFecal Pro DNA Kit (QIAGEN GmbH, Hilden, Germany), according to the manufacturers' instructions.Before extraction, the gut samples were washed twice with ice cold phosphate buffer (PBS 280-315 mOsm/kg, pH 7.4).To avoid contamination, we followed the procedure proposed by Eisenhofer et al. (2019).Gut samples from different sites were processed in a random order to avoid any systematic bias and possible batch effects.Negative controls containing sterile water without gut material were included during the DNA extraction.The quantity and quality of the purified DNAs (concentration and A260/280 ratio) were checked with a DeNovix spectrophotometer (DeNovix, Inc., Wilmington, DE, USA) according to the manufacturer's instructions.
The amplification conditions were as follows: an initial denaturation step at 95°C for 3 min, followed by 25 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, extension at 72°C for 30 s, with a final extension step of 5 min at 72°C.To remove free primers and primer dimers, PCR products were subsequently purified using the Agencourt AMpure XP PCR purification system (Beckman Coulter, Brea, CA, United States) as described in the manufacturer's protocol.To be able to multiplex the libraries, Dual Illumina indices (The Nextera XT Index Kit, Illumina, Inc., San Diego, CA, USA) were attached to the PCR products by another PCR which was performed in a 50 μL volume containing 5 μL of the purified amplicon PCR product, 25 μL KAPA HiFi HotStart ReadyMix, 5 μL of each index primer, and 10 μL of PCR grade water.The temperature profile of the amplification was as follows: an initial denaturation step at 95°C for 3 min, followed by eight cycles of 30 s at 95°C, 30 s at 55°C, 30 s at 72°C, and finally an extension step of 5 min at 72°C.Again, a negative control was included (PCR amplification control) by replacing template DNA with sterile water.The amplified fragments were subjected to another purification step using the Agencourt AMpure XP PCR purification system.The final libraries were measured with a Bioanalyzer 2100 (Agilent Technologies, Inc., Santa Clara, CA, USA) using an Agilent DNA 1000 Kit (Agilent Technologies, Inc., Santa Clara, CA, USA) according to manufacturer's protocol to verify the size and validate the quality.The normalized, pooled and denatured library samples were sequenced on the Illumina MiSeq platform (Illumina, Inc., San Diego, CA, United States) using 300 bp paired-end V3 chemistry with the MiSeq Reagent Kit v3 (600 cycle), producing on average 300,000 raw reads per sample and ∼15 Gb of sequence.
Results obtained for the negative controls (DNA extraction control and PCR amplification control) were satisfactory (without microbial reads), confirming that the experimental conditions were met to obtain robust, reliable data.Sequencing and library preparation were performed at the UD-GenoMed Medical Genomic Technologies Ltd., University of Debrecen, Hungary.

| Bioinformatic processing and statistical analysis
The processing of raw MiSeq forward and reverse paired-end reads was based on the DADA2 pipeline (Callahan et al., 2016).The quality profiles of the forward and reverse reads were inspected by a heat map of the frequency of each quality score at each base position, as well as the mean quality score at each position and the quartiles of the quality score distribution.After visualizing the quality profiles, the sequence ends were truncated and filtered using the maximum number of expected errors allowed in a read (a better filter than simply averaging quality scores, see Edgar & Flyvbjerg, 2015).The core sample inference algorithm was applied to the filtered and trimmed sequence data (Callahan et al., 2016).Subsequently, the forward and reverse reads were merged, an amplicon sequence variant (ASV) table was constructed, and chimeras (4.4% of the total reads) were removed.
Finally, the taxonomic assignment of ASVs was performed using the SILVA 16S database for bacterial sequences (Quast et al., 2013).
All statistical analyses were performed using the R program environment (version 4.2.1;R Core Team, 2022).Bioinformatic processing was done by using the dada2 package (Callahan et al., 2016).The processed bioinformatics dataset was subjected to median normalization (Dillies et al., 2012) to deal with differences in sequencing depth (Evans et al., 2018).The composition of gut bacterial microbiome from rural and urban beetles was analysed by non-metric multidimensional scaling ordination using the Bray-Curtis index of dissimilarity.Convex hulls around the points of the groups (habitat: rural vs. urban; sex: female vs. male; season: early vs. late spring) on the ordination plot were also visualized, based on the standard errors of the centroids (average) of points multiplied with the 95% confidence value from the Chi-squared distribution (with 2df).
Calculations of ordination and convex hulls were made with the vegan (Oksanen et al., 2022) and the MASS (Venables & Ripley, 2002) packages.The difference between the group (habitat: rural vs. urban; sex: female vs. male; season: early vs. late spring) centroids in the ordination plot was tested by the envfit function (by setting the number of permutations to 999) in the vegan package (Oksanen et al., 2022).
When two variables (habitat or sex or season) were used to create groups, we used the RVAideMemoire package (Hervé, 2021) to determine the differences between centroids.
The diversity of the gut bacterial microbiome was evaluated using the Rényi diversity (entropy) function (Henderson & Southwood, 2016;Maurer & McGill, 2011;Rényi, 1961), calculated by the vegan package (Oksanen et al., 2022).The Rényi diversity, HR(α) is defined as: where p i is the relative frequency of the i-th species, S the total number of species and α is the scale parameter (α ≥ 0, α ≠ 1).
At α = 0, the value of the Rényi diversity is the logarithm of the number of species of the community.At small values of the scale parameter (α), the Rényi diversity value is influenced by the rare species, while as the scale parameter increases, the diversity value is more and more influenced by the common species.Thus, this scalable diversity approach gives a diversity profile instead of a single number, allowing a more sophisticated assessment of diversity (Maurer & McGill, 2011).
We investigated the effects of urbanization level (rural vs. urban habitat) and sex (female vs. male) on the diversity of gut bacterial microbiome using linear mixed models (LMM) with the lme4 package (Bates et al., 2015).We tested the best-fitting probability distribution for the Rényi diversity values using the car (Fox & Weisberg, 2019) and MASS packages (Venables & Ripley, 2002).Based on the results, the Rényi diversity values were modelled with a normal error distribution (Zuur et al., 2009).Habitat type (rural vs. urban), sex and their interaction were fixed effects.The nested design (sampling sites nested within sampling areas) was also considered in the models.When LMM revealed a significant difference between means, the Tukey test was performed for multiple comparison among means (Zuur et al., 2009).

| Sampled beetles
From the end of March to mid-May 2021, we collected 106 C. convexus individuals in the four rural and four urban habitats.The mandibles of all beetles were little worn, indicating that they belonged to the overwintered cohort in their first breeding season.Eighty-nine beetles (49 females and 40 males) were collected at the rural, while 17 individuals (8 females and 9 males) from the urban sites.

| Sequencing data
We sequenced the hypervariable V3 and V4 regions of the 16S rRNA gene from 48 gut microbial community samples of 15 rural females and 16 males, 8 urban females, and 9 males.The rural samples were randomly selected from the total, so that the gut microbiome samples were roughly comparable.Our dataset of 48 samples yielded 3839 different ASVs with a total read count of 300,888 (for ASVs: mean = 126.52,SD = 64.72).From the 31 rural gut samples, 2893 ASVs with 229,006 reads (for ASVs: mean = 142.26,SD = 65.18), while from the 17 urban samples 1222 ASVs with 71,882 reads (for ASVs: mean = 97.82,SD = 54.62) were obtained.

| Phylogenetic structure of the gut bacterial microbiome
A total of 11 bacterial phyla were recorded in C. convexus gut samples.On the basis of their average relative abundance, the most representative phyla in both habitat types were: a de novo phylum (average 44.90% in all individuals, while average 44.83% in rural and 45.01% in urban beetles), Firmicutes (average 24.66% in all individuals, 22.97% in rural and 27.74% in urban beetles), and Proteobacteria (average 23.48% in all beetles, 25.43% in rural and 19.93% in urban individuals), amounting to 93.04% (93.24% in rural, while 92.68% in urban) of the bacterial community.
A total of 59 bacterial families were identified from the C. convexus gut samples.The most abundant families were a de novo family (average relative abundance in all samples: 62.61%), Mitochondria (9.69%), and Lactobacillaceae (5.28%).Thirty-eight families were present in both habitat types, while 14 families were recorded only in beetles from rural and seven only in urban habitats.
The identified taxa corresponded to 119 bacterial genera.The most abundant of them were: a de novo genus (average 77.42% in all individuals), Lactobacillus (average 4.77% in all individuals), Carnobacterium (1.79%), Megasphaera (1.74%), Clostridium sensu stricto 1 (1.51%), and Vagococcus (1.11%).Five genera were present in at least 50% of the samples (present in ≥24 individuals), so they can be considered as the obligate members of the gut bacterial microbiome.In contrast, 75 genera were recorded ≤10% of the samples (<5 beetles), indicating that they were likely to be facultative members.We observed little difference in the average abundance of the common genera between the rural and urban, as well as females versus males.However, the Carnobacterium genus was very abundant in a single rural female, while Vagococcus was present in a single urban female (Figure 1).

| Diversity of the gut bacterial microbiome
The Rényi diversity of the gut bacterial microbiome in rural individuals was significantly higher than in urban ones (except for a marginal (p = .0593)significance at α = 0), indicating that apart from the number of ASVs, the bacterial gut microbiome was unequivocally more diverse in rural beetles (Table 1, Figure 3).For all C. convexus individuals combined, LMM models revealed no significant sexspecific difference in microbial diversity (Table 1).Similarly, the habitat × sex interaction was not a significant predictor (Table 1), indicating that the urbanization level (rural vs. urban) had no sexspecific effect.The diversity of the gut bacterial microbiome was not significantly different between female and male beetles in either habitat (Figure S5).

| Gut bacterial microbiome of C. convexus
The 16S rRNA sequencing data revealed a mean of 142.26 and 97.82 ASVs present per rural versus urban beetles, respectively.This number is roughly in line with previous next generation sequencing studies based on ASVs, where 98.6 ASVs per beetle was detected in the gut of three ground beetle species from the subfamilies Harpalinae and Brachininae (Silver et al., 2021).
Our data revealed the presence of 11 bacterial phyla in the alimentary tract, comparable to the 19 gut-associated bacterial phyla found in 218 insect species in 21 taxonomic orders (Yun et al., 2014).Firmicutes and Proteobacteria were the predominant phyla, representing more than 40% of the gut bacterial community These phyla are commonly associated with the insect's gut (Colman et al., 2012;Yun et al., 2014), including ground beetles (Do et al., 2022;Giglio et al., 2021;McManus et al., 2018;Silver et al., 2021).The main families in the gut of C. convexus were Abundance per gut sample of the bacterial genera from common (>1% average abundance in all samples) families in rural and urban Carabus convexus females and males.
The only previous study focusing on Carabus species (Carabus arboreus Lewis, 1882 and Carabus albrechti Morawitz, 1862) in Japan, however, found that Enterococcaceae were the most abundant (Kudo et al., 2019), emphasizing that host taxonomy may considerably influence the composition of gut bacterial community (Colman et al., 2012).

| Measuring the diversity of the gut bacterial microbiome
Gut microbiome may contain numerous ASVs of low abundance even after the removal of chimeras and noises (Engel & Moran, 2013).
Commonly used diversity indices (the number of OTUs/ASVs, Shannon diversity) are considerably influenced by the presence of rare species, and provide a partially distorted assessment (Maurer & McGill, 2011) but this disadvantage is overcome by the use of scalable diversity comparisons (Lövei, 2005).Potential candidates for scalable diversity comparisons include the Hill numbers because they are parameterized by a parameter value, q, which determines the sensitivity of the measures to the relative abundances of taxa (Hill, 1973).Recently, a few studies have used this approach for the analysis of the gut bacterial microbiome (e.g.Brunetti et al., 2022;Littleford-Colquhoun et al., 2019).Another powerful scalable diversity method, incorporating several classical diversity measures and measuring diversity across the continuous abundance-spectrum of taxa, is based on the Rényi diversity (entropy) function (Maurer & McGill, 2011).However, to the best of our knowledge, this method has not yet been used to characterize and compare the diversity of the gut bacterial microbiomes.Using this scalable diversity approach for the first time to compare the diversity of gut bacterial communities, we demonstrated its outstanding performance and strongly recommend its use in further studies.

| Effects of urbanization on gut bacterial microbiome of C. convexus
The impact of urbanization on ground beetles is manifested at all levels of biological organization, from the individual (population) level, through the community and ecosystem levels (Magura & Lövei, 2021).At the ecosystem level, changes in biotic interactions are of great importance and the microbiomes can profoundly modulate these (Douglas, 2022).However, our knowledge on the role of microbiomes in biotic interactions in ground beetles is currently too fragmentary.
Our study is the first to evaluate the impact of urbanization on Microorganisms in the host environment represent a potential pool of the gut microbiome (Borer et al., 2013;Engel & Moran, 2013;Kudo et al., 2019).As urbanization is usually associated with biodiversity loss, the decreased microbiome diversity in the urban C. convexus beetles may be due to the reduced diversity of the soil and/or soil-surface bacteria in the urban environment.However, evidence for an impoverished microbial diversity in urban soils is controversial.Wang et al. (2017) found that the soil microbial diversity was significantly higher in rural than urban or suburban soils in Southeastern China, while another study (Gao et al., 2023) reported the opposite.Soil bacterial diversity was significantly the lowest in the rural and highest in urban forest soils along a Finnish urbanization gradient (Scholier et al., 2023).Another comprehensive study, Our urban sites were characterized by higher soil and soil surface temperatures, higher soil pH values, less decaying wood and higher concentration of some metals than in comparable rural sites (Bogyó et al., 2015;Simon et al., 2016).Furthermore, our urban forest patches were considerably fragmented and even isolated by impervious surfaces (Magura, Tóthmérész, et al., 2008).It seems that under these conditions, soil bacterial diversity is reduced.
the gut bacterial microbiome in a member of the Carabinae subfamily.Our results show that urbanization significantly decreases the diversity of the gut bacterial microbiome in this species.Similar phenomena were found in adult house sparrows (Passer domesticus (Linnaeus, 1758);Teyssier et al., 2020) and bumblebee queens (Bombus terrestris (Linnaeus, 1758);Bosmans et al., 2018).The composition of gut bacterial microbiome of dragonfly adults did not differ by the levels of urbanization, whereas the same in nymphs was significantly influenced by the urbanization level(Nobles & Jackson, 2020).The above results highlight that the composition and diversity of the gut bacterial microbiome in rural versus urban individuals may differ even in closely related species, or different developmental stages of the same species.Our study was based on adult C. convexus from a well-defined region, and studying other populations, and even larvae would allow us to draw more robust conclusions.

F
I G U R E 2 Non-metric multidimensional scaling ordination using the Bray-Curtis index of dissimilarity of the amplicon sequence variants of the gut bacterial microbiome from Carabus convexus adults at the habitat-level (open green circles represent gut bacterial samples from rural beetles, while open red circles represent samples from urban beetles).Ellipses (convex hulls) were based on the standard errors of the averages (centroids) of points multiplied with the 95% confidence value from the Chi-squared distribution (number of permutations: 999).Final stress value: 21.3930.investigatingsoil microbial communities from five cities on three continents found no significant difference in the number of bacterial OTUs between undisturbed rural and disturbed urban sites (EppSchmidt et al., 2017).
Thus, as our results suggest, urban C. convexus individuals can acquire fewer bacteria than their rural conspecifics, contributing to the decline in diversity of the gut bacterial microbiome in urban beetles.However, the β-diversity (expressed as the size of the convex hull on the ordination plot) of the gut bacterial communities of urban C. convexus males and females was higher than rural ones, indicating TA B L E 1 Summary of LMM results and post hoc tests on the Rényi diversity at different values of the scale parameter (α) of the gut bacterial microbiome in female and male beetles from rural and urban forested habitats (p values in bold denote significant (p < 0.05) effects, while p values in Italic denote marginally significant (p < 0.1) effects).Results of the Tukey test indicate which habitat differs significantly (p < 0.05) or marginally significantly (p < 0.1) from the other; for example, R > U indicates higher diversity in rural versus urban individuals, while R (>) U indicates marginal significance.

F
The Rényi diversity profiles (mean ± SE) of the amplicon sequence variants of the gut bacterial microbiome of Carabus convexus adults from rural and urban habitats.Signs above the means indicate the results of the Tukey tests (▪ denotes marginal significance (p = 0.0593), while asterisks indicate significant (*: p < 0.05, **: p < 0.01) difference in the diversity values).Note that only the Rényi diversity values at α = 0 to α = 8 are shown.