dataesr/scientific-paragraphs-categorization huggingface dataset:
In case the distance between two recordings taken on the same day of the same call type was smaller than 500 m, we used the recording with the largest number of individuals of the analysed call type. To account for geographic variation in the availability of tree species and test for differences in food resources, we divided the study area into five regions with similar tree species and tested these regions separately (see Supplementary material Appendix 1 Table  S2  for more details about these regions). We used data from the seasons 2013/14-2017/18 and of all call types with at least 150 decisions to feed on a specific tree species in the data set. Tree species with less than ten decisions of red crossbills to feed their seeds were removed. To investigate what influenced the composition of red crossbill call types to feed on seeds of a specific tree species, we fitted a generalized linear mixed model (GLMM; Baayen 2008) with binomial error structure and logit link function  (McCullagh and Nelder 1989)  in R (<www.
To control for potential differences among seasons with regard to the overall frequency of occurrence of a given call type, we included recording season as an additional random effect. Furthermore, we included the interaction between tree species, site and season as a random effect ('observation level random effect') to test the possibility that the preferences of certain call types for certain tree species vary among sites and/or seasons. As an overall test for the influence of tree species and site, we used a permutation test  (Adams and Anthony 1996, Bolker 2008). To this end, we randomized the response within seasons and used the sum of the estimated contribution of the two random effects (standard deviations estimated for tree species and site) as a test statistic. We conducted 1000 permutations into which we included the original data as one permutation and determined the p-value as the proportion of permuted data sets revealing a test statistic at least as large as the original data.
Localization of call types' core breeding areas Differentiation following geographical isolation should be recognizable in a spatial analysis of call-type occurrence data. We evaluated our data (see Recordings and classification) in a grid with 2° of resolution in northern-southern and eastern-western directions within the World Geodetic System 1984 (WGS84). Within each grid cell, we counted the number of individuals of each given call type (18 call types of the red crossbill described in  Martin et al. 2019  and four additional call types (Supplementary material Appendix 1 Section 10. 3)) that were present, as well as the total number of red crossbills recorded and calculated the quotient of both. For each grid cell, the number of crossbill seasons with available sound recordings was noted, as a high proportion of a call type during several seasons is more indicative for a core breeding area than data from a single season (due to irruptions, call types might be outnumbered by other call types in individual seasons).
Our data are well distributed over Europe's western half, fragmentary in eastern Europe, and there are only some data points in the eastern Palearctic (Supplementary material Appendix 1 Fig. S2 ). To avoid multiple sampling of the same individual, we only used recordings which were separated from all others of the same call type by at least two kilometres or, in case they were made within smaller distances, when more than 100 d elapsed between them (red crossbills take about 39-47 d from nest building to fledging (Glutz von Blotzheim and Bauer 2001)). If there were several recordings of the same call group within this geographical distance (2 km) and time span (100 d), we used the one with the largest number of individuals (same method as in  Martin et al. 2019). A core breeding area was defined as a single limited area, where a selected call type was present with a proportion of more than 25% of the red crossbills for more than one season.
Relationship of conifers' seed production within the presumed core breeding area and the occurrence of the respective call type in the invasion area We analysed the relation between fructification within the presumed call types' core breeding area and the occurrence in the presumed invasion area of the respective call type. We collected conifer fructification data in Europe (Supplementary material Appendix 1 Table  S5  for the sources) and compared these with call-type distribution data determined as described in 'Localization of call types' core breeding areas'. Red crossbill occurrence data from May, June and July were excluded as this is the typical period for movements of red crossbills and it is not possible to know whether birds are leaving an area and therefore should be counted for the previous season or whether they are arriving and should be counted for the next season.
The fructification extent was divided into the four categories (Supplementary material Appendix 1 Table  S6 ) used by the ICP (International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests; <www. icp-forests. net>). We only used data of trees older than 35 years, as younger trees rarely have a good cone crop  (Philipson 1990 ) and are rarely utilized by red crossbills for foraging (pers. obs. ). To obtain a polygon of the core breeding and invasion area of each call type, a Kernel-Analysis with the package 'adehabitatHR'  (Calenge 2006)  for R was used. Within the presumed core breeding area (data within the 50% kernel), we searched for available data of conifer fructification. Fructification data were only available for the presumed core breeding areas of three call types (N1, N7, N11) for more than two seasons.
We compared these fructification data with the occurrence of the respective call type in its invasion area (75-95% Kernel). To account for varying recording effort, we did not use total numbers of recorded birds but used the ratio of the selected call type of all call types within its invasion area in a given season. We tested for the relation between fructification within the core breeding area and the occurrence within the invasion area (Supplementary material Appendix 1 Fig. S3 ) using a generalized linear mixed model (GLMM) with a binomial error structure and logit link function. We included fructification and call type (factor with three levels) as fixed effects in the model. An observation-level random effect (i. e. the combination of season, call type and area) and season were included as random effects.
For each season, we calculated the fructification index for each core breeding area of the three tested call types. The model was implemented in R (ver. 3. 3. 2, <www. rproject. org>) using the function 'glmer' of the package 'lme4'  (Bates et al. 2015). We determined the R 2 -value using the R-package 'MuMIn'  (Barton 2016) , following the method of  Nakagawa and Schielzeth (2013). We checked the model for overdispersion (dispersion parameter 0. 43) and also for the distribution of the best linear unbiased predictors (BLUPs; Baayen 2008); both were not an issue. Model stability was assessed by comparing the estimates obtained from a model based on all data with those obtained from models with the levels of the random effect 'season' excluded one at a time. This revealed the model to be stable (Supplementary material Appendix 1 Table  S7 ).
S2  for information about these areas and the occurring conifer species. We had no zeros in the dataset, however, as data were collected opportunistically at different sites, missing data also hint to an absence of any call types in the respective conifer species. An empty field means, the conifer species was not available at the recording sites. South of 44°N and west of 60°E, 15 coniferous tree species visited by red crossbills for foraging grow on a regular basis (Picea omorika occurs in a very small area, P. abies, P. orientalis, Pseudotsuga menziesii and Larix decidua are just spread in small areas as well; additionally, there are Pinus brutia, P. halepensis, P. mugo, P. nigra, P. peuce, P. pinaster, P. pinea and P. sylvestris  (Andersson 2005 , Farjon and Filer 2013 , San-Miguel-Ayanz et al.
Acknowledgements -We thank everyone who recorded crossbills for us, sent us recordings or uploaded them to online platforms (Supplementary material Appendix 1 Section 10. 6). Without their extraordinary help and effort, it would have been impossible to realize this project. We also thank  ICP  ( International Co-operative Programme on Assessment and Monitoring of Air Pollution Effects on Forests ; <www. icp-forests. net>)and the institutions mentioned in Supplementary material Appendix 1 Table S9 for the fructification data of conifers they provided. We want to thank especially  Johannes Honold  for discussions and motivation in the beginning of the project,  Severin Hauenstein  for discussions about the methodical approach of the study, and  Dagmar Bednárová ,  Anton Burkhardt ,  David Jardine ,  Władysław Kantorowicz ,  Yann Kolbeinsson ,  Anita Nussbaumer ,  Örn Óskarsson  and  Ron Summers  for discussions and for helpful information.
(2015). Therefore, R S is an integrated reproduction R scaled by the total number of eggs possible to be produced during the imago stage ET E0. Note, that the value of E 0 can be used to calculate the number of eggs, however, it does not have any influence on the amount of energy stored in the reproduction buffer. According to equation (  22 ), when the reproduction buffer (energy) density threshold is reached ðE R ðtÞ ¼ E T Þ the integrated scaled reproduction R S is equal to 1 and pupation can begin. However, to use this value we have to know when it would be reached. The time needed to obtain critical value ½E j R can be estimated from the pupation dataset. The threshold for pupation is reached most likely when most of the larvae pupate (if we started an experiment with the larvae approximately at the same age and size).
Mass loss of PS To be sure that larvae eat PS, we measured the mass loss of PS chips during the experiment. Results of mass loss in each sample variant are presented in Fig. 2A  and in Table  S5. The differences in the accumulated PS mass loss among the four groups (PS, PS large, PS þ oat, and PS þ oat. þvit. ) were probably due to the differences in larval survival rates (see section Survival). This could be shown when comparing data after, e. g. , 50 days of the experiment (Table  S5 ). The highest number of larvae, 46 specimens, was alive in the PS þ oat. group, and the highest mass loss was noted e 51. 4%. In the PS þ oat. þ vit. group, containing 45 specimens, after 50 days we observed lower mass loss e 39%. The lowest mass loss was noted for PS large e 32. 6% for 39 specimens, and for PS e 16.
The length of these assemblies is in range of that predicted by a kmer analysis with the K-mer Analysis Toolkit (KAT) (Supplementary Figure  1 )  [60] , which predicted A. ervi at 142. 83Mbp and L. fabarum at 99. 26Mbp. However, the L. fabarum assembly is larger than the estimate from KAT; we suspect that this may be due to duplications in the assembly, and future work should address these duplications. These assembly lengths are also within previous estimates of 110-180Mbp for braconids, including A. ervi  [61, 62]  and are on par with those predicted in other hymenopteran genomes (Table  2 ). Both genomes were screened for potential contamination (Supplementary Figures  2  and 3 , Supplementary Table  6 , Additional files 1 and 2) based on BLAST  [63]  matches to host aphids and results of the program blobtools  [64] , which jointly examines GC content and sequencing depth.
In addition to identifying likely bacterial scaffolds (A. ervi: 35 scaffolds/ 106Kbp removed, no scaffolds removed from L. fabarum), blobtools revealed one outlier scaffold in L. fabarum with high coverage and low GC content (tig00001511, 10,205 bp, 11. 1% GC). A BLASTn search against the NCBI nt database matched this to the mitochondrial genome of Aphidius gifuensis. In this and other parasitoids, the mitochondrial genome has been shown to be highly enriched with AT repeats, with GC contents that are nearly as low as the 11. 1% found in this L. fabarum scaffold (13. 5-17. 5%) [65]. The assemblies are available in NCBI (PRJNA587428, SAMN13190903-4) and can be accessed via the BioInformatics Platform for Agroecosystem Arthropods (BIPAA, https://bipaa. genouest. org), which contains the full annotation reports, predicted genes, and can be searched via both keywords and BLAST.
We constructed linkage groups for L. fabarum using phased SNPs from the haploid sons of a single female wasp from a sexually reproducing population. This placed the 297 largest scaffolds (> 50% of the nucleotides, Supplementary Table  7 , Supplementary Figure  4 , Additional file 3) onto the expected six chromosomes  [52]. With this largely contiguous assembly, we identified stretches of syntenic sequence between the two genomes, with > 60 k links in alignments made by NUCmer  [66]  and > 350 large syntenic blocks that match the six L. fabarum chromosomes to 28 A. ervi scaffolds (Supplementary Figures  5  and 6 ). The Maker2 annotation pipeline predicted coding genes (CDS) in both genomes separately, and these were functionally annotated against the NCBI nr database  [67] , gene ontology (GO) terms  [68, 69] , and predictions for known protein motifs, signal peptides, and transmembrane domains (Supplementary Table  5 ).
In A. ervi there were 20,328 predicted genes comprising 24. 7Mbp, whereas in L. fabarum there were 15,203 genes across 21. 9Mbp (Table  1 ). Matches to the BUSCO (Benchmarking Universal Single-Copy Orthologs) genes assessed completeness against the Insecta database genes at both the nucleotide level (A. ervi: 94. 8%, L. fabarum: 76. 3%, Supplementary Table  4 ) and protein level in the predicted genes (A. ervi: 93. 7%, L. fabarum: 95. 9%). These protein level matches are close to those found in other assembled parasitoid genomes, which report between 96 and 99% total coverage of BUSCO genes  [32] [33] [34] [35] [36] [37]. In both species, there was also high transcriptomic support for the predicted genes (77. 8% in A. ervi and 88. 3% in L. fabarum). A survey of transposable Elements (TEs) identified a similar overall number of putative TE elements in the two assemblies (A.
Gene family expansions To examine gene families that may have undergone expansions in association with functional divergence and specialization, we identified groups of orthologous genes that have increased and decreased in size in the two genomes, relative to one another. We identified these species-specific gene-family expansions using the Orthologous MAtrix (OMA) standalone package  [74]. OMA predicted 8817 OMA groups (strict 1:1 orthologs) and 8578 Hierarchical Ortholog Groups (HOGs, Additional file 8). Putative gene-family expansions would be found in the predicted HOGs, because they are calculated to allow for > 1 member per species. Among these, there were more groups in which A. ervi possessed more genes than L. fabarum (865 groups with more genes in A. ervi, 223 with more in L. fabarum, Supplementary Figure  14 , Additional file 8).
Lysiphlebus fabarum Library quality was also checked with FastQC  [164]. Illumina reads were filtered using Trimmomatic ver 0. 33 to remove low quality sequences (Q < 25, 4 bp window), to trim all Illumina primers, and to discard any sequence shorter than 50 bp or without its mate-pair. NextClip was used to remove all improperly oriented read pairs. Raw PacBio reads were error-corrected using the quality filtered Illumina data with the program Proovread  [170]. These error-corrected reads were then used for de novo assembly in the program canu v1. 0  [171]. Since our PacBio reads were expected to have approximately 30X coverage (based on the presumed size of 128Mbp), canu was run with the recommended settings for low coverage data (corMhapSensitivity = high corMinCoverage = 2 errorRate = 0.
Genome completeness and synteny Completeness of the two assemblies was assessed by identifying Benchmarking Universal Single-Copy Orthologs (BUSCOs) using the BUSCO v3. 0. 2 pipeline in genome mode  [187]. We identified single copy orthologs based on the insecta_db9 (1066 genes, training species: Nasonia vitripennis). Synteny between the two genomes was assessed using the NUCmer aligner, which is part of the MUMmer v3. 23 package  [66]. For this, we used the L. fabarum chromosomes as the reference, and included the scaffolds not incorporated into chromosomes (total 1407 pieces). The A. ervi assembly was mapped to this using the default settings of NUCmer.
Predictive gene annotation For both assembled genomes, gene predictions were generated using MAKER2  [188]. Within MAKER2, predictive training was performed in a three step process. A first set of genes was predicted by similarity to known proteins or contigs from RNAseq in the same species (described below). This gene set was used thereafter for training both Augustus  [189]  and SNAP  [190] , in two steps, with the results of the first training re-used to train the software in the second round. Transcriptomic evidence was provided separately for each species. For A. ervi, six separate de novo transcriptome assemblies from Trinity  [191]  were constructed, one each for the adults reared on different hosts (NCBI PRJNA377544)  [160]. For each transcript, we only included variants based on filtering with RSEM v 1.
2. 21 using the option -fpkm_cutoff 1. 0, --isopct_cutoff = 15. 00. This resulted in 452,783 transcripts. For L. fabarum, we utilized a joint transcriptome, built using RNAseq data (NCBI PRJNA290156) collected from adults  [42]  and 4-5 day old larvae  [157]. Peptide evidence came from the Hymenoptera genomes database (http://hymenopteragenome. org, Acromyrmex echiniator v3. 8, Apis mellifera v3. 2, Nasonia vitripennis v1. 2), from the Bio-Informatics Platform of Agroecosystems Arthropod database (https://bipaa. genouest. org, Hyposoter didymator v1. 0), and Drosophila melanogaster (http://flybase. org, v6. 13), and SwissProt (October 2016) databases. Summary statistics were generated with GAG  [192]. Transcriptomic support for the predicted genes was estimated by mapping available transcriptomic data (same as above) to the respective genomes using STAR  [193]  in the "quantMode".
Functional annotation The putative functions of the proteins predicted by the above pipeline were identified based on BLASTp (v2. 5. 0) matches against Genbank nr (non-redundant GenBank CDS translations+PDB + SwissProt+PIR + PRF) release 12/2016 and interproscan v5 against Interpro (1. 21. 2017). GO terms associations were collected from BLAST nr and interproscan results with blast2GO (v2. 2). Finally, transmembrane domains were identified with Hidden Markov Models (HMM) in tmhmm v2. 0c, and peptide signals with signal (euk v4. 1)  [194, 195].
Transposable elements Transposable elements (TE) were predicted using the REPET pipeline  [196] , combining de novo and homology-based annotations. Repetitive elements were identified de novo across all scaffolds larger than the scaffold N50 for each species, following recommendafrom the developers of REPET for draft genome assemblies. Within these, repetitive elements were identified using a BLAST-based alignment of each genome to itself followed by clustering with Recon  [197] , Grouper  [198]  and Piler  [199]. For each cluster, a consensus sequence was generated by multiple alignment of all clustered elements with MAP  [200]. The resulting consensus was then scanned for conserved structural features or homology to nucleotide and amino acid sequences from known TEs (RepBase 20. 05)  [201, 202]  using BLASTER (tblastx, blastx)  [196]  or HMM profiles of repetitive elements (Pfam database 27.
GC content and codon usage We examined several measures of nucleotide composition, at both the nucleotide and protein level. Whole genome GC content was calculated by totaling the numbers of A, C, T, and G in the entire assembly. In the predicted coding sequences, this was also calculated separately for each predicted gene and third position GC composition was calculated separately in the predicted coding sequences. In all cases, this was done with the sscu package in R  [210]. Relative Synonymous Codon Usage (RSCU) was extracted from the entire CDS using the seqinR package in R  [211] , and visualized with a PCA (R packages factoextra, reshape, and ggplot2)  [212] [213] [214]. To examine GC content in coding genes of other insects, we downloaded the 118 available CDS in the RefSeq database of NCBI (date: October 2018) and again calculated per-gene GC content.
Orphan genes We identified orphan genes as those for which we could not find orthologs in any other sequenced genomes (Supplementary Table  11 ). To do this, we first used OrthoFinder  [216]  to generate clusters of orthologous and paralogous genes among the predicted genes (CDS) from the genomes of A. ervi and L. fabarum, as well as five other sequenced parasitoids (Diachasma alloeum, Fopius arisanus, Macrocentrus cingulum, Microplitis demolitor and Nasonia vitripennis). OrthoFinder produces a set of genes that were not assigned to any orthogroup. We identified species specific genes, which we are calling orphan genes, by removing all genes that had hits to any other genes in the nt, nr, and swissprot NCBI database (June 2019). Within these putative orphans, we only retained those with transcriptomic support.
Gene family expansions We examined gene families that have expanded and contracted in A. ervi and L. fabarum relative to one another using the OMA standalone package (v2. 2. 0, default values)  [74]. OMA was used to compute orthologs (OMA groups) and Hierarchical Orthologous Groups (HOGs) for the predicted proteins of L. fabarum (OGS1) and A. ervi (OGS3): 15,203 and 20,344, respectively. While OMA groups consist of strict 1:1 orthologs between OGS1 and OGS3, HOGs contain all orthologs and paralogs of a given predicted gene family. HOGs were parsed with a custom Perl script to identify all gene families in which one of the wasp species contained more members than the other. We focused on only the groups that contained more than 20 genes (10 groups, Supplementary Figure  15 ). These were identified by BLASTx against the nr database in NCBI.
Venom proteins The L. fabarum venom proteomic analysis was performed from 10 extracted glands (Supplementary Figure  17 ). The 16 most visible bands in 1D gel electrophoresis were cut, digested with trypsin and analyzed by mass spectrometry. All raw data files generated by mass spectrometry were processed to generate mgf files and searched against: (i) the L. fabarum proteome predicted from the genome (L. fabarum annotation v1. 0 proteins) and (ii) the L. fabarum de novo transcriptome  [42]  using the MASCOT software v2. 3  [217]. The mass spectrometry proteomics data are deposited in the ProteomeXchange Consortium (proteomecentral. proteomexchange. org) via the PRIDE partner repository  [62] , with the ID PXD015758. Sequence annotation was performed based on BLAST similarity searches. Signal peptide prediction was performed with SignalP  [194, 195].
Desaturases Desaturase genes in both genomes were automatically identified and annotated with GeMoMa  [223]  using desaturase gene annotations from Diachasma alloeum, Fopius arisanus, and Microplitis demolitor, retrieved from NCBI's protein database as queries (retrieved May 2017). Additionally, all desaturase genes were manually inspected. To determine the cuticular hydrocarbon (CHC) profiles in A. ervi, wasps were freeze-killed and stored separately by sex at -20 °C. For CHC extraction, single individuals were covered with 50 μl of MS pure hexane (UniSolv) in 2 ml GC vials (Agilent Technologies,) and swirled for 10 min on a Thermo-shaker (IKA KS 130 Basic, Staufen). The hexane extracts where then transferred to a fresh conical 250 μl GC insert (Agilent Technologies), where the hexane was completely evaporated under a constant flow of CO 2.