Bumblebees (B. terrestris) were collected from the rearing room at the Institute of Apicultural Research, Chinese Academy of Agricultural Science, Beijing, China. Bumblebees were sampled from six independent colonies that were raised in an artificial breeding room (in constant darkness with a temperature of 28 ± 0.5°C and 60 ± 5% relative humidity) and fed fresh frozen pollen and a 50% sugar solution every other day (Gurel and Gosterit 2008). The new emerged queens were aged for 7 days, and half of them were mated with the males, to ensure different stages of virgin and mated queens were collected at the same days old. The queens were kept at room temperature until they were dissected for the collection of the spermathecae. The spermathecae is empty and almost translucent if the bumblebee queen is unfertilized. After mating, the sperm in the spermathecae form an obvious opaque white sphere, which is easy to see (Figure 1) (Duvoisin et al., 1999). The spermathecae of virgin and mating bumblebee queens to be used for morphological observation were dissected in 1 × PBS (10 mM NaH₂PO4/Na₂HPO4, 175 mM NaCl, pH 7.4, Solarbio, Beijing, China), and fixed in 4% paraformaldehyde for 25 min at room temperature, washed with 0.1% Triton X-100 in 1 X PBS (1 X PBT) for 2 × 5 min, and then mounted in mounting medium (90% glycerol). The images were performed with a Leica EZ4W microscope. The spermathecae to be used for RT-qPCR were dissected from the abdomen, and immediately frozen in liquid nitrogen. Tissues from ten bumblebees were pooled as a biological replicate. Three biological replicates were performed for each mating status. All samples were stored at −80°C until they were used for RNA extraction.

Total RNA was extracted from the spermathecae of virgin and the mated queens with TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s instructions. The purity of the RNA was assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) at 260/280 nm, and RNA integrity was screened by 1.5% (w/v) agarose gel electrophoresis. The first strand of cDNA was synthesized according to the instructions of the Reverse Transcription kit (Takara, Dalian, China). The reaction conditions were as follows: 42°C for 30 min, 99°C for 5 min, and 5°C for 5 min, and then the products were stored at −20°C until use.

RNA-sequencing was performed with three biological replicates consisting of three pooled spermatheca samples (n = 10 virgin or mated queens). Following the manufacturer’s recommendations, RNA-sequencing libraries were generated with the NEBNext Ultra™RNA Library Prep Kit from Illumina (NEB, Ipswich, MA, United States), and index codes were added to attribute the sequences to each sample. Using poly-T oligo-attached magnetic beads, mRNA was purified from total RNA. Fragmentation was carried out using divalent cations under high temperature in NEBNext First-Strand Synthesis Reaction Buffer (5 ×). First- and second-strand cDNA was successfully synthesized. The remaining overhangs were converted into blunt ends via exonuclease/polymerase treatment. The short fragments and adapters were linked together. The suitable fragments were chosen as templates, and the subsequent PCR amplification was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and an Index (X) Primer to obtain Index-coded samples. Finally, the PCR products were purified, and the quality of library was assessed on an Agilent Bioanalyzer 2100 system. Index-coded samples were prepared on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v4-cBot-HS (Illumina). The library preparations were sequenced on the Illumina NovaSeq platform.

Adaptors were removed from the raw reads, and low-quality reads were filtered out to obtain the clean reads with NGS QC Toolkit (version 2.3) and Cutadapt (version 1.16). After the quality control, the clean reads were mapped to the reference genome of B. terrestris (version 1.1.1) through Hisat2 software. Each sample was quantified with StringTie. The R package ballgown was used to acquire the gene expression levels. The gene expression level of each transcript was estimated with the fragments per kilobase of transcript per million mapped reads (FPKM) method. FPKM values were directly used to compare the differences of gene expression between various samples. The transcripts with a p-value ≤ 0.05 and the absolute value of the |log2 fold change| ≥ 1 were considered as differentially expressed genes (DEGs).

To annotate the DEGs, Blast2GO software was used to search against the nonredundant (NR) database in NCBI, Clusters of Orthologous Groups of Proteins (COG), Clusters of Orthologous Groups for Eukaryotic Complete Genomes (KOG), and Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups (eggNOG) databases. Furthermore, Gene Ontology (GO terms) and Kyoto Encyclopedia of Genes and Genomes (KEGG pathway) analyses were performed with the default parameters.

Fifteen DEGs were selected as candidates to analyze their expression differences in the spermathecae of queens that response to mating by RT-qPCR using a Stratagene Mx3000 real-time PCR system (Agilent, United States). The primers were designed with Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi), and the primer sequences are shown in Supplementary Table S1. First-strand cDNA samples were diluted (1:10 v/v) with DEPC-treated water. Amplification was carried out in a 20 µL reaction volume containing 10 µL of 10 × TB Green Master Mix (Takara, Dalian, China), 2 µL cDNA, and 0.5 µL of each primer at 10 µm. Quantitative measurements were normalized using β-actin and RP49. The qPCR conditions were as follows: 95°C for 30 s, followed by 40 cycles of 95°C for 5 s, and 63°C for 1 min. RT-qPCR was performed in duplicate on each of three independent biological replicates. All results are presented as the mean ± SEM of the biological replicates. The relative quantities of transcripts were calculated using the comparative Ct method (Livak and Schmittgen, 2001).

Immunostaining was performed as previously described (Li et al., 2013). Briefly, the spermathecae of mating and virgin bumblebee queens were dissected in 1 × PBS and fixed in 4% paraformaldehyde for 25 min at room temperature. Samples were rinsed, washed with 1 × PBT for 2 × 5 min and blocked in 3% BSA in 1 × PBT for 20 min. The primary antibody rabbit anti-PLOD antibody (ProteinTech Group, 1:400) was detected with a fluorescent-conjugated secondary antibody. Incubation with the secondary antibody was performed for 2 h at room temperature, and DAPI (Sigma; 0.1 mg/ml) and phalloidin (Cell Signaling Technology, CST; 8878S) were added during secondary antibody staining. The samples were mounted in mounting medium [70% glycerol containing 2.5% 1, 4-diazabicyclo (2.2.2) octane]. Confocal fluorescence imaging was performed with a Leica SP8 laser-scanning microscope (Leica). For quantification of PLOD1, fluorescence intensity was measured in ImageJ software (Figures 5, 6A).

Spermathecae were homogenized in 1 × PBS. The supernatant was recovered by centrifugation at 3,000 g and 4°C for 10 min. The enzymatic activities of phenoloxidases (POs) in the spermathecae of bumblebee virgins or mated queens were determined with an insect PO ELISA kit (Abmart, AB-3369B, Shanghai, China), and the optical density of each well at 450 nm was measured with an ELISA reader (Molecular Devices SpectraMax i3, United States). The PO activities were calculated according to the calibration curve generated according to the manufacturer’s instructions. Three independent biological replicates were performed for each treatment.

Statistical analysis was performed using the Mann-Whitney U test or one-way ANOVA, and data are presented as the mean ± SEM. The graphs were created with R project software (version 4.0.5). The values and error bars presented in figures represent the means and standard errors of biological replicates.

The raw data quality statistics showed that the Q30 percentages value of the six samples was ranged from 89.07 to 93.65%. The percentages of rRNA were between 0.04 and 0.30% with an average 0.21%. The number of raw reads ranged from 23, 000, 077 to 31, 124, 050 with an average 25, 998, 530 reads. After the mapping of reads to the genome, 11, 069 genes in virgin and 10, 862 genes in mated queens were identified. The number of differentially expressed genes between virgin and mated queen spermathecae was 176. Among the DEGs, 110 (62.5%) were upregulated and 66 (37.5%) were downregulated (Figure 2A and Table 1).

The clustering analysis based on the scale of DEGs FPKM values showed a satisfactory biological replication in the mated and virgin groups (Figure 2B). These genes, such as TIL, LRR, Def, hnRNP, LOC100643810, l-lactate dehydrogenase (L-LDH), LOC100644683, Tret1, nose resistant to fluoxetine protein 6 (NRF), V-type proton ATPase subunit E (NtpE), PO, V-type proton ATPase subunit D (VATPD), serine protease easter-like (SPc), PLOD1, waprin-Phi1-like (WAP), LOC100648738, were dramatically upregulated at 24 h post mating. Two novel genes, LOC100643810 (ncRNA) and LOC100644683 (mRNA), unique to the bumblebee B. terrestris, also showed significantly increased expression at 24 h post mating (Figures 2B, 3).

The GO terms analysis showed that among the 110 upregulated DEGs, those related to nuclear protein-containing complex (12), chromatin (10), and intracellular anatomical structure (10) in cellular component terms were the most abundant, followed by those related to the response to stimulus (12) in biological process term, and the organic cyclic compound binding (10) and catalytic activity (6) in molecular function terms, (Figure 2C). In turn, the 66 downregulated DEGs were associated with the intracellular anatomical structure (20), endomembrane system (18), and supramolecular complex (16) terms (Figure 2C).

The DEGs were also mapped to canonical KEGG pathways to identify possible active biological pathways. Most genes in those pathways were upregulated. Specifically, the gluconeogenesis and neutrophil degranulation categories contained 13 upregulated DEGs, and the thiamine metabolism and purine metabolism were the next most enriched categories. The proton buffering model category contained more than 11 upregulated DEGs, molybdenum cofactor biosynthesis contained 9 upregulated DEGs, and hepatocellular carcinoma and pathways in cancer contained 5 upregulated DEGs, respectively (Table 2).

The transcript levels of 15 DEGs were selected for RT-qPCR analyses to confirm the validity of mating-related genes in the RNA-seq results, which revealed the upregulation of genes in one of the examined tissues (spermathecae of mated and virgin queens). The amplification efficiencies used for correction in all normalized fold-expression analyses ranged from 95.22 to 105.67%, which was within the acceptable range 90–110% (Bustin et al., 2009).

The transcript levels of gene NRF-6, LRRC70, SPc, Def, LOC100644683, hnRNP-1, hnRNP-2, WAP-1, WAP-2, PO2, TIL, and four transcripts of facilitated trehalose transporter genes (Tret1-1, Tret1-2, Tret1-3, and Tret1-4) were significant increased post mating at 24 h (Figure 3). These genes were chosen for RT-qPCR validation. The results of RNA-seq and RT-qPCR showed that 93.3% of these selected genes had consistent expression (Figures 2–4), with the exception of the PLOD1 gene. The correlation analysis showed that RT-qPCR and RNA-seq data was significantly correlated. The RT-qPCR data confirmed some of the differences in mRNA levels first identified in the RNA-sequencing data.

The transcripts of these genes in spermathecae were compared between virgins and mated queens of bumblebee, B. terrestris L., at 1-, 3-, 5-, and 7-days post mating. A total of ten mating-related DEGs (verified in RT-qPCR results at 24 h post mating) were used to verify the expression profiles at different time points within 7 days. The expression levels of the remaining genes gradually decreased. Interestingly, within 7 days, the expression levels of LRRC70, LOC100644683, and PO2, in mating females were gradually decreased post mating, and reached levels similar to those of virgin females, especially at 7 days (Figure 4). However, SPc, Def, Tret1-1, and Tret1-3 were always upregulated in mated females in comparison with virgin females (Figure 4). Although the expression level of TIL, hnRNP-1, and hnRNP-2 tended to fluctuate, they were still always higher than those of virgin females at 1, 3, and 5 days (Figure 4).

A slight decrease in PLOD1 was observed in the mated queens spermathecae, as determined by RT-qPCR (Supplementary Figure S1). To further confirm the expression of PLOD1, we examined the PLOD1 levels in the spermathecae of virgin and mated bumblebee queens by IF, at 24 h post mating. Consistent with the RT-qPCR results, the level of PLOD1 in the mated queens spermathecae (Figure 5F) was dramatically decreased compared to that in the virgin (Figure 5B). The PLOD fluorescence intensity in spermathecae of mated and virgin queen was exhibited in Figure 6A.

The activity of PO in spermathecae of the mated queens significantly increased compared with the virgins (Figure 6B). POs are rate-limiting enzymes, and PO-mediated melanization plays a critical role in the insect immune system (Yassine et al., 2012; Binggeli et al., 2014; Ma et al., 2020). POs were also required for the pea aphid’s defense against bacterial and fungal infections (Xu et al., 2019). Therefore, the function of PO in spermathecae may protect sperm against infection.