The reasons why some individuals are solitary and others gregarious are the subject of ongoing debate as we seek to understand the emergence of sociality. Recent studies suggest that the expression of aggregation behaviors may be linked to the gut microbiota of the host. Here, we tested this hypothesis in females of the European earwig. This insect is ideal for addressing this question, as adults both naturally vary in the degree to which they live in groups and show inter-individual variation in their gut microbial communities. We video-tracked 320 field-sampled females to quantify their natural variation in aggregation and then tested whether the most and least gregarious females had different gut microbiota. We also compared the general activity, boldness, body size, and body condition of these females and examined the association between each of these traits and the gut microbiota. Contrary to our predictions, we found no difference in the gut microbiota between the most and least gregarious females. There was also no difference in activity, boldness, and body condition between these two types of females. Independent of aggregation, gut microbiota was overall associated with female body condition, but not with any of our other measurements. Overall, these results demonstrate that a host's gut microbiota is not necessarily a major driver or a consequence of aggregation behavior in species with inter-individual variation in group living and call for future studies to investigate the determinants and role of gut microbiota in earwigs.

Animal sampling and laboratory rearing

We field sampled 1000 F. auricularia females and males in a peach orchard near Valence, France (Lat 44.9772790, Long 4.9286990) at the end of June 2022. All these individuals were sampled as mated adults and belonged to F. auricularia species “A” [1, 2]. Upon collection, the 1000 individuals were first mixed in a single plastic container and then randomly distributed among four other plastic containers of 100 females and 100 males each (called “test containers”) and two additional plastic containers of 100 females each (called “attraction containers”). The test container held females for which we then measured the aggregation level, general activity, boldness, fresh weight, and gut microbial diversity (see below). The attraction containers held females that served as stimuli for measuring aggregation levels. All plastic containers were grounded with moistened sand and cardboard shelters. Throughout the experiment, the animals were fed with carrot chunks ad libitum, which were changed twice a week. Containers were maintained under a controlled 12:12 light-dark cycle at 20°C and 18°C, respectively [3]. All measurements proceeded three days after laboratory conditions had been set up.

Aggregation measurements

To investigate the association between the level of aggregation of a female and her gut microbial diversity, we first measured the aggregation level of 320 females from the “test containers” (Fig S1). To this end, we used 3D-printed arenas consisting of four linearly aligned circular chambers (diameter 4 cm) (Fig S1; Van Meyel et al., 2022). Three of the chambers were connected by 0.5 cm wide corridors allowing earwigs to move between chambers. The fourth chamber (called the “attraction chamber”) had a reduced corridor of 0.15 cm to prevent animals from reaching the other chamber while allowing the circulation of odors and antennal contacts between individuals on both sides [5].

The measurement of each aggregation score started by placing a randomly selected female from the two test containers in the circular chambers at one end of the system, and three females taken randomly from the two attraction containers in the attraction chamber at the other end of the system (we alternated the orientation of the system between trials). All these females had access to the food source until they were used in the 3D-printed arenas. To avoid injuries during animal handling, all these individuals were previously anaesthetized with CO₂ for 30 s. One hour later (to allow recovery after anesthesia), we started to record whether the focal female was in the chamber next to the group of females (yes or no) and then repeated this measurement by taking pictures every hour for 72 hours using infrared cameras and the software pylon Viewer v5.1.0 (Basler©, Ahrensburg, Germany). For each female, we thus calculated an aggregation score, which was defined as the total number of pictures in which a female was in the chamber next to the group of females for 72 hours. No female was used as both a focal and a non-focal individual. We measured the aggregation score of 66 females per week for five consecutive weeks (10 died during the experiment). Note that the aggregation score (AS) calculated over 24h was consistent over the three days of the experiment (Linear models; ASday2~ASday1: F_1,310 = 175,52, Adjusted R² = 0.358, P < 0.001; ASday3~ASday1: F_1,310 = 52.85, P < 0.001, Adjusted R² = 0.132; and ASday3~ASday2: F_1,310 = 206.29, P < 0.001, Adjusted R² = 0. 387) and did not change over the five weeks of the experiment (ASday1: F_4,310 = 1.25, P = 0.292; ASday2: F_4,310 = 0.08, P = 0.990; ASday3: F_4,310 = 0.25, P = 0.912; all interactions P > 0.141) (Fig S2A).

At the end of the 72 hours, we isolated each tested female in a Petri dish (5 cm diameter) lined with moistened sand for further measurements (see below), while we returned the three attractor females to the attraction containers. To avoid cross-contamination with biological material between trials (e.g., due to chemical signatures or feces), we grounded each plastic arena with a paper sheet that was renewed at each trial and cleaned the glass plates with ethanol between each trial.

For the rest of the experiment, we selected the 20 females (out of 320) with the lowest aggregation scores (AS from 0 to 15, n =20; later called low-aggregation females) and the 20 females with the highest aggregation scores (AS from 49 to 72, n =20; later called high-aggregation females) calculated over 72 hours. These females were evenly distributed over the five weeks of the experiment, i.e. 4 females per week and category (Fig S2B).

Behavioral and morphological measurements

We then measured the general activity and boldness of the 40 high- or low-aggregation females (Fig S1). At the end of the aggregation test (i.e. in the early afternoon), we gently transferred each female to a circular arena (diameter 7.8 cm) held between two glass sheets and maintained on an infrared light table. We then video-recorded females’ activity for 15 minutes in the dark and under infrared light (BASLER BCA 1300, Germany; Media Recorder v4.0, Noldus Information Systems, Netherlands) – as this species is nocturnal and lucifugous. We defined the level of general activity as the total distance walked by a female during these 15 minutes [6]. This distance was automatically extracted from our videos using the video-tracking software EthoVision XT 16 (Noldus Information Technology©, Wageningen, Netherlands). We measured boldness using the same videos and the same software by extracting the time spent by each female in the central area (diameter 5.8 cm) of the tested arena. This measurement of boldness is standard in rodents [7, 8] and applies to earwigs, as they are highly thigmotactic and typically avoid open areas [9]. At the end of the behavioral measurements, we analyzed the 72 pictures sampled per female to select the four and four females (of the given session/week) with the least and highest aggregation scores, respectively.

The day after (morning), we measured the body condition of the 40 females by first anaesthetizing each of them with CO₂ and then measuring their eye distance and fresh weight [10]. Using these two measurements, we calculated the initial body condition for each female based on the “scaled mass index” [11]. In brief, this index standardizes body mass at a fixed value of a linear body measurement based on the scaling relationship between these measures [12]. Accordingly, this index indicates which mass a particular female would have at the average eye distance. We measured eye distance to the nearest 0.01 mm using a binocular scope coupled to the Leica Application Suite software (Leica Biosystems©, Wetzlar, Germany). We weighed each female using an isoCaL Quintix® precision balance to the nearest 0.01 mg (Sartorius©, Göttinger, Germany). Note that the 8 females from the last week were dissected before weight measurement and were thus not included in this measurement.

DNA extraction and amplification protocol

Immediately after morphological measurement (i.e., less than 24 hours after the last aggregation score recording), we extracted the gut of the 40 selected females to investigate whether high or low aggregation females had different gut microbial communities. The remaining females (and males) were used for other experiments not shown in this study. Following the protocol detailed in Van Meyel et al. (2021), each female was CO₂ anaesthetized and dissected under sterile conditions. We immediately flash-freeze each gut with liquid nitrogen and stored them individually in Eppendorf tubes at -80°C until DNA extraction. Total genomic DNA was extracted using the NucleoMag® Tissue extraction kit (Macherey-Nagel™, Düren, Germany) and the V3-V4 region of the 16S rRNA gene was amplified with the prokaryotic primers 343F (5′- ACGGRAGGCAGCAG – 3′) and 784R (5′- TACCAGGGTATCTAATC – 3′) [14] using the Taq Polymerase Phusion® High-Fidelity PCR Master Mix with GC buffer (Qiagen, Hilder, Germany). PCR amplification steps involved an initial denaturing step for 30 sec at 98°C, followed by 22 amplification cycles (denaturation for 10 sec at 98°C; annealing for 30 sec at 61.5°C; extension for 30 sec at 72°C), and ended by a final extension step of 10 min at 72°C. We checked the amplification success and probe specificity by electrophoresis on a 1.5% agarose gel. We amplified all samples in duplicates and equally pooled them for a final product of 30 µL. Each step involved several blank samples. We verified the inter-sample contamination by testing the material used during the dissection and extraction steps (i.e., dissection tools washed in T1 buffer, sample tubes empty with T1 buffer). These samples were placed in the same extraction plate of the other samples and followed the same PCR protocol. Finally, blanks with molecular water were also used during the amplification steps. As all blanks were negative, they were not sequenced (Figure S3). We sequenced each amplicon by 2x250bp Illumina MiSeq technology at the Bio-Environnement platform (University of Perpignan, France).

Bioinformatic processing and diversity metrics

The obtained sequence reads were trimmed and filtered using the quality profiles from the DADA2 algorithm v1.24.0 [15], cleaned for errors, dereplicated, and inferred towards Amplicon Sequence Variants (ASVs) [16]. Chimaeras were removed and taxonomy assignment was performed on a count table where forward and reverse ASVs were merged and pooled, using the SILVA reference database (release 138) [17]. The multiple alignment of the sequences was provided with the MAFFT program [18] and we inferred the phylogenetic tree with the FastTree 2 tool [19] and Phangorn package v2.8.1 [20]. The table was then transformed into a phyloseq object using the phyloseq package v1.40.0 [21] of which we removed the organites and eukarya sequences. Finally, we normalized all samples at 10 180 sequences (i.e., the minimum sample sum of the dataset) and obtained a final dataset constituted by 397 020 sequences (1 593 ASVs) for 20 high- and 19 low-aggregation females. For each sample, we calculated taxonomic and phylogenetic alpha diversity with the Shannon and the Hill1 indices (q = 1) [22]. We then measured the taxonomic and phylogenetic dissimilarities computed on relative abundances of ASVs between pairs of gut microbiota using Bray-Curtis and Weighted Unifrac distances, respectively [23]. The community homogeneity (i.e. microbial dispersion) within low- and high-aggregation status is calculated for each beta diversity metric used with the function betadisper from the package vegan v2.6.2 [24].

Statistical analyses

We conducted all statistical analyses using R v4.2.0 [25] loaded with car v3.1 [26], lme4 v1.1.29 [27], veganv2.6.2 [28] and ape v5.6.2 [29]. We ran six independent statistical models where the response variable was either Shannon (Linear Model, LM), Hill1 (LM), Bray-Curtis dispersion (LM), Weighted Unifrac dispersion (LM), Bray-Curtis dissimilarity matrix (Permutational Multivariate Analyses of Variances; PERMANOVA) or Weighted Unifrac dissimilarity matrix (PERMANOVA). In each of these models, we entered the same 5 explanatory variables, namely aggregation status (high or low), general activity, boldness, body condition (i.e., scaled body mass index), and sampling week. We first included the interaction between sampling week and each variable in the different models and then, where appropriate, removed it after model simplification by AIC comparison. The PERMANOVA on the beta-diversity dissimilarity matrices included 999 permutations and was run using the adonis2 function from the package vegan v2.6.2 [28]. We performed post-hoc pairwise Adonis tests between sampling weeks when it was significant in the PERMANOVA. For each statistical model, we checked that all assumptions were fulfilled using the DHARMa R package v0.4.6 [24]. There was no correlation between the continuous variables entered in each statistical model (Pearson tests; P > 0.05 and |R| ≤ 0.28 for all pairs).

To identify potential microbial discriminants of females with high or low-aggregation status, we used a Linear discriminant analysis with Effect Size (LEfSe, [30]) based on microbial Genera conducted on the Galaxy platform.

Finally, to test whether high- and low-aggregation females differ in their other life-history traits, we performed four t-tests on the general activity, boldness, and body condition (i.e., scaled body mass index).

References

1. Wirth T, Le Guellec R, Vancassel M, Veuille M (1998) Molecular and reproductive characterization of sibling species in the European earwig (Forficula auricularia). Evolution 52:260–265. https://doi.org/10.1111/j.1558-5646.1998.tb05160.x
2. González-Miguéns R, Muñoz-Nozal E, Jiménez-Ruiz Y, et al (2020) Speciation patterns in the Forficula auricularia species complex: cryptic and not so cryptic taxa across the western Palaearctic region. Zoological Journal of the Linnean Society 190:788–823. https://doi.org/10.1093/zoolinnean/zlaa070
3. Meunier J, Wong JWY, Gómez Y, et al (2012) One clutch or two clutches? Fitness correlates of coexisting alternative female life-histories in the European earwig. Evol Ecol 26:669–682. https://doi.org/10.1007/s10682-011-9510-x
4. Van Meyel S, Devers S, Meunier J (2022) Earwig mothers consume the feces of their juveniles during family life. Insect Science 29:595–602. https://doi.org/10.1111/1744-7917.12941
5. Van Meyel S, Meunier J (2022) Costs and benefits of isolation from siblings during family life in adult earwigs. Animal Behaviour 193:91–99. https://doi.org/10.1016/j.anbehav.2022.09.003
6. Merleau L-A, Larrigaldie I, Bousquet O, et al (2022) Exposure to pyriproxyfen (juvenile hormone agonist) does not alter maternal care and reproduction in the European earwig. Environ Sci Pollut Res 29:72729–72746. https://doi.org/10.1007/s11356-022-20970-z
7. Archer J (1973) Tests for emotionality in rats and mice: A review. Animal Behaviour 21:205–235. https://doi.org/10.1016/S0003-3472(73)80065-X
8. Stratton JA, Nolte MJ, Payseur BA (2021) Evolution of boldness and exploratory behavior in giant mice from Gough Island. Behav Ecol Sociobiol 75:65. https://doi.org/10.1007/s00265-021-03003-6
9. Rankin SM, Palmer JO (2009) Dermaptera. In: Encyclopedia of Insects. Elsevier, pp 259–261
10. Thesing J, Kramer J, Koch LK, Meunier J (2015) Short-term benefits, but transgenerational costs of maternal loss in an insect with facultative maternal care. Proc R Soc B 282:20151617. https://doi.org/10.1098/rspb.2015.1617
11. Kramer J, Meunier J (2016) Maternal condition determines offspring behavior toward family members in the European earwig. BEHECO 27:494–500. https://doi.org/10.1093/beheco/arv181
12. Peig J, Green AJ (2010) The paradigm of body condition: a critical reappraisal of current methods based on mass and length: The paradigm of body condition. Functional Ecology 24:1323–1332. https://doi.org/10.1111/j.1365-2435.2010.01751.x
13. Van Meyel S, Devers S, Dupont S, Dedeine F (2021) Alteration of gut microbiota with a broad-spectrum antibiotic does not impair maternal care in the European earwig. Peer Community in Evolutionary Biology
14. Muyzer G, de Waal EC, Uitterlinden AG (1993) Profiling of complex microbial populations by denaturing gradient gel electrophoresis analysis of polymerase chain reaction-amplified genes coding for 16S rRNA. Appl Environ Microbiol 59:695–700. https://doi.org/10.1128/aem.59.3.695-700.1993
15. Callahan BJ, McMurdie PJ, Rosen MJ, et al (2016) DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13:581–583. https://doi.org/10.1038/nmeth.3869
16. Glassman SI, Martiny JBH (2018) Broadscale ecological patterns are robust to use of Exact Sequence Variants versus Operational Taxonomic Units. mSphere 3:e00148-18. https://doi.org/10.1128/mSphere.00148-18
17. Quast C, Pruesse E, Yilmaz P, et al (2012) The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research 41:D590–D596. https://doi.org/10.1093/nar/gks1219
18. Katoh K (2002) MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Research 30:3059–3066. https://doi.org/10.1093/nar/gkf436
19. Price MN, Dehal PS, Arkin AP (2010) FastTree 2 – Approximately maximum-likelihood trees for large alignments. PLoS ONE 5:e9490. https://doi.org/10.1371/journal.pone.0009490
20. Schliep KP (2011) phangorn: phylogenetic analysis in R. Bioinformatics 27:592–593. https://doi.org/10.1093/bioinformatics/btq706
21. McMurdie PJ, Holmes S (2013) phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8:e61217. https://doi.org/10.1371/journal.pone.0061217
22. Chao A, Chiu C-H, Jost L (2010) Phylogenetic diversity measures based on Hill numbers. Phil Trans R Soc B 365:3599–3609. https://doi.org/10.1098/rstb.2010.0272
23. Lozupone C, Knight R (2005) UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol 71:8228–8235. https://doi.org/10.1128/AEM.71.12.8228-8235.2005
24. Hartig F (2022) DHARMa: Residual diagnostics for hierarchical (multi-Level / mixed) regression models. R package version 0.4.6
25. R Core Team (2022) R: A language and environment for statistical computing
26. Fox J, Weisberg S, Fox J (2019) An R companion to applied regression, 3rd ed. Sage, Thousand Oaks, Calif
27. Bates D, Mächler M, Bolker B, Walker S (2015) Fitting linear mixed-effects models using lme4. J Stat Soft 67:. https://doi.org/10.18637/jss.v067.i01
28. Oksanen J, Simpson G, Blanchet F, et al (2022) vegan: community ecology package
29. Paradis E, Schliep K (2019) ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35:526–528
30. Segata N, Izard J, Waldron L, et al (2011) Metagenomic biomarker discovery and explanation. Genome Biol 12:R60. https://doi.org/10.1186/gb-2011-12-6-r60