1. IntroductionAbout 200,000 species of venomous animals belonging to many different taxa (e.g., cnidarians, arthropods, reptiles, fishes, mammals) are currently known to science [1]. Among them, snakes are arguably the most notorious ones. Of more than 3900 snake species currently recognised [2], about 300, mainly belonging to the families Viperidae (e.g., adders, rattlesnakes, palm pit vipers), Elapidae (e.g., cobras, coral snakes, sea snakes), and Colubridae (e.g., vine snakes, twig snakes, boomslangs) are considered of medical importance by the World Health Organization (WHO) [2,3,4].Snake venom is a complex mixture of peptides, proteins, small organic molecules, and salts [5,6], able to disrupt the general homeostasis of the envenomated organism, affecting it in different ways and with different levels of specificity and potency [7,8,9]. Snake venom has long been and remains in the spotlight for researchers from all over the world, mostly because of the medical importance of snakebites for human health [10,11,12]. For many years, it has been studied through traditional biochemical and pharmacological approaches, often focusing mainly on abundant toxins present in venoms produced by the most common medically relevant species (e.g., Daboia russelii [13], Bothrops jararaca [14]), and neglecting the venoms of a large number of rarer and/or generally less studied species [15,16,17].Advances in the so-called “omic” technologies, defined as the application of high-throughput methodologies [18], and their inclusion in the conventional analysis procedures, completely revolutionised snake venom studies. The term “venomics” currently describes the comprehensive study and characterisation of the whole venom profile of a toxic organism by the means of integrated “omic” methodologies, namely proteomics, transcriptomics, and genomics [5,18,19]. Specifically, modern snake venomics allow for the identification of venom components both directly, through the analysis of the protein content of crude venom (i.e., proteomics), and indirectly, through the sequencing of venom gland mRNA or cDNA (i.e., transcriptomics), or of the full genome (i.e., genomics) of the analysed species [19]. By allowing the rapid characterisation of all venom components of a growing number of snake species from all over the world, the deployment of snake venomics is gradually filling the knowledge gap left by the application of previous, less comprehensive approaches. The remarkable amount of information produced by venomics about composition and properties of different snake venom phenotypes, together with its implementation in functional studies, is helping to elucidate the processes underlying snake venom origin and evolution, and is driving the development of powerful therapeutic tools to be used to mitigate the global burden of snakebite [20,21] and successfully treat different types of diseases (e.g., [22,23] and references within).Despite the growing attention snake venom studies have received since the rise of venomic approaches, it has been noted that, in this line of research, some topics tend to be investigated more frequently than others [5,24,25]. Questions regarding snake ecology (e.g., interactions between snakes and their prey), for example, are seldom addressed in venomic studies, despite being relevant for both evolutionary biology and the development of effective snakebite mitigation strategies [15,26,27]. Similarly, research efforts seem to be greatly biased towards species belonging to the families Elapidae and Viperidae, whereas other snake families (e.g., Psammophiidae, Pseudoxyrophiidae, Atractaspididae, Homalopsidae) tend to be overlooked [15]. Although these tendencies in snake venom studies have long been recognised (e.g., [15,28]), they have never been formally quantified.

In the present work, we aim to provide formal categorisation and quantification of the current biases in snake venom research. To this end, we (i) present an estimation and description of the prevailing trends in snake venom studies published between 1964 and 2021, (ii) analyse whether and how the focus of the retrieved studies changed in terms of topics and taxa covered across the defined time frame, and (iii) test whether potential biases in terms of number of articles dedicated to each retrieved snake species could be related to specific factors (e.g., taxonomy, biogeographic realm of origin). We expect our findings to uncover the taxonomic and topic imbalances present in this field of study, and potentially help identify their origin and define the directions to follow to redress them.

4. ConclusionsWe acknowledge the possible presence of methodological limitations in this work. Specifically, the exclusive use of Google Scholar for article search, in combination with the article selection criteria applied, likely excluded some potentially relevant publications from our analysis. Nevertheless, our results are in line with trends and biases in snake venom studies already reported in the previous literature [15,28,38], and are thus to be considered reliable. We detected an overall positive trend, with a consistent increase since the early 2000s in the number of snake venom studies published yearly, even more evident over the last decade. Nevertheless, our analysis also highlighted a consistent neglect of snake families of supposedly minor medical relevance (e.g., Homalopsidae, Psammophiidae, Pseudoxyrhophiidae), an apparently limited focus on some of the areas most impacted by snakebite (i.e., Asia and Africa), and potentially minor interest in the ecological and functional context of snake venom. The study and characterisation of the venoms produced by other venomous snake taxa, excluding the typically more studied families (i.e., Colubridae, Elapidae, Viperidae), should be implemented in future snake venom studies in order to increase the knowledge about snake venom evolution and composition, and help widen the spectrum of treatable snakebite envenomations. Additionally, more effort should be put into developing studies focusing on species originating from the areas where snakebite incidence is high and the economical level low, such as the Indomalayan and the Afrotropic realms, which appear to be underinvestigated. Finally, snake venom should be analysed taking into account ecological and functional contexts of the species producing it, in order to pave the way to obtaining a more detailed and comprehensive view of the driving forces behind snake venom evolution and variation. We hope that by providing a qualitative and quantitative estimate of the taxonomic, geographic, and topic biases present in snake venom research, our work will be useful to define a road map for future efforts aiming at focusing on the most glaring knowledge gaps in this field of study. 5. Materials and Methods 5.1. Article SelectionPublications considered for the current study were gathered and organised using the Google Scholar (Google Inc. (Menlo Park, CA, USA)) web search engine (https://scholar.google.com), between the months of December 2018 and March 2022. To perform the search, the following query was used, applying every possible combination of the ten selected keywords: ( OR ) AND ( OR OR ) AND ( OR OR OR OR ). The evaluation timeframe we defined went from 1964 to 2021. Search results were sorted by relevance following Google Scholar default search options, with the quality of the result-search match being higher on top of the result list and progressively decreasing. We thus reviewed for consideration the first 200 articles obtained for each keyword combination searched, checking their suitability for inclusion in the final dataset. Articles focusing on the study of snake venom composition and variation, presenting either a protein-centred venom approach or an indirect approach based on different techniques (e.g., transcriptomics, bioinformatics, toxicity assays) were taken into account for analysis. Articles not investigating whole snake venoms (e.g., reviews, publications focusing only on single venom fraction analysis, single toxin studies), and/or not published in refereed, impacted journals were not considered.

The following information was recorded from each article: (i) publication year, (ii) taxonomy of the analysed species, (iii) country and biogeographic realm of origin of the analysed specimens, and (iv) topics covered.

5.2. Taxonomic InformationIn order to assess what the most studied and most represented snake taxa were, information about family, subfamily, genus, and species of the specimens analysed in each article was collected. Due to phylogenetic uncertainty within the family Elapidae [55,56,57], we did not consider subfamilies for this group, but instead divided it into two main categories widely used in the literature, irrespective of concerns over monophyly [58,59,60]: (i) Old World and American elapids and (ii) Australo-Papuan and marine elapids. The retrieved taxonomic information was updated mainly following the taxonomy reported by The Reptile Database [2], based on information about species names and sampling localities of the specimens. When insufficient locality and taxonomic information did not allow the unambiguous identification of the analysed snake species, we kept the specific IDs as reported in the original articles. 5.3. Hazard CategoriesIn order to test whether the harmful potential of a species’ venom could influence eventual biases in terms on number of studies dedicated to it, we developed a hazard index based on the existing bibliography (e.g., [11,40]), WHO guidelines (e.g., [4,33,34]), and authors’ opinion. We classified the snake species considered in the retrieved studies into four categories, based on the severity of the envenomation they can cause: (i) category 1—“critical clinical relevance”: envenomations have a generally high chance to cause death or significant disability if professional medical care is not obtained; (ii) category 2—“high clinical relevance”: envenomations usually cause significant illness, hospitalisation is required, death and/or disability are unlikely but possible if professional medical care is not obtained; (iii) category 3—“moderate clinical relevance”: envenomations are unpleasant but typically not life-threatening, significant disability is exceptional, typically treated symptomatically; (iv) category 4—“low clinical relevance”: envenomations likely cause only very mild symptoms (e.g., local swelling, itching, limited blistering), generally not interfering with normal activities and not being life-threatening, and professional medical care rarely necessary. Species we could not assign to any of the abovementioned categories were classified as “unknown” and not included in the analyses. 5.4. Origin of the SpecimensInformation about the country where each snake species that produced the analysed venom samples came from, and the corresponding biogeographic realm, was also gathered and used to assess possible geographical biases in snake venom studies. Country and biogeographic realm of origin of specimens for which information about the place of origin was ambiguous or unavailable (e.g., captive specimens, pooled venoms) were considered as “unknown” and not included in the analyses. Biogeographic realms were identified following the RESOLVE Ecoregions 2017 website [61]. 5.5. Topics Covered

In order to identify the most investigated research topics in the retrieved articles, we gathered information about the research topics covered in the reviewed publications, and grouped them into eight categories: (i) “venom characterisation”: defining the composition of the venom of snake species through the application of one or more techniques, from basic venom fractionation to “omic” approaches (i.e., proteomics, transcriptomics, genomics); (ii) “antivenomics and neutralisation”: evaluating immunological mechanisms in model animals and/or efficacy of one or more antivenoms against the venom of the analysed snake species; (iii) “biological activity”: assessing the enzymatic, toxic, and/or lethal (i.e., LD50) activity of the venom produced by the analysed snake species; (iv) “envenomation symptoms”: description of envenomation symptoms in humans resulting from snakebite accident; (v) “geographic venom variation”: comparing venom profiles, components, and/or biological activity between individuals belonging to the same snake species but coming from different populations and/or habitats across their natural range; (vi) “individual venom variation”: comparing venom profiles, components, and/or biological activity between individuals of the same snake species, with a focus on venom variation related to differences in age (i.e., ontogeny), sex, and/or diet; (vii) “interspecific venom variation”: comparing profiles, components, and/or activity of venoms produced by snakes belonging to different species; (viii) “prey specificity”: testing efficacy and/or efficiency of the venom of the analysed snake species against the preferred natural prey.

5.6. Chronological Trends

Information about the publication year of each analysed article was gathered in order to define the total number of publications per year, and thus identify the most and least productive years in terms of published articles. Using this information, we built cumulative curves in order to identify trends in terms of studied families, subfamilies, and research topics varied across the retrieved articles along the considered timeframe. The data obtained this way allowed to assess patterns of chronological variation in the above-mentioned categories.

5.7. Statistical AnalysesWe performed chi-squared (χ2) tests to assess the significance of the differences in terms of article coverage detected between snake taxa (i.e., family, subfamily, genus), countries, biogeographic realms, and topic categories. To investigate the presence of significant relationships between number of publications on snake venom and years from 1964 to 2021, we tested the following regression models: (i) 1st order polynomial, (ii) 2nd order polynomial, and (iii) 3rd order polynomial. We ranked the models on the basis of the corrected Akaike’s Information Criterion (AICc) [62], ultimately applying the model with the lowest AICc score considered as the best-fitting one. We applied the same method to also choose the best model to test the presence of significant relationships between the number of years that passed from 1964 to 2021 and the number of yearly papers covering each of the eight topic categories defined.

To investigate whether family, hazard category, and biogeographic realm of origin of the snake species retrieved from the analysed articles could be correlated with the number of articles dedicated to each one of them, we used Generalised Linear Models (GLM) assuming a Poisson distribution for the response variable. Country of origin and subfamily were excluded from the used predictors because the retrieved information relative to them was often fragmentary and ambiguous, and because they were nested in the predictors “biogeographic realm” and “family”, respectively. Collinearity between the three predictors considered (i.e., family, hazard category, biogeographic realm) was low (Variance Inflation Factors (VIF) always <5.11), thus we included all of them in the regression models generated. We built the models using the number of articles dedicated to each species as response variable, and all possible combinations of the three predictors considered. The produced models were ranked on the basis of their AICc score, considering the model with the lowest AICc score as the best-fitting one.

Polynomial regression models were generated using the software SPSS (version 13.0. [63]). All other analyses were performed in R environment (version 4.1.1 [64]). We used the packages vegan [65] and MuMIn to build the full set of Generalised Linear Models [66].