Subjects and sample collection

Our analyses used three males of six passerine species representing several families within the Passerida clade. Five species were free-living, while males from the captive population of zebra finches housed at the Institute of Vertebrate Biology bird breeding facility were used as the reference samples. Males of free-living species were released, immediately after which a sperm sample was obtained. Sampling was performed in May-June 2014-2015 in southern Bohemia, Czech Republic. Sperm samples were obtained by a standard massage of the cloacal protuberance (Wolfson 1952). Ejaculates (approx. 1 µl) were immediately dissolved in 20 µl of sterile phosphate-buffered saline (PBS). An aliquot (3 µl) was stored in 10% formalin for further inspection of sperm morphology and ejaculate quality (presence of red blood cells or debris), and the rest was deep frozen in liquid nitrogen. Only clean samples containing only sperm cells (based on the inspection of aliquotes) were used in further proteomic analysis and sperm measurements.

Sperm morphology

Formalin-fixed sperm samples were transferred with an automatic pipette onto a slide, smeared, air-dried and rinsed in dH₂O to wash off impurities. The images were captured on a digital camera (DP71, Olympus, Japan) mounted on a light microscope (BX51, Olympus, Japan). The measurements of three sperm parameters (head length, midpiece length and tail length) were made using the imaging software QuickPHOTO Industrial (Olympus, Japan). From these measurements, we calculated the total sperm length and flagellum length (the sum of midpiece length and tail length). For each individual, ten spermatozoa with normal sperm morphology were analysed, and the measured values were averaged for each species.

Proteomic analysis

Sperm samples (n = 18) were defrosted at room temperature, vortexed and a total of 10 µl transferred to new 0.5 ml tubes to which we added 20 µl of PBS. The samples were vortexed (1 min) and centrifuged (14,000 RCF) to obtain sperm cells (pellets: n = 18) and seminal fluids (supernatant: n = 18). Next, we added 100 µl of isolation buffer containing 8 M urea, thiourea, 1% Chaps, 20 mM DTT and dH₂O to each tube and left it at room temperature for one hour. Protein samples were precipitated with the ice-cold acetone (1:4) and centrifuged (14,000 RCF). This procedure was followed by a re-suspension of dried pellets in the digestion buffer (1% SDC, 100 mM TEAB – pH = 8.5). Protein concentration of each lysate was determined using the BCA assay kit (Fisher Scientific, Waltham, MA, USA). Cysteines in 20µg of proteins were reduced with a final concentration of 5 mM TCEP (60 °C for 60 min) and blocked with 10 mM MMTS (i.e. S-methyl methanethiosulfonate, 10 min at room temperature). Samples (n = 36) were cleaved with trypsin (i.e. 1/50, trypsin/protein) at 37 °C overnight. Peptides were desalted on a Michrom C18 column. Eluting peptide cations were converted to gas-phase ions by electrospray ionisation and analysed on a Thermo Orbitrap Fusion (Q-OT-qIT; Thermo Fisher, Waltham, MA, USA) using the same conditions and setup as previously described (Stopka et al. 2016, Kuntova et al. 2018). For example, survey scans of peptide precursors from 400 to 1,600 m/z were performed at 120 K resolution (at 200 m/z) with a 5 × 105 ion count target. Tandem MS was performed by isolation at 1.5 Th with the quadrupole, HCD fragmentation with normalised collision energy of 30, and rapid scan MS analysis in the ion trap. The MS2 ion count target was set to 104 and the max injection time was 35 ms. Only those precursors with charge states 2-6 were sampled for MS2. The dynamic exclusion duration was set to 45 s with a 10 ppm tolerance around the selected precursor and its isotopes. Monoisotopic precursor selection was turned on. The instrument was run in top speed mode with 2 s cycles.

Fig. 2.

Diverging patterns of passerine proteomes. Sparse Partial Least Squares Discriminant Analysis (sPLS-DA) revealed that all the studied passerines significantly separated from TG. Sperm proteomes are in A, B, seminal fluid proteomes are in C, D. Note that the phylogeny of the two major clades is well demonstrated by the second component (Comp2 – y-axis) in all discrimination. AP – Acrocephalus palustris; CC – Cinclus cinclus; CH – Chloris chloris; HR – Hirundo rustica; TG – Taeniopygia guttata.

Bioinformatics

All data were quantified using MaxQuant software version 1.6.34. The false discovery rate (FDR) was set to 1% for identifying all peptides and proteins. We set a minimum peptide length of seven amino acids. Enzyme specificity was set as C-terminal to Arg and Lys, allowing cleavage at proline bonds 98 and a maximum of two missed cleavages. Dithiomethylation of cysteine was selected as fixed modification, and N-terminal protein acetylation and methionine oxidation as variable modifications. The ‘match between runs’ feature of MaxQuant was used to transfer identifications to other LC-MS/MS runs based on their masses and retention time (maximum deviation 0.7 min). Quantifications were performed using the label-free algorithms 65 with a combination of unique and razor peptides. The Andromeda search engine was used for the MS/MS spectra search against the Uniprot database (UP000007754, 2021), with all duplicates being removed. Our raw dataset was LFQ normalised (Cox et al. 2014) with MaxQuant and is available as Table S1. From this point, all data manipulations, filtering, statistical analyses and plots were generated in R software (Crawley 2007). First, we removed all protein identifications that were detected but not quantified (zero abundances), all proteins with less than two unique peptides, and all contaminants from MaxQuant output files (porcine trypsin, bovine albumins, human keratins, etc.). Next, we removed singletons (rare cases) from our data such that only those proteins produced by all three individuals in at least one species (e.g. TG) were passed to further analyses; otherwise, the whole row was deleted. This data format served for data exploration with the Sparse Partial Least Squares Discriminant Analysis (sPLS-DA) within the ‘mixOmics’ package (Rohart et al. 2017). To make sure that interspecific variation is higher than the variation between replicates, we performed Spearman’s rank correlations between individual replicates and between averaged values per species. Spearman’s rank correlation was also used to find relationships between midpiece length and the abundance of individual proteins. Gene ontology (GO) searches were performed using Gene Set Enrichment Analysis (GSEA) within the ‘clusterProfiler’ package (Yu et al. 2012) from Bioconductor, while significant GO terms were searched by using p-value cut-off P = 0.1 (P-value < 0.05 revealed only a few terms) within the latest chicken database (August 2022) – org.Gg.eg.db. We recognise that there are limitations to using the chicken database for passerine species. If we used a passerine database, we would struggle with another limitation: the presence of zebra finches in our samples and pure annotations. So, we opted for chicken, which is similarly phylogenetically distant to all our studied species. Also, this approach does not look at species-specific proteins due to data reduction. To add, we identified 97 significantly correlated proteins in sperm samples, which successfully mapped to 72 GO terms (i.e. 74%). Our final dataset has an approximately similar number of proteins in all the species: 1,198 and a reasonable number of unique peptides (> 2). With this approach, we probably omitted many reproductive proteins that evolve faster, but the data set thus enabled comparison between species, at least in the proteins they share. All plots and figures were generated in R using ‘ggplot2’ and related packages (Wickham 2016).

Characterisation of sperm and seminal fluid proteomes

The proteomic dataset was generated with nLC-MS/ MS from 36 samples containing 18 sperm samples and 18 fluid samples, obtained from six species with three replicates per species. The whole matrix extracted from label-free quantification (Cox et al. 2014) was LFQ normalised. The resulting dataset is based on 144,449 high-quality spectral matches to 2,384 protein identifications (Table S1). Of these, 88.2% were identified by multiple unique peptides shared between species (on average 7.2). These searches were conducted against T. guttata proteome (UP000007754, 2021) in the absence of the annotated genomes in the other species. We removed all rows with fewer unique peptides than two, all potential contaminants and all proteins identified but not quantified. Next, the matrix was split into the sperm matrix (n = 1,203) and seminal fluid matrix (n = 1,130), while both matrices contained a similar number of identifications. After the data reduction, we used ‘normalizerDE’ in R to perform all eight most common normalisations (e.g. log2, VSN, Quantile, GI, etc.) and used the Quantile normalisation simply because the within triplicate variation was the lowest. For the subsequent analyses, we used two datasets, the complete proteome and the so-called core proteome, which used only the proteins detected in all the studied samples/individuals.

The hierarchical clustering of sperm and seminal fluid proteomes revealed a surprisingly similar branching structure (entanglement = 0.01). This structure is, however, different from the phylogenetic branching in Fig. 1A. For example, genetically different TG and HR cluster together in Fig. 1B but not in Fig. 1A. Thus, the phylogeny was to some extent overridden by as yet unspecified ecological factors that have driven proteome evolution. Correlations between individual replicates were high (r > 0.8). In contrast, correlations between averaged TG replicates with other species were much lower (approx. 0.5). Particularly, Spearman's rank correlation between TG and other passerines was relatively low (SPERM: r_CH = 0.51, r_CC = 0.29, r_AP = 0.46, r_PC = 0.48, r_HR = 0.52; FLUIDS: r_CH = 0.51, r_CC = 0.41, r_AP = 0.46, r_PC = 0.44, r_HR = 0.46; P < 2.2e-16), while the correlation between biological replicates was generally high (SPERM: TG – 0.77-0.85, CH – 0.77-0.86, CC – 0.75-0.85, AP – 0.79-0.86, PC – 0.79-0.83, HR – 0.86-0.89; FLUIDS: TG – 0.78 for all comparisons, CH – 0.73-0.88, CC – 0.56-0.71, AP – 0.80-0.85, PC – 0.84-0.88, HR – 0.85-0.89; P < 2.2e-16). High correlations between biological replicates are also demonstrated in a tanglegram (Fig. 1B), where individuals cluster together within a particular species.

Sources of proteomic variation

We used the Sparse Partial Least Squares Discriminant Analysis (sPLS-DA) to assess the potential sources of variation in proteomes. In Fig. 2, we demonstrate using the Area Under Curve analysis (AUC) that all the five passerine species in this study have patterns significantly diverging from TG (e.g. for Fig. 2A – Comp1: AUC = 1, P = 0.008; Comp2: AUC = 1.0, P = 0.008). A similarly diverging pattern observed in sperm samples (Fig. 2A, B) was obtained with sPLS-DA on the data from seminal fluids (Fig. 2C, D) (P-values for all comparisons are in Table S1). Note that sPLS-DA is separated by the second component (Comp2), the two major clades from Fig. 1A. TG, CH, and CC are well separated from AP, HR, and PC, reflecting species phylogenetic relatedness. Next, we investigated the key proteins responsible for separating given species and whether these proteins represent biologically relevant features such as particular sperm metabolic pathways or molecular functions. For this purpose, we visualised the top 15 loading values and thus the importance of proteins contributing to each principal component (Comp1 and Comp2) of the sPLS-DA model (Fig. 3). Proteins are identified as important based on the mean abundancy value in each species. Positivity and negativity of importance values (i.e. positive and negative contribution on each component) indicate species in which the protein is highly abundant; however, proteins with positive and negative importance values contribute equally to species separation. The most important result from this analysis revealed that sperm and seminal fluids have different proteins responsible for the differentiation, which thus provides evidence that the two fractions from the ejaculate samples were not biased by extensive cross-contamination. There is yet another message in the visualisations of loadings values; some species differ from others by lower protein expressions while others have higher expressions when compared to other species (do not mistake with negative and positive expression indicated by negative and positive values of importance, respectively). The same colour code representing species is used as in Fig. 2.

Fig. 3.

The contribution of particular proteins to species differentiation is shown with loadings from sPLS-DA analysis. The colour coding is similar to in Fig. 2 and indicates the species in which the protein has the maximum mean value. Sperm proteomes are in A and B, while seminal fluid proteomes are in C and D (complete proteomes left, core proteomes right). The importance of proteins for separation between species is ranked from the bottom (the greatest) to the top. AP – Acrocephalus palustris; CC – Cinclus cinclus; CH – Chloris chloris; HR – Hirundo rustica; TG – Taeniopygia guttata.

Fig. 4.

Correlation of complete proteomes with sperm midpiece length. We performed Spearman's rank correlation between midpiece length and sperm proteomes (A) and fluids proteomes (B) averaged over species. Significantly correlated proteins are scaled from green (P < 0.05) to blue (P < 0.01).

Highly important sperm proteins play important roles in sperm metabolism, e.g. ATP-dependent 6-phosphofructokinase (PFKP), phosphoglycerate mutase 1 (PGAM1), aldehyde dehydrogenase 6 family member A1 (ALDH6A1), glycerol-3-phosphate dehydrogenase NAD(+) (GPD1L), glycerol-3phosphate dehydrogenase (GPD2), and in energetics, e.g. cytochrome c, somatic (CYCS), NADHubiquinone oxidoreductase subunit, mitochondrial (NDUFS1), NADH dehydrogenase iron-sulfur protein 2, mitochondrial (NDUFS2), cytochrome b-c1 complex subunit Rieske, mitochondrial (UQCRC1), cytochrome b-c1 complex subunit 2, mitochondrial (UQCRC2) and complex III subunit 9 (UQCR10) (Fig. 3A, B). Like spermatozoa, seminal fluids contain proteins with a high importance involved in energetics, e.g. CYCS, NDUFS1, metabolism, e.g. aconitate hydratase, mitochondrial (ACO2), acetyltransferase component of pyruvate dehydrogenase complex (DLAT), PFKP, fructose-bisphosphate aldolase (ALDOC), GPD1L, fertilisation, e.g. tetraspanin CD9, and in protein degradation, e.g. 26S proteasome nonATPase regulatory subunit 1 (PSMD1) (Fig. 3C, D).

Fig. 5.

Enriched gene ontology (GO) terms revealed with Gene Set Enrichment Analysis (GSEA). The colour of the dots represents the level of significance adjusted with Benjamini-Hochberg correction for multiple testing. The size of the dots indicates the number of proteins enriched in that term. Activated proteins (left) are those positively correlated with the midpiece length. Note that many of them are related to mitochondrial metabolism. Suppressed proteins (negatively correlated) are on the right.

Fig. 6.

GSEA analysis showing the prevalence of GO terms associated with mitochondria. The distribution of activated and suppressed proteins of the complete sperm (A) and seminal fluid (B) proteomes and their association with one or more GO pathways (see the intersection of black dots) is captured with UpSet plots. Network analysis in C shows the complex association between proteins and multiple annotation terms.

In detail, PFKP has positive importance values for both sperm datasets and the core proteome of seminal fluids in PC, a highly promiscuous species (EPP = 48% in (Lifjeld et al. 2019)). In the same species, CYCS is important for separating the core sperm proteome. On the other hand, proteins DLAT and heat shock protein (HSPD1), which are both mitochondrial, are positively separating seminal fluids of HR from other species regardless of the dataset used. Further, a similar contribution is observed for protein CYCS in differentiation of AP species. In CC, monogamous species with the shortest sperm tails examined is cytosolic GPD1L responsible for separation in complete proteomes, although with negative values of importance. Taken together, although the sperm design looks pretty much the same in all the six studied species, there has been a selective pressure that has driven biochemical sperm components to species-specific states.

Interestingly, mammalian acrosomes also contain mitochondrial and proteasomal proteins (Guyonnet et al. 2012, Otčenášková et al. 2023), presumably because redundant mitochondria are removed from the cytosol after midpiece formation, and the resulting proteins are ubiquitinated and deposited in acrosomes. Here, we provide new evidence that this could be also true for avian species. Thus, the presence of mitochondrial proteins in seminal fluids practically means that the similarity between sperm and fluid dendrograms in Fig. 1B is highly influenced by proteins occurring in both fractions. On the other side, the presence of acrosomal proteins, e.g. tetraspanin CD9, acrosin-binding protein (ACRBP) and zona pellucida-binding protein 2 (ZPBP2), in seminal fluids might provide important information on the acrosomal content. Noteworthy, these proteins, together with acrosin (ACR), tetraspanin CD81, hyaluronidase PH-20 (SPAM1) involved in sperm-oocyte interactions, have been recently detected in seminal fluid proteomes in chicken and two closely-related sparrow species (Borziak et al. 2016, Li et al. 2020, Rowe et al. 2020). Although the abundance of acrosomal proteins in seminal fluid is often explained by acrosome damage either during sample preparation (which might also be our case) or as an indicator of poor semen quality (Li et al. 2020), they still might be of biological significance.

Broader context

Taking a broader view, we touched upon a fundamental question in reproductive biology and evolution: the relationships between the phylogeny, varying levels of sperm competition, sperm morphology (size), and proteomic variation. We found that each species differs from the other by species-biased groups of proteins and that these groups represent biologically relevant features such as mitochondrial energy production, metabolism, protein degradation and sperm-egg fusion. These features represent important building blocks of sperm cells in the view of sperm competition because mitochondria are the drivers of proper sperm movement while the acrosome-mediated gamete interaction and fusion are the key processes during fertilisation. Particular roles of these proteins were corroborated by gene ontology analysis, where most proteins that correlate with sperm midpiece length are involved in mitochondrial energetic metabolism. Although this result appears expected, it is not trivial. This is because some mammalian species (e.g. muroid rodents) have variable apical hooks containing proteins involved in fertility and acrosome reaction (Breed 2004, Immler et al. 2007, Johnson et al. 2007, Clift et al. 2009, Sandera et al. 2013, Hook et al. 2021) and, for example, resulting proteomes of mouse spermatozoa and acrosomes contain proteins enriched in energetic metabolism (Vicens et al. 2017, Otčenášková et al. 2023). In this study, the fertilisation process was not detected by GO analysis. However, it was detected with sPLS-DA techniques, which revealed that combining several techniques to extract relevant biological features is helpful. Although this study brings just facets of sperm divergence on the proteomic level, the results of fifty years of solving the enigma of sperm competition in animals (Parker 1970, 1992, Dixson 1993, Dixson & Anderson 2004, Gomendio et al. 2006, Immler 2008, Pizzari & Parker 2009, Vicens et al. 2017, Birkhead & Montgomerie 2020) are now explored with new OMICs techniques, which helps to detect important candidate proteins and their involvement in particular biological processes. These could be further explored on larger datasets to reconstruct the selective forces that shaped the evolution of sperm competition in animals.