The “F4” strain of C. lectularius was sampled in London (UK) in 2008 (Thongprem et al., 2020). Isofemale lines were created immediately after sampling and were reared on human blood in Oliver Otti’s lab (Dresden, Germany) until 2021. For our experiments, we used three of these isofemale lines, provided by O. Otti: F4-V1, F4-V6, and F4-V48. These lines were confirmed by PCR to be negative for torix Rickettsia and are positive for BEV-like. We fed bed bugs on the Hemotek^® system with human blood provided by the Etablissement Français du Sang [blood group: O or A with 17 IU.mL Sodium Heparin (BD Vacutainer); temperature: 36.5 ± 0.5°C; membrane: parafilm^®]. Colonies were maintained in round plastic jars containing corrugated cardboard harborage shelters (41 mm × 30 mm), at 24°C, 60% relative humidity, and with a photoperiod of 12 L:12D.

Couples from either F4-V1, F4-V6, and F4-V48 lines were fed weekly on 1 mL of human blood for 30 min. The cohorts were generated as follows: fed couples were isolated in 24-well plates to mate and oviposit for 5 days, then males were removed. After 2 weeks, females were removed, and nymphs were gathered in round plastic jars containing corrugated cardboard harborage shelters. For each cohort and each nymphal stage, unfed (UF) nymphs (1 to 3 nymphs per cohort) were sampled from the first nymphs hatched (#1 in Supplementary Figure S1) and the remaining nymphs were fed. Nymphs that did not feed were discarded. Sampling was completed at each nymphal stage by collecting nymphs one-day post-feeding (1DPF, #2 in Supplementary Figure S1), and five-days post-feeding (5DPF, #3 in Supplementary Figure S1). Molting to the next instar occurred 1 week after the last meal; nymphs that did not molt were discarded, and samplings were performed (UF, 1DPF, 5DPF) as described previously for the following stages (second, third, fourth, and fifth instar). Bed bugs were frozen immediately after collection and kept at −20°C until DNA extraction. A total of 168 nymphs were sampled, originating from 6 cohorts (1 F4-V1, 2 F4-V6, and 3 F4-V48).

Nymphs at the fifth instar were taken from either F4-V1, F4-V6, and F4-V48 lines and were fed on 1 mL of human blood for 30 min. To quantify Wolbachia density over the nymph-to-adult transition, 5DPF fifth instars were sampled for both sexes, the sex being determined at that stage by a combination of approaches described previously in Langer et al. (2020). One week after the fifth instar feeding, residual nymphs were discarded.

To determine the effect of sex, feeding, and aging on Wolbachia density, 1DPF male and female adults were sampled, and populations were divided equally into two round plastic jars. One jar was fed weekly (condition: weekly feeding), while the other one was not fed (condition: starvation after feeding on day 1). Each week, for 4 weeks, males and females of both jars were sampled at the date corresponding to 1DPF for the fed population (Supplementary Figure S2). Bed bugs were frozen immediately after collection and kept at −20°C until DNA extraction. A total of 146 fifth instars and adults were sampled, originating from 6 cohorts (1 F4-V1, 2 F4-V6, and 3 F4-V48).

DNA was extracted using the MACHEREY-NAGEL NucleoSpin^® 96 Tissue kit, following the protocol for genomic DNA from tissue, except for small individuals (i.e., from first to third instars). Small instar samples were placed in 0.2-mL tubes with three 1.5-mm stainless steel beads while large instar samples were placed in 2-mL tubes with one 5-mm diameter bead. All samples were left for 15 min at −80°C before being ground with a Tissue Lyzer (Retsch, Qiagen) for 2 min at 20 Hz. The following steps of the protocol followed the supplier’s recommendations, except for the small individuals. Indeed, for small individuals, the volume of lysis and wash buffers was divided by 3 compared to the recommendations and compared to the samples considered as large (i.e., from the fourth nymphal stage to the adult stage). For all samples, pre-lysis was carried out for 2 h at 56°C, washes were performed under vacuum, while elution was performed by centrifugation at 5600 g for 2 min with 100 μL of elution buffer. Eluted DNA was stored at −20°C until qPCR quantification.

Conventional PCR was used to verify the presence/absence of symbionts. The wCle 16S rRNA gene [INT2F, INT2R; 136 bp (Sakamoto and Rasgon, 2006)], the γ-proteobacteria 16S rRNA gene [BEV_F, BEV_R; 420 bp (Hosokawa et al., 2010)], the torix Rickettsia gltA gene [RiGltA405_F, RiGltA1193_R; 786 bp (Pilgrim et al., 2017)], and the C. lectularius ribosomal protein (RPL18) gene [RPL18F, RPL18R; 137 bp (Fisher et al., 2018)] were PCR-amplified as follows on a subset of unfed fifth instars (n = 10) and on one control sample known to be infected by the three symbionts. A 25-μL reaction containing 2.5 μL 10X DreamTaq^® Green Buffer, 0.5 μL of dNTP (10 μM), 0.1 μL of DreamTaq^® DNA Polymerase (5 U/μL), 0.5 μL of each primer set (10 μM) (primer sequences in Supplementary Table S1), and 2 μL of template DNA (previously diluted to 1/20) was prepared. All PCR reactions were performed in a Bio-Rad C1000 Touch™ thermal cycler with the following program: an initial denaturation step at 95°C for 5 min, followed by 35 cycles of 95°C for 30 s, Tm °C for 30 s and 72°C for 1 min. A final extension step of 72°C for 5 min was included. To visualize the amplicons, 5 μL of the PCR products were electrophoresed at 100 mV for 25 min (TBE 0.5X, 1% agarose, 0.5% BIOTIUM GelRed^® Nucleic Acid Gel Stain).

A real-time quantitative PCR (qPCR) duplex assay targeting the wCle-16S rDNA and the bed bug RPL18 gene was used to obtain relative quantification of wCle in each bed bug. The 10-μL reaction mix contained 5 μL SsoAdvanced^™ Universal Probe Supermix (Bio-Rad), 500 nM of each forward and reverse primer, 300 nM of each probe, and 2 μL of extracted DNA previously diluted to 1/20. We used primers and probes previously described (Supplementary Table S1). The probes were provided by IDT DNA Technologies and included two quenchers, Iowa Black^® at the 3′ end, and another internal quencher called ZEN^™. Each probe also had a specific fluorophore at the 5’ end, FAM for wCle, and HEX for the RPL18 gene of C. lectularius. Real-time qPCR was performed in QuantStudio 6 Flex^™ (Applied Biosystems) with the following program: an initial denaturation step at 95°C for 30 s, followed by 40 cycles of 95°C for 10 s, 60°C for 30 s and 72°C for 30 s. The qPCR assay efficiency (primer efficiency in Supplementary Table S1) was determined and confirmed in every run by constructing a standard curve using serial dilutions of a purified and quantified amplicon of each target. In addition, the specificity and absence of inhibitors in the samples were assessed upstream using dilutions of DNA extracts produced under the same conditions.

Each qPCR measurement was made in duplicate. Values for which the delta of quantification cycles (Cq) between replicates was greater than 0.5 cycles or the values of Cq were higher than 35 cycles were discarded [wCle-16S: 7.14% discarded (n = 12), RPL18: 5.36% discarded (n = 9; 8 in common with wCle-16S)]. Most of these discarded samples were young nymphal stages (71.42% of discarded are first instars).

The relative quantity (RQ) of Wolbachia within the insect was calculated based on equation one of Pfaffl (2001) where the ratio between the quantity of insect gene (reference) and wCle gene (target) relied on the PCR efficiency (E) and the number of Cq:

We used the R software (version 4.1.0) for all analyses (R Core Team, 2024) and the ggplot2 R package (version 3.4.4) for all plots (Wickham, 2016). To test the effect of different factors (i.e., development, feeding, sex) and their interactions on relative wCle quantities, we analyzed log₁₀-transformed data using a linear mixed effect model with the lme4 R package (Bates et al., 2015).

For nymphal (N) samples, we tested the following model: lmer [log₁₀(RQ) ~ Age_N + Feeding_State + Age_N:Feeding_State + (1|Cohort)], where Age_N is a quantitative variable associated to nymphal development (weeks), Feeding_State is a qualitative factor associated to blood feeding status (UF/1DPF/5DPF), and Cohort is a random factor linked to replicate sampling (A to F). Unfed (UF) first instars are considered as references in the model.

For adult (A) samples, because the feeding factor varies only after the adult’s first blood meal, we first tested on fifth instars 5DPF and adults 7DPF of both sexes (n = 26) the effect of development (fifth instar to adult) on the relative quantity in wCle in both sexes, with the following model: lmer [log₁₀(RQ) ~ Age + Sex + Sex:Age + (1|Cohort)], using the fifth instar females as reference. We then tested on adults only the effect of aging, blood feeding, and sex, with the following model: lmer [log₁₀(RQ) ~ Age_A + Sex + Feeding:Age_A + Sex:Age_A + Feeding:Sex:Age_A + (1|Cohort)], using the weekly fed-first week females as reference. Age and Age_A are quantitative variables associated with development and aging, respectively; Sex is a qualitative factor (female, male); Feeding is a qualitative factor associated to blood feeding status (fed, unfed), whose dynamics depends on the aging/duration of starvation; and Cohort is a random factor linked to replicate sampling (G to L).

Normality and homoscedasticity were checked graphically on residuals for each fitted model. Residuals were also checked for homogeneity of variance. Statistics for global effects of factors and their interactions are reported in Supplementary Tables S2, S4, S5. Model coefficients are reported with their mean ± confidence intervals (0.95). Interaction effects can be interpreted as an additive effect compared with the reference. Because calculating p-values from linear mixed effect model is complex and subject to controversy (Wasserstein and Lazar, 2016; Wasserstein et al., 2019), p-values estimated using the lmerTest R package (Kuznetsova et al., 2017) are only specified in Supplementary Tables S2, S4, S5.

We documented the localization of wCle for each nymphal stage after molting (10 nymphs per instar) in the F4-V48 line, using the FISH technique (Thongprem et al., 2020). Tissues were dissected in 1X PBS and preserved immediately in Carnoy solution (chloroform: ethanol: glacial acetic acid = 6:3:1) overnight. All samples were cleared by incubation in 6% H₂O₂ in ethanol until the body was transparent, i.e., for at least 24 h for first and second instars, and up to 5 weeks for third to fifth instars. We then used a tungsten micro-needle to make micropores in the nymph cuticle to allow the fluorescence probes to diffuse through the cuticle during the hybridization step. The samples were hybridized by incubating the tissues overnight in a hybridization buffer (20 mM Tris–HCl pH 8.0, 0.9 M NaCl, 0.01% Sodium dodecyl sulfate 30% formamide) with 10 pmol/mL of rRNA specific probes for wCle (Hosokawa et al., 2010; probe sequences in Supplementary Table S1). After incubation, tissues were washed in buffer (0.3 M NaCl, 0.03 M sodium citrate, 0.01% sodium dodecyl sulfate) and mounted onto a slide using fluoro Gel with DABCO (Electron Microscopy Science) as a mounting medium. Slides were then observed under Leica THUNDER Imager 3D Assay epifluorescence microscope.

Length of bacteriome and thorax were measured on ImageJ 1.53c (Schneider et al., 2012) taking the maximum length (L_b) and width (W_b) of both bacteriomes, and the length (L_p) and width (W_p) at the middle of the pronotum. Relative bacteriome sizes were calculated by the ratio between pronotum estimated area (L_p ×W_p) and the mean estimated area (π × (1/2) L_b × (1/2) W_b) of both bacteriomes. We used the R software (version 4.1.0) for all analyses and plots (R Core Team, 2024).

To test the effect of age on absolute area of the bacteriome, we analyzed data using a Linear Model (LM). We tested the following model: lm (log₁₀(Mean_area_bact) ~ Stage), where Stage is a qualitative factor associated to nymphal development (first to fifth instar) and Mean_area_bact is the mean of the estimated area of both bacteriomes in a nymph.

To test the effect of the age on the relative area of the bacteriome, we analyzed data using the following model: lm (Mean_ratio_area ~ Stage), where Mean_ratio_area is the mean of the estimated area of both bacteriomes relatively to the estimated area of the thorax in the nymph.

Normality and homoscedasticity were checked graphically on residuals for each fitted model. Residuals were checked for homogeneity of variance. LM statistics were given for global effects of factors. To compare the means of either the area or the relative size of the bacteriome between stages, we performed post-hoc Tukey tests (see Supplementary Table S3 for the complete output of each model).

We performed the following manipulations at the Centre Technologique des Microstructures (CTμ, Université Claude Bernard Lyon 1, France). To avoid losing the samples during the procedure, bacteriomes of fifth instar females were dissected in PBS together with the proximal ovary and the surrounding cuticula. Samples were placed into specimen carriers previously covered with 0.5% lecithin in chloroform. Carriers were directly loaded into the HPM100 high pressure freezer and fast frozen. After freezing, samples were placed into a Leica EM AFS II Freeze substitution machine and incubated at -90°C for 30 h in substitution buffer (1% osmium tetroxide, 0.25% uranyl acetate, 0.5% glutaraldehyde, 1.5% H₂O in acetone) and then gradually heated to -30°C (+5°C/h) and maintained at that temperature for 24 h. Samples were washed with -30°C cold acetone and then gradually heated to 20°C (+10°C/h). Samples were placed in 25% Epon^™ (Epoxy substitute embedding medium kit from Sigma-Aldrich®) in acetone for 3 h, 50% Epon^™ in acetone for 17 h, 75% Epon^™ in acetone for 3 h, then in three different baths of pure Epon^™ for a total of 24 h. Infiltration and embedding in 1.7% benzyl dimethyl amine in Epon™ resin were performed during 4 days at 60°C.

Sectioning was performed on Leica EM UC7 ultramicrotome using a diamond knife (Diatome) and mounted on uncoated copper grids. Ultra-thin sections (70 nm) for Transmission Electro-microscopy (TEM) were observed at 120 keV using a JEOL 1400 electron microscope. Electron micrographs of standard sections were taken with GATAN DigitalMicrograph software (Pleasanton, CA) and further analyzed using ImageJ 1.53c.

To analyze the dynamics of the obligatory endosymbiont wCle during nymphal development, we measured by qPCR duplex assay the relative quantity of wCle in 157 nymphs at different stages of development. wCle quantity was below the level of detection in 2.55% of the samples (n = 4), which were all first instars. We tested the effect of development and feeding status on wCle relative quantity on the remaining 153 samples, considering the variability observed between cohorts as a random effect in the model (Supplementary Figure S3).

The log₁₀-transformed relative quantity of wCle linearly increased with age (Figure 1) with a slope equal to 0.68 ± 0.11 (confidence interval at 0.95), equivalent to an increase by 4.68 of the relative quantity of wCle at each nymphal stage; the log₁₀-transformed absolute quantity of wCle, as measured by quantification of its 16S rDNA, increased ~2 times more strongly than that of the bed bug gene RPL18 genes (Supplementary Figure S4). This result confirms an increase in the wCle quantity relative to the number of host cells over nymphal development.

The blood feeding status also impacted the relative wCle quantity, with an increase of 0.73 ± 0.45 on log₁₀-transformed relative quantity (i.e., equivalent to an increase by 5.37 of the relative quantity) between unfed and 5DPF (Days Post Feeding) nymphs. This positive impact of blood feeding on the wCle relative quantity dimmed as the nymph was growing (see interaction Age_N*Feeding State 5DPF Figure 1B).

To determine if the increase in wCle relative quantity is associated with a change in the bacteriome shape or an over-growth of endosymbiont within the bacteriome, we performed Fluorescence In situ Hybridization (FISH) using probes targeting the endosymbiont 16S rRNA. wCle was detected in bacteriomes of each sex at all instars during nymphal development (Figures 2A–F). All the bacteriocytes observed within a bacteriome were infected, and their cytoplasm was packed with wCle (Figures 2A–E; Supplementary Figure S5). TEM observation performed on fifth instar females confirmed that wCle exhibits a three-layer membrane presumably composed of the two membranes of the symbiont and an individual vacuolar membrane of insect origin similar to what has been described 50 years ago (Chang and Musgrave, 1973) (Figure 2G).

To determine if the increase in wCle relative quantity is associated with an increase in symbiotic organ size, we measured the size of the bacteriome in 39 nymphs (n = 6–10/stage). While the bacteriome absolute area doubled at each stage between the first and fourth instars (Figure 2H), its relative area to the thorax remained unchanged between the different developmental stages (Figure 2I).

We first aimed to determine whether the increase in Wolbachia relative quantity during nymphal development persists after the last molt (Figures 3A–D). We thus measured the relative quantity of wCle in fifth instars (5DPF) and newly emerged adults (7DPF) a few days after feeding (data used symbolized by * on Figure 3C (females) and Figure 3D (males); statistics: Figure 3E). We showed a significant increase with development, as log₁₀(RQ) of adult females increased by 0.63 ± 0.42 (equivalent to an increase by 4.57 of the relative quantity) compared to their fifth instar counterparts. No significant effect of sex on bacterial relative quantity was detected in this developmental window.

We then focused on adults and determined the effect of sex, aging, and feeding on the relative quantity of wCle. We chose a weekly feeding protocol to mimic the natural feeding rate of bed bugs (Reinhardt and Siva-Jothy, 2007; Saveer et al., 2021). We thus measured the relative quantity of wCle weekly during 28 days after emergence, in weekly fed or starved adults of each sex (total = 119 adults) and analyzed data using a global statistical model (Figure 3F). The relative quantity of wCle slightly decreased with age in females for both feeding conditions (slope on log₁₀-transformed data = −0.11 equivalent to a decrease of 0.78 each week; Figures 3A,C). In males, while the relative quantity of wCle also decreased with time in starved males (slope on log₁₀-transformed data = 0.16–0.28 = −0.12, equivalent to a decrease of 0.76 each week; Figure 3D), it increased with time in weekly-fed males (slope on log₁₀-transformed data = 0.27–0.11 = 0.16, equivalent to an increase of 1.44 each week; Figure 3B), showing complex interactions between sex, age, and feeding conditions. Overall, males had nevertheless significantly less wCle than females, as log₁₀(RQ) of males was 0.62 ± 0.31 times less the one of females (equivalent to a 4,2-fold difference in relative quantity, Figures 3A–D,F). The log₁₀-transformed absolute quantities of Wolbachia showed similar patterns (Supplementary Figure S6).