Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications

Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications View ORCID Profile Kelley Leung , Frank Becker , Peter Šima , Roland Kreskóci , Elzemiek Geuverink , Thomas V. M. Groot , View ORCID Profile Bart A. Pannebakker , View ORCID Profile Joost van den Heuvel doi: https://doi.org/10.1101/2025.08.12.669960 Kelley Leung 1 Laboratory of Genetics, Wageningen University and Research , PO Box 16, 6700 AA Wageningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Kelley Leung Frank Becker 1 Laboratory of Genetics, Wageningen University and Research , PO Box 16, 6700 AA Wageningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site Peter Šima 2 Koppert s.r.o. , Komárňanská cesta 13, SK-94001 Nové Zámky, Slovakia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Roland Kreskóci 2 Koppert s.r.o. , Komárňanská cesta 13, SK-94001 Nové Zámky, Slovakia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Elzemiek Geuverink 3 Groningen Institute for Evolutionary Life Sciences, University of Groningen , PO Box 11103, 9700 CC, Groningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site Thomas V. M. Groot 4 Koppert Biological Systems , Veilingweg 14, 2651 BE Berkel en Rodenrijs, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site Bart A. Pannebakker 1 Laboratory of Genetics, Wageningen University and Research , PO Box 16, 6700 AA Wageningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Bart A. Pannebakker Joost van den Heuvel 1 Laboratory of Genetics, Wageningen University and Research , PO Box 16, 6700 AA Wageningen, The Netherlands Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Joost van den Heuvel For correspondence: joost.vandenheuvel{at}wur.nl Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract All Hymenoptera (bees, sawflies, ants, wasps) are haplodiploid. Haploid males develop from unfertilized eggs and diploid females from fertilized eggs. Many species also have single locus complementary sex determination (sl-CSD): hemizygosity of a CSD gene results in normal males and heterozygosity in normal females, but homozygosity typically yields sterile diploid males. Populations that lose CSD allelic diversity to critical level, increasingly produce detrimental diploid males. Emblematic of this problem are the bumblebees. Bumblebees are one of the fastest declining animal groups worldwide, with CSD diversity loss playing a likely role, especially in isolated populations. Presently utilized and potential future species in agriculture are also challenged by diploid males in production. Population management of both requires exact knowledge of the CSD locus. Here we report the identity of Bombus terrestris CSD, a bumblebee model that is also a major commercial pollinator. As CSD must be reciprocally homozygous in diploid males and heterozygous in diploid females, we whole-genome-sequenced these for four colonies. From this, we identified B. terrestris CSD (BtCSD) on Chromosome 2: a 26 Kb region with a long non-coding RNA region (lncRNA) and a proximate hypervariable region. Comparative analyses indicated conservation of this region across infraorder Aculeata, including other bumblebees, mason bees, and ants. This supports a deep-rooted ancestral origin for BtCSD despite a lack of homology with con-familial honey bee csd and rapid evolutionary divergence of insect sex determination cascades. Our analyses do not support the lncRNA region as functional CSD, as proposed for other taxa, because of inadequate variation to constitute distinct alleles. Rather, a downstream region is implicated for its massive structural variation and candidate coding genes including Creld and two copies of aCOP . We discuss how this identification of BtCSD has critical implications for hymenopteran evolution and bumblebee population management. Introduction There are more than 150,000 described species in Hymenoptera (sawflies, ants, wasps and bees), with total diversity possibly exceeding a million species ( Forbes et al., 2018 ). They all have haplodiploid sex determination. Haploid males arise from unfertilized eggs, and diploid females from fertilized eggs. Furthermore, numerous taxa including all Aculeata (bees, ants, stinging wasps) have complementary sex determination (CSD) ( Whiting, 1943 ; Asplen et al., 2009 ; Leung & van der Meulen, 2022 ). In CSD, hemizygosity of any allele results in the development of functional males. Heterozygosity results in the development of functional females. Homozygosity however results in the development of aberrant diploid males that are typically sterile ( Heimpel & de Boer, 2008 ). Hence any loss of CSD allelic diversity (e.g. through bottlenecks typical to isolated populations) increases diploid male production. These diploid males typically cannot produce functional female offspring, so a population is predicted to become smaller and more male-biased each generation, even possibly going extinct (diploid male vortex, or DMV) ( Zayed & Packer, 2005 ; Leung & van der Meulen, 2022 ). Measurement and management of CSD parameters to mitigate diploid male production is impossible without exact CSD loci to elucidate taxon-specific effects. Identifying CSD loci however has proven difficult across Hymenoptera, with candidates known in only a few species (ants Vollenhovia emeryi ( Miyakawa & Mikheyev, 2015 ), Linepithema humile ( Pan et al., 2024 ), Ooceraea biroi ( Lacy et al., 2025 ); and orchard bee Osmia bicornis ( Rönneburg et al., 2025 )). A specific csd gene and its mechanism is only identified for A. mellifera (honey bee) ( Beye et al., 2003 ; Hasselmann et al., 2008 ; Otte et al., 2023 ; Suzuki et al., 2023 ). Complicating matters further, species-specific characteristics can obscure under what circumstances the DMV might occur, e.g. mitigative factors like multiple mating and failures of diploid male development and reproduction ( A. mellifera has both; Ratnieks, 1990 ) versus exacerbative factors like monogamous mating and functional diploid males ( Leung & van der Meulen, 2022 ). Therefore, while often cited as a major risk to population survival in insect breeding and conservation, the mechanistic validity of a DMV is unknown for any CSD species. A group uniquely positioned for resolving gaps of CSD knowledge, including roles in population health, are the bumblebees. Like all bees, the ~250 species of widely distributed bumblebees ( Bombus spp.) are prescribed to be sl-CSD as well. ( Duchateau et al., 1994 ; Leung & van der Meulen, 2022 ). Bumblebees are undergoing faster decline than other bees ( Goulson et al., 2008 ; Williams & Osborne, 2009 ; Cameron & Sadd, 2020 ) (e.g. 3 times faster species-richness than other bees in the Netherlands ( Van Dooren, 2019 ). The contribution of taxon-specific aspects of its sl-CSD mechanism is unclear. While there have been suggestions to use diploid males are population health indicator, it has not been incorporated in broader practice for bumblebee management ( Whitehorn et al., 2009 ). For instance, a conservation model for bumblebees incorporates habitat loss, pesticides, and pathogens, but not genetic variation nor sex determination ( Becher et al., 2018 ). However, bumblebees are not known to have compensatory mechanisms for diploid males. Queens are usually monogamous ( Goulson, 2010 ). Unlike some other CSD species, Bombus queens also do not mate-discriminate against more closely related haploid males that would produce sterile diploid males ( Bogo et al., 2018 ), nor against the diploid males themselves ( Lecocq et al., 2017 ). Matings with diploid males very rarely result in any viable female offspring ( Ayabe et al., 2004 ), and workers do not kill diploid male larvae or diploid-male producing queens, as in some bee species ( Woyke, 1963 , 1980 ; Duchateau et al., 1994 ; Alves et al., 2011 ; Vollet-Neto et al., 2017 ). Bumblebees therefore represent fewer complicating variables for direct measurement of CSD diversity impacts on population health, once a specific locus is known. Meanwhile, a few bumblebee species are the most heavily produced group of agricultural pollinators after honey bees ( Osterman et al., 2021 ). European buff-tailed bee Bombus terrestris has been imported worldwide; it and North American counterpart B. impatiens pollinate over 240 crops ( IPBES, 2016 ). Diploid males are a known challenge in commercial breeding. Hives in which they occur are unusable because they displace pollinating female workers. This creates greater cost to breeders to replace these hives ( Duchateau et al., 1994 ; Di Pietro et al., 2022 ). Such applied relevance means that the presence of sl-CSD in B. terrestris ( Duchateau et al., 1994 ) It has also resulted in B. terrestris ’s status as a model system with large, well-annotated genomics datasets, facilitating investigation for CSD locus identification. Here, we report a single genomic region as Bombus terrestris CSD, BtCSD. CSD must be reciprocally homozygous in diploid males and heterozygous in diploid females. A single candidate matching these criteria was detected with de-novo whole genome sequencing and analyses of other available genomics datasets. Phylogenetic investigation indicated elements of remarkable evolutionary conservation across distantly related hymenopterans such as ants and Osmia mason bees, but not more closely related honey bees. In the discussion, we advise on applying knowledge of BtCSD towards future research on re-evaluating patterns of hymenopteran CSD evolution; conservation of wild bumblebee populations; managing present and future commercial populations; and further investigating Bombus CSD molecular mechanism(s). Methods Bumblebee material A commercial producer supplied four, specially selected, early male-producing B. terrestris colonies. Typical outbred colonies do not produce males before the “switch point” ( Duchateau & Velthuis, 1988 ), so all males present in these selected colonies were presumed diploid males ( Di Pietro et al., 2022 ). Female workers, males, and queens were stored at 4°C in absolute ethanol until DNA extraction. Whole genome sequencing Two hindlegs were collected from ten males and the queen of three B. terrestris colonies. Hindlegs of each individual were placed in separate tubes and the ethanol evaporated. They were then frozen in liquid nitrogen and mechanically ground using LabTIE Zirconia beads (Molgen, Veenendaal, The Netherlands). DNA was extracted using a standard cetyltrimethylammoniym bromide (CTAB, Sigma Life Sciences, St. Louis, MO, USA) method ( https://openwetware.org/wiki/Maloof_Lab ). DNA library preparation was done with Hackflex ( Gaio et al., 2022 ). Briefly, DNA samples were quantified using SYBR green (Invitrogen, USA) and diluted to 0.2-0.25 ng/μL. 2 μL of diluted DNA was tagmented and adapter-ligated in a 3 μL reaction using an Illumina Nextera DNA Flex Library Prep Kit (Illumina, San Diego, CA, USA). Libraries were amplified using Illumina TruSeq primers and Robust 2G enzyme (Kapa Biosystems, Cape Town, South Africa). Libraries were size-selected for 300-500bp fragments using solid phase reversible immobilization beads (SPRI, SpeadBead Magnetic Carboxylate Modified Particles, Cytiva, Marlborough, MA, USA). Selected libraries were paired-end sequenced by Novogene Europe (Cambridge, United Kingdom), targeting an average depth of 15X whole genome coverage per diploid male, and 50X coverage per queen. Bioinformatic analyses In addition to our sequenced libraries, we used data from the Short Read Archive to quantify variation in diploid workers (PRJEB52013 ( Barribeau et al., 2022 ), PRJEB49221 ( Kardum Hjort et al., 2022 )) and haploid drones (PRJNA628944 ( Colgan et al., 2022 )) (Supplemental File 1). Raw reads were trimmed using Trimmomatic (v 0.33, adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50, using adapter sequences CTGTCTCTTATACACATCTGACGCTGCCGACGA and GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG ( Bolger et al., 2014 )). Trimmed reads were then mapped to the reference genome of B. terrestris (iyBomTerr1.2, ( Crowley et al., 2023 )) using bwa-mem2 ( Li & Durbin, 2009 ). Both single and paired-end reads were mapped after trimming. Mapped reads filtered for quality and alignments sorted for position resulted in a sorted bam file using samtools (v. 0.1.20, -q 20, ( Li et al., 2009 )). Duplicates were removed with Picard tools (v 2.8.2, (Broad Institute, 2019 )) and Freebayes-parallel was used to call variants (-F 0.25, ( Garrison & Marth, 2012 )). Repeats were removed using bedtools intersect ( Quinlan & Hall, 2010 ). Population genetic parameters We identified heterozygous variants with the worker datasets and false positive sites of heterozygosity with the drone dataset. We used the R library vcfR ( Knaus & Grunwald, 2017 ) in R (4.0.2; ( Team, 2024 ) to calculate minor allele frequencies (MAF), variant quality, and mean depth over all samples. Variants were kept if MAF>0.02, 10 < mean depth 5 ( Kardum Hjort et al., 2022 ) and MAF>0.02, 5 < mean depth 5 ( Barribeau et al., 2022 ). Heterozygosity for all variants was conservatively counted for the drone dataset without MAF of read depth cutoffs. We masked the heterozygosity in both datasets if variants were in 10 bps proximity of drone false heterozygous sites. A variant was called heterozygous if two alleles differed from each other, using the GT field from the function extract.gt(vcf) in the vcfR package. For the diploid males, we calculated pairwise genetic distance as the difference in genotype calls in windows of 100 adjacent variants that were called heterozygous in the queens, using genotype calls from the same GT field (3 < mean depth 1). For the MDS, genetic distance was calculated using alternative reads frequencies. These frequencies were obtained from the read depths, where alternative read depth was taken as the ‘AO’ field from extract.gt(), and reference read depth was taken from the ‘RO’ field from extract.gt(). Once a distance matrix was produced, the MDS plot was produced with cmdscale (R base, 4.0.2;( Team, 2024 ). De novo assembly of drones We used whole genome data from drones ( Colgan et al., 2022 ) to de-novo assemble the CSD region (see Supplemental File 2 for code). First, sorted and deduplicated bam files were used to extract raw reads aligning around the putative CSD locus using “samtools view bams/$pop.sort.bam “NC_063270.1:16252000-16284000” | awk ‘{print $1}’ | sort | uniq”. Seqtk ( https://github.com/lh3/seqtk ) was used to extract a list of reads from raw fastq files. Reads were assembled using spades (v3.13.0, ( Bankevich et al., 2012 ) We used blastn ( Camacho et al., 2009 ) to identify scaffolds as the CSD locus and only kept those that showed a hit with NC_063270.1:16250000-16290000 as a reference. We then produced an assembly using these reads. We used this 40 Kb CSD region to search for reads in the raw data, using a kmer search (bbduk.sh, k=31, hd. https://sourceforge.net/projects/bbmap/ .) Those hits were also used for spades assembly to produce an alternative assembly. These two assemblies were merged using quickmerge ( Chakraborty et al., 2016 ) with the first assembly as ‘reference’ and the bbduk second assembly as query genome (labelled “merged1”). The resulting merged assembly was again filtered using blastn using the 40 Kb CSD region, but the output was filtered for 90% identity and hits of at least 1000 bp length. As many assemblies were still discontinuous, we isolated sequence on the ‘edges’ of the scaffolds by taking the first and last 300 bp of each scaffold (seqkit –r 1:300 and –r −300:−1) and used these sequences to search for matching raw data (bbduk.sh, same parameters above). These reads were assembled for a third assembly using spades. The resulting scaffolds were merged with the previous “merged1” assembly to produce a final merged assembly. From this final assembly, only sequences that contained both the upstream gene LOC125387186 and the downstream gene Creld (LOC100650350) were kept. In total this resulted in 17 assemblies from 51 samples (Supplemental File 4). Kmer synteny blocks Kmer synteny blocks were identified using a custom R script (Supplemental File 3). Briefly, all kmers were identified in the 17 haploid assemblies and in the putative CSD region of the reference genome (Supplemental File 4). Only unique kmers were kept. Common kmers were identified between assemblies and only those which were found more than once (but not in excess in the number of assemblies) were kept, to remove repeats. Kmers were grouped into blocks if the distance between them was less than 20 basepairs. To each common kmer, sample-specific block numbers were added. Common kmers were grouped into numbered blocks, with numbers indicating proximity on the genome. Blocks with more than 400 common kmers were kept as synteny blocks. The relative position of these blocks (start to end) was then plotted. Long distance PCR Using the 18 haploid assemblies, two primer pairs were developed using primer-blast: yF1 – GACCGTCTGCGTTGAAATCG, hr1 – GTGATAGCGCATGGGAGGAA and yF2 – AGCCTTCGGAAATCAGTCCG, hr2 – ATTTCCTCCATGTGCGGCTT. Pair 1 is positioned at 16,275,511-16,283,018 and pair 2 at 16,275,487-16,283,337 in the reference genome. PCR was performed using PCRBIO VeriFi™ Mix. The PCR mix was 12.5 µL Verifi Mix, 1 µL (10µM) forward primer, 1 µL (10µM) reverse primer, 1 µL DNA template and 9,5 µL MQ. The PCR cycle was 1 min of 95°C degrees denaturation and 35 cycles of 95°C denaturation (15s), 60°C degrees annealing (15s), and 72°C degrees extension (5.5 minutes). Ethidium bromide 1% gel electrophoresis was used to visualize PCR results. Following SPRI bead cleanup, PCR products were sent to Plasmidsaurus for long read sequencing using Oxford Nanopore Technology. Amplicon sequencing The structural variation might have obscured or produced false heterozygosity sites in the alignment of the short reads. Therefore, we used amplicon sequencing of four conserved regions that contained variants in queens for Colonies 1 to 3 using 1.8 Kb – 2.5 Kb length products. We amplified these products for each queen, 10 workers and 12 diploid males of each family. We also did this analysis for the fourth family, which was not included in the whole genome sequencing. For amplicon sequencing, primers were developed using the reference genome (regions of low variation in the sequenced populations visualized in IGV ( Robinson et al., 2011 )). Four primer pairs were developed (Ampl1F: Ftail-ATCCTGTGGCTAGACTCGGT, Ampl2F: Ftail-GCCAATTTAGCACGCGAGAG, Ampl3F: Ftail-ATCACCTGTGCACCCATCAC, Ampl4F: Ftail-ATGATTAGGTGGCTTGGGGC, Ampl1R: Rtail-CCTCTTTAGCTCCCTGTGGC, Ampl2R: Rtail-CGAATTTCGCTGCACCTACG, Ampl3R: Rtail – TCACGAAATGCGACACGTCA, Ampl4R: Rtail –TGTTCAGCGTTTCAATCACG). Primers contained Ftail TTTCTGTTGGTGCTGATATTGC and Rtail ACTTGCCTGTCGCTCTATCTTC to enable library prep according to manufacturer protocol with the Ligation sequencing amplicon prep kit (V14, SQK-LSK1114) (Oxford Nanopore Technologies, Oxford, UK). Each of the 4 PCR reactions was 4 µL 5x GoTaq buffer, 0.8 µL DNTPs, 0.2 µL GoTaq polymerase (Promega, Madison, WI, USA), 0.8 µL (10µM) forward primer, 0.8 µL (10µM) reverse primer, 2 µL DNA template, and 11.6 µL MQ. The PCR cycle was 5 min of 95°C denaturation, and 33 cycles of 95°C denaturation (30s), 62°C degrees annealing (30s), and 72°C degrees extension (2.5 minutes). To start library prep with roughly equimolar PCR product, electrophoresis gels were inspected. Because bands were weaker for amplicons 3 and 4, we used 4 µL for amplicons 1 and 2, and 6 µL for amplicons 3 and 4. We then measured and then equalized DNA concentrations of all PCR product pools with PicoGreen (Thermo Fisher Scientific, Waltham, MA, USA). After barcode ligation, DNA concentrations were measured and normalized. The last end prep was sequenced on a Mk1C Nanopore device (Oxford Nanopore Technologies, Oxford, UK), which resulted in 965 Mb of data of fastq passed bases attributed to specific barcodes. The resulting fastq files were aligned to a partial reference genome (NC_063270.1:16,000,000 – 16,500,000) using minimap2 (v2.1, ( Li, 2018 ). Samtools (v. 0.1.20, -q 20, ( Li et al., 2009 )) was used to sort reads and write out in bam format. Picard tools (v 2.8.2 (Broad Institute, 2019 )) was used to add read groups. Variant calls were made with freebayes-parallel ( Garrison & Marth, 2012 ). Variants were filtered using vcfR ( Knaus & Grunwald, 2017 ). Only variants that were heterozygous in >5 individuals (0.4 < allele frequency 8 were kept. Only individuals with an average coverage over 200x for these variants were considered. Homology blast search We used blast to investigate CSD loci region similarity among different hymenopteran taxa. We used sequences NC_063270.1: 16,210,000-16,290,000 ( B. terrestris ), NC_090131.1: 26,510,000-26,590,000 ( Linepithema humile ), NC_060216.1: 14,375,000-14,425,000 ( Osmia bicornis ) and NC_063270.1|:16210000-16290000 ( Ooceraea biroi ). We used blastn ( Camacho et al., 2009 ) to blast all sequences to each other and identify regions larger than 40 basepairs (-word_size 10). A homology plot was made to identify these regions. To detect the CSD loci region across Hymenoptera, we used peptide sequences of the genes flanking B. terrestis CSD (BtCSD) (THUMPD3 XP_003393947.1, CRELD XP_048267559.1 and aCOP XP_048267555.1) in tblastn against available high-quality genome assemblies. The selected databases (Table S1) included all major groups of the infraorder Aculeata as well as known sl-CSD species in the infraorder Proctotrupomorpha. Results Whole genome sequencing identified variation within three colonies of one queen and 10 diploid males. For each queen there was a large degree of heterozygosity (423,527, 392,101 and 421,972 heterozygous variants for Colonies 1, 2 and 3 respectively). Calculated genetic distance ( Supplementary Figure S1 ) indicated high relatedness between diploid males and their queens. The diploid males also had a high average number of heterozygous variants (109,578, 24,100 and 77,511 for each of the three colonies). This confirmed diploid (excluded haploid) status for the male samples. Lower realized sequencing coverage resulted in fewer heterozygous sites in diploid males relative to queens, and in queen 2 relative to other queens. For all heterozygous variants in each queen, the pairwise genetic variation (π) between the 10 diploid males in each family identified BtCSD ( Figure 1 ). That is, in diploid males, the CSD locus should be of low pairwise genetic variation relative to the rest of the genome. One such region was recovered for Colonies 1 and 3 on Chromosome 2, and two such regions were recovered for Colony 2 on Chromosomes 2 and 16. Download figure Open in new tab Figure 1. Genetic variation among 10 DMs within three colonies of B. terrestris . Y axis indicates mean pairwise genetic variation p. Numbers on the X axis indicate chromosome. Mean pairwise genetic variation is quantified for 100 consecutive positions that were heterozygous for the queen, with steps of 25 positions (75% overlap between neighboring windows). Regions of average low heterozygosity in the diploid males were identified in each colony (heterozygosity<0.05, grey area Figure 2 ), and those overlapped between position 16,261,994 and 16,287,800 on Chromsosome 2 (green area Figure 2 ) for all colonies. Average heterozygosity 2.5 Kb windows for two studies that whole genome sequenced B. terrestris workers ( Barribeau et al., 2022 ; Kardum Hjort et al., 2022 ) was highest between 16,278,500 and 16,281,500 on Chromsosome 2, further validating this region as BtCSD. In addition to this region’s high average heterozygosity, the number of variants is also highest there for both worker studies ( Figure 2 ), as expected for a CSD locus. Furthermore, our identified BtCSD locus is homologous to the CSD locus previously reported in the ants L. humile ( Pan et al., 2024 ) and Ooceraea biroi ( Lacy et al., 2025 ) and in another bee, Osmia bicornis ( Rönneburg et al., 2025 ). This corroborates this region as CSD in B. terrestris and indicates that a CSD locus is evolutionary conserved among bees and ants. Download figure Open in new tab Figure 2. Genetic variation around the CSD locus. Panel A-C: heterozygosity in Colonies 1 to 3. Each dot represents the average heterozygosity per variant, the grey area indicates a region with average low heterozygosity among 50 adjacent variants. Panel D-E: heterozygosity in 2.5 Kb windows in workers, from ( Kardum Hjort et al., 2022 ) and ( Barribeau et al., 2022 ) who performed whole genome sequencing. Panel F-H: Number of variants in this study (F) and in ( Kardum Hjort et al., 2022 ) (G) and ( Barribeau et al., 2022 ) (H). The overlapping region of low heterozygosity in Colonies 1-3 overlap with a highly variable region in the two studies of ( Kardum Hjort et al., 2022 ) and ( Barribeau et al., 2022 ), both for heterozygosity per locus and the number of variants present within a 2.5 Kb window. Although the average heterozygosity was low within the 26 Kb candidate region, Family 1 and Family 3 showed some individual loci with high heterozygosity. Visual inspection of these regions for the three colonies ( Supplemental Figure S2 ) showed structural variation within BtCSD. This was indicated by loss of coverage in the alignments; a relatively high abundance of reads that only partially aligned; and reads with abnormal sequence length between left and right read per paired-end. De-novo assembly on shorts reads of whole genome sequenced haploid males ( Colgan et al., 2022 ) better reconstructed this variation. Comparison of synteny blocks in the 17 assembled CSD loci and the reference genome sequence (iyBomTerr1.2) ( Figure 3 ) clearly indicated structural variation. For instance, individual SRR11640395 showed an inverted sequence of synteny blocks ‘rmly’. Some synteny blocks were found in only a few individuals, such as ‘woe’, ‘i’ or ‘g’. In our whole genome sequence alignment, missing coverage was potentially due to these insertions and deletions. Lastly, synteny blocks differed greatly in their relative distances ( Figure 3 ). Download figure Open in new tab Figure 3. Structural variation at the CSD locus. Variation of the CSD locus of de-novo assembled sequences of 17 haploid males is shown (data taken from ( Colgan et al., 2022 )). A synteny block indicates a region that is highly similar between at least two individuals. The upstream and downstream regions are highly conserved. Gene models on top of the figure are indicated by gene or locus names, where a number is listed, it corresponds to an Ensembl gene number. Individuals are ordered along the Y axis according to their genetic distance, based on kmer dissimilarity (Jaccard dissimilarity of presence / absence of kmers), clustered using hierarchical clustering. Samples are named after their NCBI accession number. Heterozygosity and number of variants are indicated (per 1 Kb window) above the gene models. Variation peaks just before the Creld gene. Because short reads were used, we needed to exclude this variation being due to incomplete (mis)assembly. Therefore, to validate the length differences between the synteny blocks, we performed long-distance PCR of synteny block ‘y’ to ‘h’ ( Figure 3 ) on the queen and two diploid males from an additional Family 4. Surprisingly, only one relatively short band was recovered (~7 Kb in the queen, ~10.5 Kb in the two diploid males) ( Supplemental Figure S3 ). Because the two diploid males were sons of the queen, the longer DNA molecule should have also been present in the queen. We hypothesized that the larger band was not visible using gel electrophoresis because the shorter product amplified more efficiently. Long-read sequencing and alignment to identify similar regions between the two alleles found that, as expected, the first 800 and last 2000 base pairs aligned well, while the other sequences could not be aligned (Supplemental File S1). Both PCR products started in synteny block ‘y’ and ended in synteny block ‘h’, although the length of the product differed considerably. Sequencing recovered the long fragment in the queen sample as well, but at lower coverage than the shorter product. The long-distance PCR therefore corroborated large length differences in the BtCSD locus across individuals. We then amplified four additional regions in 86 individuals within BtCSD to further validate differences in heterozygosity between diploid males and workers in four colonies. Similar to the short-read analysis, amplicon sequencing revealed that diploid males had very low heterozygosity levels, while queens and workers showed high levels of heterozygous sites for BtCSD ( Figure 4 ). The positions of those heterozygous sites in diploid males were also similar, compared to the short-read analysis. For instance, in colony 1, heterozygous sites were found upstream (at LOC125387186) while for colony 3 heterozygous variants were identified downstream of BtCSD (at Hcf ). The few heterozygous sites found in diploid males within the CSD region indicates that some variation exists between functional similar alleles. This analysis further corroborates with region as BtCSD, as workers show heterozygosity where diploid males lack it. Download figure Open in new tab Figure 4. Amplicon sequencing results. From left to right 90 loci are indicated spread over 4 PCR products indicated by the grey boxes above. Green boxes indicate genes around the CSD locus, while numbers along that line indicate positions on Chromosome 2 in Mb. For each sample types (W=worker, DM = diploid male) we called genotypes as either homozygous reference (orange), heterozygous (black) or homozygous alternative (light blue). As expected, DMs show much lower degrees of homozygous genotypes, while workers are more likely heterozygous. As the CSD candidates of L. humile, Ooceraea biroi, Osmia bicornis and B. terrestris are in the same genomic region, it was of interest to identify homologous genetic elements among these species. Blasting the four regions against each other notably showed that the larger structural variant in B. terrestris lacked any similarity with these species. The relative distance between the upstream and downstream elements is also much larger for B. terrestris ( Figure 5 ). The genes Creld (LOC100650350) and aCOP (LOC100648009 and LOC125386243, 2 copies) are within the BtCSD region and highest heterozygosity levels were found just downstream of Creld (based on its orientation), similar to the variation found in O. biroi and O. bicornis . Blasting of peptide sequences of the flanking genes of BtCSD ( THUMPD3 (LOC100650110)/ Creld / aCOPI the only coding elements) detected wider conservation of this locus in a range of Hymenoptera genome assemblies. The region is conserved across all Aculeata, except the non-CSD ant Cardiocondyla obscurior , the three tested Vespinae species and one halictid bee (Table S1). This broad conservation is not detected in the more basal infraorder Proctrupomorpha. The synteny of THUMPD3 / Creld / aCOP has been lost and the location of the homologs is at least 0.6Mb apart (Supplemental Table S1). The exception is Bracon brevicornis , where all three genes are located in a region of 106 Kb. Download figure Open in new tab Figure 5. Homology among the CSD loci of B. terrestris, Osmia bicornis, Linepthema humile and Ooceraea biroi . Lines and blocks indicate places in the locus that are shared among all four species. Green box indicates the positions identified by overlapping low heterozygosity levels at BtCSD (same position as green box figure 2 ). Numbers on top of the figure indicate positions on Chromosome 2 in Mb. Discussion The evolutionary implications of an ancestral hymenopteran sex determination mechanism Despite universally controlling maleness versus femaleness, insect sex determination cascades undergo rapid evolution and diversification ( Verhulst et al., 2010 ; Beukeboom & Perrin, 2014 ). Con-generics and even con-specifics do not always have the same mechanism. For example, Musca domestica houseflies have a base XY system, but different populations have different numbers of male sex-determiner genes on different autosomes ( Li et al., 2022 ). Instructor signals (such as any CSD locus) are particularly variable and fast-evolving ( Guerra & Verhulst, 2025 ). Our study detected conservation between BtCSD, annotated CSD regions of ants ( Miyakawa & Mikheyev, 2015 ; Pan et al., 2024 ; Lacy et al., 2025 ) and an orchard bee ( Rönneburg et al., 2025 ), and homologous regions of diverse hymenopteran taxa including other bumblebees (Supplemental Table S1). Furthermore, this Bombus CSD candidate shares no homology with the honey bee ( A. mellifera) csd , a close relative of bumblebees. In A. mellifera, csd is a paralog of the conserved insect transducing sex signal transformer (called feminizer in honey bees), with a hypervariable region corresponding to high allelic diversity ( Beye et al., 2003 ; Beye, 2004 ). While there are related genes in B. terrestris ( transformer and its paralog feminizer1 ), we and others ( Sadd et al., 2015 ) have ruled these out as CSD candidates because they are not universally heterozygous in diploid females and do not have a hypervariable region, as the region we report here is. This collectively suggests that BtCSD has deep evolutionary ancestry and represents a conserved Aculeata CSD mechanism ( Figure 5 ). Apis CSD, and other possible candidates of Bracon brevicornis ( Ferguson et al., 2020 ) and Lysephebius fabrum ( Matthey-Doret et al., 2019 ) parasitoid wasps, would then be independent branches of CSD evolution. This reshapes the evolutionary framework of Hymenopteran CSD, as it means inferences from Apis as the best-studied system are not broadly applicable. The rapid evolution of new csd alleles in the Apis hypervariable region corresponds to hundreds of variants worldwide, which is partly due to variation in the number of asparagine/tyrosine repeats ( Lechner et al., 2013 ; Zareba et al., 2017 ; Seiler & Beye, 2024 ; Ihle et al., 2025 ). While this has some superficial similarity to the structural variation of BtCSD, BtCSD’s hypervariable region needs to be assessed for whether it has comparable, or lower de-novo mutation rate and shared allelic diversity across populations. The latter would greatly reduce the potential of genetic rescue in diploid-male producing bottlenecked populations. However, our study’s finding of a distinct phylogenetic origin of BtCSD from Apis csd indicates that taxonomic proximity should not be considered predictive of CSD conservation, nor similarity of downstream consequences for population evolution and structure. Factoring in Bombus CSD in wild bumblebee conservation There are multiple hypotheses as to why bumblebees are suffering more severe losses than other insect and bee groups. One is that many species are ground nesters, and these sites have been greatly reduced due to industrial agriculture ( Cameron & Sadd, 2020 ). Another considers how some stable species such as B. terrestris and B. lucorum are short-tongued generalist foragers; the loss of longer-tongued species might thus be linked to the specialization to more specific range of flowers, e.g. native species ( Goulson & Darvill, 2004 ). Given its leverage on the detrimental effects of inbreeding in bumblebees, CSD may add to the impacts of these and other factors. Identification of a Bombus CSD here paves the way to assessing its role in conservation status, including remedial possibilities ( Lecocq et al., 2017 ). Since our study suggests likely conservation of Bombus CSD, a foremost objective is measuring the number of CSD alleles in threatened species, as its diversity predicts survival prospects. Documentation of high diploid incidence in populations of once-common North American B. affinis strongly implicates CSD diversity loss in its decline for example ( Mola et al., 2024 ). CSD diversity measurement will identify populations for which conservation efforts should be prioritized due to low CSD diversity, and donor populations needed for novel CSD allele infusion (accounting for possible outbreeding depression effects; ( Gerloff & Schmid-Hempel, 2005 ). Unfortunately, it may also identify populations that cannot be saved because the minimum level of CSD diversity to sustain a population no longer exists. These designations need to be made carefully, however, as historical context is needed to establish what sustainable levels of CSD diversity are. For example, general inbreeding in B. veteranus in the Netherlands, where it is now critically endangered ( Slikboer et al., 2023 ), was indeed confirmed using museum specimens ( Maebe et al., 2013 ). Targeted sequencing, e.g. some marker-level assay for Bombus CSD, would give a higher resolution answer on how CSD diversity impacts population health over time. An important observation of our study however is that the colonies inbred for BtCSD were outbred (had high heterozygosity) throughout the rest of genome ( Figures 1 and 2 ). This underscores the importance of knowing CSD loci and measuring it directly, as random genome-wide markers are a poor proxy for inferring diversity for this single gene with outsized impact on population health. Factoring Bombus CSD in producing bumblebees Currently, commercial bumblebee production of commonly used B. terrestris and B. impatiens seems to be stable ( Osterman et al., 2021 ). The exact costs of unusable diploid male-producing colonies are not disclosed by commercial breeders. However, by destroying these colonies, they are performing a human-mediated form of CSD diploid male vortex mitigation. In future, commercial breeders may electively evaluate the CSD-related aspects of population health in their stocks by measuring and monitoring the amount of CSD allelic diversity. Even with diploid-male colony destruction, if the initial pool of CSD alleles is insufficient and crosses not sufficiently variable, an isolated captive breeding stock may have a limited “lifespan”. Critically, knowledge of this Bombus CSD candidate will be essential to developing new local bumblebee strains or species for large-scale agricultural pollination. The rapidly growing interest in this industry is motivated by a desire to reduce reliance on B. terrestris imports, which at times has led to resource competition with native bees including other bumblebees ( Ishii et al., 2008 ; Williams & Osborne, 2009 ; Goka, 2010 ); disease spillover of pathogens and parasites ( Smith-Ramírez et al., 2023 ); and introgression into native populations ( Kardum Hjort et al., 2022 ; Franchini et al., 2023 ). The last endangers native genetic diversity and its linkage to local adaptation; it can also produce sterile hybrids ( Kondo et al., 2009 ) However, poor management of CSD diversity may prevent successful development of new commercial populations. This candidate may be used to design new breeding protocols that ensure sufficient CSD alleles in a starting population; crosses that maintains this diversity and prevents any costly diploid male production; and identifies sources of distinct alleles for genetic rescue if needed (e.g. different wild populations). Functional investigation of Bt CSD A precursor to any of the above fundamental and applied avenues for this Bombus CSD candidate is more mechanistic knowledge. Importantly, BtCSD still requires functional validation. As CSD is a feminization signal, preventing its expression in heterozygous embryos through RNA interference would divert diploid female development to male development (e.g. as with other hymenopteran sex determiners: ( Zou et al., 2020 ; Pan et al., 2024 ). The exact target(s) in the large 2Kb region are still to be determined. While a Bombus CSD candidate was previously suggested using a B. terrestris linkage map, it was derived from randomly distributed RAPD markers ( Wilfert et al., 2006 ; Stolle et al., 2011 ) and could not be mapped back here to narrow a region of interest. However, a hint may come from the ant L. humile , as targeting the non-coding region in this ant resulted in only a weak knockdown effect (expression of male splice transformer in ~10% of treated heterozygous eggs) ( Pan et al., 2024 ). This long non-coding RNA does not show any homology to the BtCSD and CSD locus of Osmia bicornis . This suggests that lncRNA candidates are linked to CSD loci rather than being CSD outright. As such, the most promising candidate targets could be the downstream coding genes Creld and the two copies of a COP in the hypervariable region. Creld is positioned with its C-terminal region extending into the structurally variable region ( Figure 3 ). A proper functional testthat elucidates the genetic element for functional sex determination in BtCSD (i.e. a knockdown-or-knockout of a candidate that results in heterozygous male development) will allow for a better comparison between the bumblebee and the ant CSD locus. Once the functional gene(s) have been identified, a major aspect of the functionality of CSD is the amount of sequence difference needed to qualify different alleles. For example, the amount of difference needed between alleles for heterozygosity (heterotrimer formation) versus homozygosity (homotrimer formation) is still debated even within the well-studied and high-CSD diversity Apis system ( Otte et al., 2023 ). There is a complex model requiring both differences of length and a threshold number of amino acid mismatches in the hypervariable region ( Lechner et al., 2013 ); and simpler models arguing that a minimum of 3-5 amino acid differences ( Beye et al., 2013 ), or even one ( Mroczek et al., 2022 ), is sufficient. Testing this for bumblebees would require controlled crosses with alleles falling across a gradient of relatedness. This in turn needs sequencing of CSD candidate to cross individuals with the right allelic combinations. As a synthesis, identification of Bombus CSD is the launching point for many branches of bumblebee research. This ranges from to better definition of the extinction risk of CSD in Bombus population dynamics and remedial possibilities; to improved breeding of existing and to-be-founded commercial populations; to further delineation of the mechanism itself. Author Contributions Joost van den Heuvel, Kelley Leung and Bart A. Pannebakker conceived this study. Elzemiek Geuverink, Peter Šima, Thomas V. M. Groot and Roland Kreskóci advised on its design. Peter Šima, and Roland Kreskóci provided specimens. Frank Becker and Joost van den Heuvel collected data. Joost van den Heuvel and Elzemiek Geuverink analyzed data. Kelley Leung secured funding. Kelley Leung and Joost van den Heuvel wrote the initial draft of the manuscript. All authors discussed results and contributed to the writing of other versions. Data availability statement All raw data have been uploaded to the NCBI Short Read Archive in Bioproject PRJNA1296802. Supplemental figures Download figure Open in new tab Figure S1. Genetic distance between sampled B. terrestris colonies. Multidimensional scaling of genetic distance among the 33 samples, including 1 queen and 10 diploid males per family. Dimension 1 explains 32.1% of the total variation, while 25.6% is explained by dimension 2. Different colonies are indicated by varying symbols and color. The larger symbols indicate queens. This plot is indicative of all family members being genetically relatively close to each other. Download figure Open in new tab Figure S2. Integrative genome viewer output for the CSD locus on Chromosome 2. From top to bottom alignments of the queens of Family 1, 2 and 3 are shown. Alignment is lacking around 16,260 Kb, but also partly around 16,265 for Queen 1 in the upstream area. In the downstream area coverage is lacking between 16,275 and 16,280 Kb. The colors in each upper panel (coverage) panels indicate variation compared to the reference genome. Very clearly, a lot of genomic variation is inferred around the regions with missing coverage. Download figure Open in new tab Figure S3. Gel electrophoresis result from long distance PCR on the queen and two diploid males from Family 4 produced one long fragment in two diploid males (dm4-1 and dm4-2) and one shorter fragment in the queen (Q4). The lengths of two separate primer pair amplified seemingly the same region. Table S1. Detection of THUMPD3/Creld/acOPI genomic region across Hymenoptera. Species were selected based on availability of high-quality genome, (if known) presence/absence of CSD mechanism and phylogenetic position. (Not included here because of size; see attached excel file) Acknowledgements This work was funded by the Dutch Research Council (NWO) with VI.Veni.232.066 awarded to Kelley Leung. We thank our commercial partner for providing materials for sequencing. The authors declare no conflict of interest. Funder Information Declared Dutch Research Council, https://ror.org/04jsz6e67 , VI.Veni.232.066 Footnotes Author emails: Kelley Leung ( kelley.leung{at}wur.nl ) Frank Becker ( frank.becker{at}wur.nl ) Peter Šima ( psima{at}koppert.sk ) Roland Kreskóci ( RKreskoci{at}koppert.sk ) Elzemiek Geuverink ( e.geuverink{at}rug.nl ) Thomas V. M. Groot ( tgroot{at}kopper.nl ) Bart A. Pannebakker ( bart.pannebakker{at}wur.nl ) References ↵ Alves , D. A. , Imperatriz-Fonseca , V. L. , Francoy , T. M. , Santos-Filho , P. S. , Billen , J. , & Wenseleers , T. ( 2011 ). Successful maintenance of a stingless bee population despite a severe genetic bottleneck . Conservation Genetics , 12 ( 3 ), 647 – 658 . doi: 10.1007/s10592-010-0171-z OpenUrl CrossRef ↵ Asplen , M. K. , Whitfield , J. B. , De Boer , J. G. , & Heimpel , G. E. ( 2009 ). Ancestral state reconstruction analysis of hymenopteran sex determination mechanisms . Journal of Evolutionary Biology , 22 ( 8 ), 1762 – 1769 . doi: 10.1111/j.1420-9101.2009.01774.x OpenUrl CrossRef PubMed ↵ Ayabe , T. , Hoshiba , H. , & Ono , M. ( 2004 ). Cytological evidence for triploid males and females in the bumblebee, Bombus terrestris . Chromosome Research , 12 ( 3 ), 215 – 223 . doi: 10.1023/B:CHRO.0000021880.83639.4b OpenUrl CrossRef PubMed ↵ Bankevich , A. , Nurk , S. , Antipov , D. , Gurevich , A. A. , Dvorkin , M. , Kulikov , A. S. , Pevzner , P. A. ( 2012 ). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing . J Comput Biol , 19 ( 5 ), 455 – 477 . doi: 10.1089/cmb.2012.0021 OpenUrl CrossRef PubMed ↵ Barribeau , S. M. , Schmid-Hempel , P. , Walser , J. C. , Zoller , S. , Berchtold , M. , Schmid-Hempel , R. , & Zemp , N. ( 2022 ). Genetic variation and microbiota in bumble bees cross-infected by different strains of C. bombi . PLOS ONE , 17 ( 11 ), e0277041 . doi: 10.1371/journal.pone.0277041 OpenUrl CrossRef ↵ Becher , M. A. , Twiston-Davies , G. , Penny , T. D. , Goulson , D. , Rotheray , E. L. , & Osborne , J. L. ( 2018 ). Bumble-BEEHAVE: A systems model for exploring multifactorial causes of bumblebee decline at individual, colony, population and community level . Journal of Applied Ecology , 55 ( 6 ), 2790 – 2801 . doi: 10.1111/1365-2664.13165 OpenUrl CrossRef PubMed ↵ Beukeboom , L. W. , & Perrin , N. ( 2014 ). The Evolution of Sex Determination : Oxford University Press . ↵ Beye , M. ( 2004 ). The dice of fate: the csd gene and how its allelic composition regulates sexual development in the honey bee, Apis mellifera . Bioessays , 26 ( 10 ), 1131 – 1139 . doi: 10.1002/bies.20098 OpenUrl CrossRef PubMed Web of Science ↵ Beye , M. , Hasselmann , M. , Fondrk , M. K. , Page , R. E. , & Omholt , S. W. ( 2003 ). The gene csd is the primary signal for sexual development in the honeybee and encodes an SR-type protein . Cell , 114 ( 4 ), 419 – 429 . doi: 10.1016/s0092-8674(03)00606-8 OpenUrl CrossRef PubMed Web of Science ↵ Beye , M. , Seelmann , C. , Gempe , T. , Hasselmann , M. , Vekemans , X. , Fondrk , M. K. , & Page, Robert E. ( 2013 ). Gradual molecular evolution of a sex determination switch through incomplete penetrance of femaleness . Current Biology , 23 ( 24 ), 2559 – 2564 . doi: 10.1016/j.cub.2013.10.070 OpenUrl CrossRef PubMed ↵ Bogo , G. , de Manincor , N. , Fisogni , A. , Galloni , M. , Zavatta , L. , & Bortolotti , L. ( 2018 ). No evidence for an inbreeding avoidance system in the bumble bee Bombus terrestris . Apidologie , 49 ( 4 ), 473 – 483 . doi: 10.1007/s13592-018-0575-1 OpenUrl CrossRef ↵ Bolger , A. M. , Lohse , M. , & Usadel , B. ( 2014 ). Trimmomatic: a flexible trimmer for Illumina sequence data . Bioinformatics , 30 ( 15 ), 2114 – 2120 . doi: 10.1093/bioinformatics/btu170 OpenUrl CrossRef PubMed Web of Science ↵ Broad Institute . ( 2019 ). Picard toolkit . In Broad Institute, GitHub repository . ↵ Camacho , C. , Coulouris , G. , Avagyan , V. , Ma , N. , Papadopoulos , J. , Bealer , K. , & Madden , T. L. ( 2009 ). BLAST+: architecture and applications . BMC Bioinformatics , 10 ( 1 ), 421 . doi: 10.1186/1471-2105-10-421 OpenUrl CrossRef PubMed ↵ Cameron , S. A. , & Sadd , B. M. ( 2020 ). Global trends in bumble bee health . Annu Rev Entomol , 65 , 209 – 232 . doi: 10.1146/annurev-ento-011118-111847 OpenUrl CrossRef PubMed ↵ Chakraborty , M. , Baldwin-Brown , J. G. , Long , A. D. , & Emerson , J. J. ( 2016 ). Contiguous and accurate de novo assembly of metazoan genomes with modest long read coverage . Nucleic Acids Research , 44 ( 19 ), e147 – e147 . doi: 10.1093/nar/gkw654 OpenUrl CrossRef PubMed ↵ Colgan , T. J. , Arce , A. N. , Gill , R. J. , Ramos Rodrigues , A. , Kanteh , A. , Duncan , E. J. , Wurm , Y. ( 2022 ). Genomic signatures of recent adaptation in a wild bumblebee . Mol Biol Evol , 39 ( 2 ). doi: 10.1093/molbev/msab366 OpenUrl CrossRef PubMed ↵ Crowley , L. , Sivell , O. , Sivell , D. , null , n. , null , n. , null , n. , null , n. ( 2023 ). The genome sequence of the Buff-tailed Bumblebee, Bombus terrestris (Linnaeus, 1758) [version 1; peer review: 2 approved, 1 approved with reservations] . Wellcome Open Research , 8 ( 161 ). doi: 10.12688/wellcomeopenres.19248.1 OpenUrl CrossRef PubMed ↵ Di Pietro , V. , Ferreira , H. M. , Van Oystaeyen , A. , Wäckers , F. , Wenseleers , T. , & Oliveira , R. C. ( 2022 ). Distinct colony types caused by diploid Male production in the buff-tailed bumblebee Bombus terrestris . Frontiers in Ecology and Evolution , Volume 10 2022 . doi: 10.3389/fevo.2022.844251 OpenUrl CrossRef ↵ Duchateau , M. J. , Hoshiba , H. , & Velthuis , H. H. W. ( 1994 ). Diploid males in the bumble bee Bombus terrestris . Entomologia Experimentalis et Applicata , 71 ( 3 ), 263 – 269 . doi: 10.1111/j.1570-7458.1994.tb01793.x OpenUrl CrossRef ↵ Duchateau , M. J. , & Velthuis , H. H. W. ( 1988 ). Development and Reproductive Strategies in Bombus terrestris Colonies . Behaviour , 107 ( 3-4 ), 186 – 207 . doi: 10.1163/156853988×00340 OpenUrl CrossRef ↵ Ferguson , K. B. , Pannebakker , B. A. , Centurion , A. , van den Heuvel , J. , Nieuwenhuis , R. , Becker , F. F. M. , Verhulst , E. C. ( 2020 ). Bracon brevicornis genome showcases the potential of linked-read sequencing in identifying a putative complementary sex determiner Gene . Genes (Basel) , 11 ( 12 ). doi: 10.3390/genes11121390 OpenUrl CrossRef ↵ Forbes , A. A. , Bagley , R. K. , Beer , M. A. , Hippee , A. C. , & Widmayer , H. A. ( 2018 ). Quantifying the unquantifiable: why Hymenoptera, not Coleoptera, is the most speciose animal order . BMC Ecology , 18 ( 1 ), 21 . doi: 10.1186/s12898-018-0176-x OpenUrl CrossRef PubMed ↵ Franchini , P. , Fruciano , C. , Wood , T. J. , Shastry , V. , Goulson , D. , Hughes , W. O. H. , & Jones , J. C. ( 2023 ). Limited introgression from non-native commercial strains and signatures of adaptation in the key pollinator Bombus terrestris . Molecular Ecology , 32 ( 21 ), 5709 – 5723 . doi: 10.1111/mec.17151 OpenUrl CrossRef ↵ Gaio , D. , Anantanawat , K. , To , J. , Liu , M. , Monahan , L. , & Darling , A. E. ( 2022 ). Hackflex: low-cost, high-throughput, Illumina Nextera Flex library construction . Microb Genom , 8 ( 1 ). doi: 10.1099/mgen.0.000744 OpenUrl CrossRef ↵ Garrison , E. , & Marth , G. ( 2012 ). Haplotype-based variant detection from short-read sequencing . arXiv preprint arxiv: 1207.3907 . ↵ Gerloff , C. U. , & Schmid-Hempel , P. ( 2005 ). Inbreeding depression and family variation in a social insect, Bombus terrestris (Hymenoptera: Apidae) . OIKOS , 111 ( 1 ), 67 – 80 . doi: 10.1111/j.0030-1299.2005.13980.x OpenUrl CrossRef ↵ Goka , K. ( 2010 ). Introduction to the Special feature for ecological risk assessment of introduced bumblebees: Status of the European bumblebee, Bombus terrestris, in Japan as a beneficial pollinator and an invasive alien species . Applied Entomology and Zoology , 45 ( 1 ), 1 – 6 . doi: 10.1303/aez.2010.1 OpenUrl CrossRef ↵ Goulson , D. ( 2010 ). Bumblebees : behaviour, ecology, and conservation (2nd ed .). Oxford [England] : Oxford University Press . ↵ Goulson , D. , & Darvill , B. ( 2004 ). Niche overlap and diet breadth in bumblebees; are rare species more specialized in their choice of flowers? Apidologie , 35 ( 1 ), 55 – 63 . Retrieved from doi: 10.1051/apido:2003062 OpenUrl CrossRef ↵ Goulson , D. , Lye , G. C. , & Darvill , B. ( 2008 ). Decline and conservation of bumble bees . Annu Rev Entomol , 53 , 191 – 208 . doi: 10.1146/annurev.ento.53.103106.093454 OpenUrl CrossRef PubMed Web of Science ↵ Guerra , F. , & Verhulst , E. C. ( 2025 ). To each their own: the diversity of primary sex determination signals in insects . Current Opinion in Insect Science , 101414 . doi: 10.1016/j.cois.2025.101414 OpenUrl CrossRef ↵ Hasselmann , M. , Vekemans , X. , Pflugfelder , J. , Koeniger , N. , Koeniger , G. , Tingek , S. , & Beye , M. ( 2008 ). Evidence for convergent nucleotide evolution and high allelic turnover rates at the complementary sex determiner gene of Western and Asian honeybees . Mol Biol Evol , 25 ( 4 ), 696 – 708 . doi: 10.1093/molbev/msn011 OpenUrl CrossRef PubMed Web of Science ↵ Heimpel , G. E. , & de Boer , J. G. ( 2008 ). Sex determination in the Hymenoptera . Annual Review of Entomology . doi: 10.1146/annurev.ento.53.103106.093441 OpenUrl CrossRef PubMed Web of Science ↵ Ihle , K. E. , Payne , A. N. , Elsik , C. G. , & Bilodeau , L. ( 2025 ). Scientific note updating allelic nomenclature standards of the highly diverse complementary sex-determiner locus in honey bees . Apidologie , 56 . doi: 10.1007/s13592-025-01161-4 OpenUrl CrossRef ↵ IPBES . ( 2016 ). The assessment report of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services on pollinators, pollination and food production .. Retrieved from Bonn, Germany : ↵ Ishii , H. S. , Kadoya , T. , Kikuchi , R. , Suda , S.-I. , & Washitani , I. ( 2008 ). Habitat and flower resource partitioning by an exotic and three native bumble bees in central Hokkaido, Japan . Biological Conservation , 141 ( 10 ), 2597 – 2607 . doi: 10.1016/j.biocon.2008.07.029 OpenUrl CrossRef ↵ Kardum Hjort , C. , Paris , J. R. , Olsson , P. , Herbertsson , L. , de Miranda , J. R. , Dudaniec , R. Y. , & Smith , H. G. ( 2022 ). Genomic divergence and a lack of recent introgression between commercial and wild bumblebees (Bombus terrestris) . Evolutionary Applications , 15 ( 3 ), 365 – 382 . OpenUrl PubMed ↵ Knaus , B. J. , & Grunwald , N. J. ( 2017 ). vcfr: a package to manipulate and visualize variant call format data in R . Mol Ecol Resour , 17 ( 1 ), 44 – 53 . doi: 10.1111/1755-0998.12549 OpenUrl CrossRef PubMed ↵ Kondo , N. I. , Yamanaka , D. , Kanbe , Y. , Kunitake , Y. K. , Yoneda , M. , Tsuchida , K. , & Goka , K. ( 2009 ). Reproductive disturbance of Japanese bumblebees by the introduced European bumblebee Bombus terrestris . Naturwissenschaften , 96 ( 4 ), 467 – 475 . doi: 10.1007/s00114-008-0495-4 OpenUrl CrossRef PubMed ↵ Lacy , K. D. , Lee , J. , Rozen-Gagnon , K. , Wang , W. , Carroll , T. S. , & Kronauer , D. J. C. ( 2025 ). Heterozygosity at a conserved candidate sex determination locus is associated with female development in the clonal raider ant (Ooceraea biroi) . BioRxiv . doi: 10.1101/2025.01.24.634795 OpenUrl Abstract / FREE Full Text ↵ Lechner , S. , Ferretti , L. , Schöning , C. , Kinuthia , W. , Willemsen , D. , & Hasselmann , M. ( 2013 ). Nucleotide variability at its limit? Insights into the number and evolutionary dynamics of the sex-determining specificities of the honey bee Apis mellifera . Molecular Biology and Evolution , 31 ( 2 ), 272 – 287 . doi: 10.1093/molbev/mst207 OpenUrl CrossRef PubMed ↵ Lecocq , T. , Gerard , M. , Maebe , K. , Brasero , N. , Dehon , L. , Smagghe , G. , Michez , D. ( 2017 ). Chemical reproductive traits of diploid Bombus terrestris males: Consequences on bumblebee conservation . Insect Sci , 24 ( 4 ), 623 – 630 . doi: 10.1111/1744-7917.12332 OpenUrl CrossRef PubMed ↵ Leung , K. , & van der Meulen , H. ( 2022 ). Revisiting the hymenopteran diploid male vortex: a review of avoidance mechanisms and incidence . Entomologia Experimentalis et Applicata , 170 ( 12 ), 1010 – 1031 . doi: 10.1111/eea.13234 OpenUrl CrossRef ↵ Li , H. ( 2018 ). Minimap2: pairwise alignment for nucleotide sequences . Bioinformatics , 34 ( 18 ), 3094 – 3100 . doi: 10.1093/bioinformatics/bty191 OpenUrl CrossRef PubMed ↵ Li , H. , & Durbin , R. ( 2009 ). Fast and accurate short read alignment with Burrows-Wheeler transform . Bioinformatics , 25 ( 14 ), 1754 – 1760 . doi: 10.1093/bioinformatics/btp324 OpenUrl CrossRef PubMed Web of Science ↵ Li , H. , Handsaker , B. , Wysoker , A. , Fennell , T. , Ruan , J. , Homer , N. , Durbin , R. ( 2009 ). The Sequence Alignment/Map format and SAMtools . Bioinformatics , 25 ( 16 ), 2078 – 2079 . doi: 10.1093/bioinformatics/btp352 OpenUrl CrossRef PubMed Web of Science ↵ Li , X. , Lin , F. , van de Zande , L. , & Beukeboom , L. W. ( 2022 ). Strong variation in frequencies of male and female determiners between neighboring housefly populations . Insect Science , 29 ( 5 ), 1470 – 1482 . doi: 10.1111/1744-7917.13017 OpenUrl CrossRef PubMed ↵ Maebe , K. , Meeus , I. , Maharramov , J. , Grootaert , P. , Michez , D. , Rasmont , P. , & Smagghe , G. ( 2013 ). Microsatellite analysis in museum samples reveals inbreeding before the regression of Bombus veteranus . Apidologie , 44 , 188 – 197 . OpenUrl ↵ Matthey-Doret , C. , van der Kooi , C. J. , Jeffries , D. L. , Bast , J. , Dennis , A. B. , Vorburger , C. , & Schwander , T. ( 2019 ). Mapping of multiple complementary sex determination loci in a parasitoid wasp . Genome Biol Evol , 11 ( 10 ), 2954 – 2962 . doi: 10.1093/gbe/evz219 OpenUrl CrossRef PubMed ↵ Miyakawa , M. O. , & Mikheyev , A. S. ( 2015 ). QTL mapping of sex determination loci supports an ancient pathway in ants and honey bees . PLoS Genet , 11 ( 11 ), e1005656 . doi: 10.1371/journal.pgen.1005656 OpenUrl CrossRef ↵ Mola , J. M. , Pearse , I. S. , Boone , M. L. , Evans , E. , Hepner , M. J. , Jean , R. P. , Koch , J. B. U. ( 2024 ). Range-wide genetic analysis of an endangered bumble bee (Bombus affinis, Hymenoptera: Apidae) reveals population structure, isolation by distance, and low colony abundance . J Insect Sci , 24 ( 19 ), ieae041 . doi: 10.1093/jisesa/ieae041 OpenUrl CrossRef ↵ Mroczek , R. , Laszkiewicz , A. , Blazej , P. , Adamczyk-Weglarzy , K. , Niedbalska-Tarnowska , J. , & Cebrat , M. ( 2022 ). New insights into the criteria of functional heterozygosity of the Apis mellifera complementary sex determining gene-Discovery of a functional allele pair differing by a single amino acid . PLOS ONE , 17 ( 8 ), e0271922 . doi: 10.1371/journal.pone.0271922 OpenUrl CrossRef PubMed ↵ Osterman , J. , Aizen , M. A. , Biesmeijer , J. C. , Bosch , J. , Howlett , B. G. , Inouye , D. W. , Paxton , R. J. ( 2021 ). Global trends in the number and diversity of managed pollinator species . Agriculture, Ecosystems & Environment , 322 , 107653 . doi: 10.1016/j.agee.2021.107653 OpenUrl CrossRef ↵ Otte , M. , Netschitailo , O. , Weidtkamp-Peters , S. , Seidel , C. A. M. , & Beye , M. ( 2023 ). Recognition of polymorphic Csd proteins determines sex in the honeybee . Science Advances , 9 ( 40 ), eadg4239. doi : doi: 10.1126/sciadv.adg4239 OpenUrl CrossRef PubMed ↵ Pan , Q. , Darras , H. , & Keller , L. ( 2024 ). LncRNA gene ANTSR coordinates complementary sex determination in the Argentine ant . Science Advances , 10 ( 22 ), eadp1532. doi : doi: 10.1126/sciadv.adp1532 OpenUrl CrossRef PubMed ↵ Quinlan , A. R. , & Hall , I. M. ( 2010 ). BEDTools: a flexible suite of utilities for comparing genomic features . Bioinformatics , 26 ( 6 ), 841 – 842 . doi: 10.1093/bioinformatics/btq033 OpenUrl CrossRef PubMed Web of Science ↵ Ratnieks , F. L. W. ( 1990 ). The evolution of polyandry by queens in social Hymenoptera: the significance of the timing of removal of diploid males . Behavioral Ecology and Sociobiology , 26 ( 5 ), 343 – 348 . doi: 10.1007/BF00171100 OpenUrl CrossRef Web of Science ↵ Robinson , J. T. , Thorvaldsdottir , H. , Winckler , W. , Guttman , M. , Lander , E. S. , Getz , G. , & Mesirov , J. P. ( 2011 ). Integrative Genomics Viewer . Nat Biotechnol , 29 ( 1 ), 24 – 26 . doi: 10.1038/nbt.1754 OpenUrl CrossRef PubMed Web of Science ↵ Rönneburg , T. , Taliadoros , D. , Olsson , A. , Magnusson , S. , Huser , L. , Fuad , M. N. I. B. N. , Webster , M. T. ( 2025 ). Genetic mapping in the red mason bee Osmia bicornis implicates ANTSR as an ancient sex-determining locus in bees and ants . BioRxiv , 2025.2007.2003.662907. doi: 10.1101/2025.07.03.662907 OpenUrl Abstract / FREE Full Text ↵ Sadd , B. M. , Barribeau , S. M. , Bloch , G. , de Graaf , D. C. , Dearden , P. , Elsik , C. G. , Worley , K. C. ( 2015 ). The genomes of two key bumblebee species with primitive eusocial organization . Genome Biol , 16 ( 1 ), 76 . doi: 10.1186/s13059-015-0623-3 OpenUrl CrossRef PubMed ↵ Seiler , J. , & Beye , M. ( 2024 ). Honeybees’ novel complementary sex-determining system: function and origin . Trends Genet , 40 ( 11 ), 969 – 981 . doi: 10.1016/j.tig.2024.08.001 OpenUrl CrossRef PubMed ↵ Slikboer , L. , Smit , J. T. , van der Jagt , L. , & Huijzers , G. W. ( 2023 ). Panda van de Nederlandse Delta’: de zandhommel Bombus veteranus in Nederland (Hymenoptera: Apidae) . entomologische berichten , 83 ( 6 ), 162 – 170 . OpenUrl ↵ Smith-Ramírez , C. , Rendón-Funes , A. , Leiva , M. , Arbetman , M. , Aizen , M. , & Agüero , L. ( 2023 ). Non-compliance with the World Trade Organization agreements by exporters of the European bumblebee, Bombus terrestris . Sustainability: Science, Practice and Policy , 19 ( 1 ), 2256173 . doi: 10.1080/15487733.2023.2256173 OpenUrl CrossRef ↵ Stolle , E. , Wilfert , L. , Schmid-Hempel , R. , Schmid-Hempel , P. , Kube , M. , Reinhardt , R. , & Moritz , R. F. ( 2011 ). A second generation genetic map of the bumblebee Bombus terrestris (Linnaeus, 1758) reveals slow genome and chromosome evolution in the Apidae . BMC Genomics , 12 , 48 . doi: 10.1186/1471-2164-12-48 OpenUrl CrossRef PubMed ↵ Suzuki , Y. , Yamada , T. , & Suzuki , M. G. ( 2023 ). In vitro comparison of sex-specific splicing efficiencies of fem pre-mRNA under monoallelic and heteroallelic conditions of csd, a master sex-determining gene in the honeybee . J Dev Biol , 11 ( 1 ). doi: 10.3390/jdb11010010 OpenUrl CrossRef ↵ Team , R. C. ( 2024 ). R: A language and enviroment for statistical computing . In. Vienna: R Foundation for Statistical Computing . ↵ Van Dooren , T. J. M. ( 2019 ). Assessing species richness trends: Declines of bees and bumblebees in the Netherlands since 1945 . Ecol Evol , 9 ( 23 ), 13056 – 13068 . doi: 10.1002/ece3.5717 OpenUrl CrossRef ↵ Verhulst , E. C. , van de Zande , L. , & Beukeboom , L. W. ( 2010 ). Insect sex determination: it all evolves around transformer . Curr Opin Genet Dev , 20 ( 4 ), 376 – 383 . doi: 10.1016/j.gde.2010.05.001 OpenUrl CrossRef PubMed Web of Science ↵ Vollet-Neto , A. , Oliveira , R. C. , Schillewaert , S. , Alves , D. A. , Wenseleers , T. , Nascimento , F. S. , Ratnieks , F. L. W. ( 2017 ). Diploid male production results in queen death in the stingless bee Scaptotrigona depilis . J Chem Ecol , 43 ( 4 ), 403 – 410 . doi: 10.1007/s10886-017-0839-7 OpenUrl CrossRef PubMed ↵ Whitehorn , P. R. , Tinsley , M. C. , Brown , M. J. F. , Darvill , B. , & Goulson , D. ( 2009 ). Impacts of inbreeding on bumblebee colony fitness under field conditions . BMC Evolutionary Biology , 9 ( 1 ), 152 . doi: 10.1186/1471-2148-9-152 OpenUrl CrossRef PubMed ↵ Whiting , P. W. ( 1943 ). Multiple alleles in complementary sex determination of Habrobracon . Genetics , 28 ( 5 ), 365 – 382 . doi: 10.1093/genetics/28.5.365 OpenUrl FREE Full Text ↵ Wilfert , L. , Gadau , J. , & Schmid-Hempel , P. ( 2006 ). A core linkage map of the bumblebee Bombus terrestris . Genome , 49 ( 10 ), 1215 – 1226 . doi: 10.1139/g06-075 OpenUrl CrossRef PubMed ↵ Williams , P. , H., & Osborne , J. , L. ( 2009 ). Bumblebee vulnerability and conservation world-wide . Apidologie , 40 ( 3 ), 367 - 387 . Retrieved from doi: 10.1051/apido/2009025 OpenUrl CrossRef ↵ Woyke , J. ( 1963 ). What happens to diploid drone larvae in a honeybee colony . Journal of Apicultural Research , 2 ( 2 ), 73 – 75 . doi: 10.1080/00218839.1963.11100063 OpenUrl CrossRef ↵ Woyke , J. ( 1980 ). Effect of sex allele homo-heterozygosity on honeybee colony populations and on their honey production I. Favourable development conditions and unrestricted queens . Journal of Apicultural Research , 19 ( 1 ), 51 – 63 . doi: 10.1080/00218839.1980.11099997 OpenUrl CrossRef ↵ Zareba , J. , Blazej , P. , Laszkiewicz , A. , Sniezewski , L. , Majkowski , M. , Janik , S. , & Cebrat , M. ( 2017 ). Uneven distribution of complementary sex determiner (csd) alleles in Apis mellifera population . Sci Rep , 7 ( 1 ), 2317 . doi: 10.1038/s41598-017-02629-9 OpenUrl CrossRef PubMed ↵ Zayed , A. , & Packer , L. ( 2005 ). Complementary sex determination substantially increases extinction proneness of haplodiploid populations . Proceedings of the National Academy of Sciences , 102 ( 30 ), 10742 – 10746 . doi: 10.1073/pnas.0502271102 OpenUrl Abstract / FREE Full Text ↵ Zou , Y. , Geuverink , E. , Beukeboom , L. W. , Verhulst , E. C. , & van de Zande , L. ( 2020 ). A chimeric gene paternally instructs female sex determination in the haplodiploid wasp Nasonia . Science , 370 ( 6520 ), 1115 – 1118 . doi: 10.1126/science.abb8949 OpenUrl Abstract / FREE Full Text View the discussion thread. Back to top Previous Next Posted August 15, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications Kelley Leung , Frank Becker , Peter Šima , Roland Kreskóci , Elzemiek Geuverink , Thomas V. M. Groot , Bart A. Pannebakker , Joost van den Heuvel bioRxiv 2025.08.12.669960; doi: https://doi.org/10.1101/2025.08.12.669960 Share This Article: Copy Citation Tools Bombus terrestris Complementary Sex Determiner (BtCSD) is identified as a conserved hymenopteran sex determination region: evolution, breeding, and conservation implications Kelley Leung , Frank Becker , Peter Šima , Roland Kreskóci , Elzemiek Geuverink , Thomas V. M. Groot , Bart A. Pannebakker , Joost van den Heuvel bioRxiv 2025.08.12.669960; doi: https://doi.org/10.1101/2025.08.12.669960 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13895) Bioinformatics (41953) Biophysics (21456) Cancer Biology (18595) Cell Biology (25521) Clinical Trials (138) Developmental Biology (13381) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24323) Genetics (15612) Genomics (22511) Immunology (17738) Microbiology (40401) Molecular Biology (17184) Neuroscience (88623) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7644) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)