We extracted accepted names of all species of the three families according to the World Checklist of Vascular Plants (WCVP) V7 (Govaerts et al., 2021), totalling 21,111 species – 6,495, 496 and 14,120 from Apocynaceae, Loganiaceae and Rubiaceae respectively. We use Genus and Family names as categorical traits.

Due to the documented link between traditional medicinal usage and bioactivity, evidenced in, for example, (Krettli, 2009), we collected binary traits documenting the presence and absence of known antimalarial usage (Antimalarial Use) and general medicinal usage (Medicinal). To compile these data we conducted a comprehensive literature review of medicinal usage in the three plant families, along with data provided by the Medicinal Plant Names Services (MPNS, 2022) and references to medicinal usage on the Plants of the World Online (POWO, 2022).

As an extension of the ethnobotanical data, we included binary traits to capture whether a plant is commonly known – which we approximated by recording the presence of a Wikipedia ¹ page (Wiki Page) and the existence of a common name (Common Name). The existence of Wikipedia pages for species is determined by searching all species, subspecies and varieties (and their synonyms). Common name data are compiled from a variety of sources, outlined in the Supplementary Material , with the majority of the data coming from MPNS and the United States Department of Agriculture Plants Database (USDA, 2022b).

There is much evidence of the pharmacological and pharmaceutical importance of plant-derived alkaloids (Cordell et al., 2001; Dey et al., 2020; Howes et al., 2020; Daley and Cordell, 2021) and we have therefore collected binary traits on their presence/absence. These data were collected through a comprehensive literature review as well as metabolite data compiled from KNApSAcK (Afendi et al., 2012).

Though the coverage of the alkaloid data is relatively good (980 species with reported presence or absence), for the vast majority (97%) of these species, reports indicate a presence of alkaloids compared to 3% where alkaloids are absent. This may be the result of reporting bias, where publications are focused on species found to contain alkaloids and absences of alkaloids are not published. To assess the prevalence of the reporting bias, we contacted 11 authors of papers after the year 2000 that solely reported presences to ask if they had found any absences which they did not publish. We received responses from three authors detailing three species where alkaloids had been tested for and not found. Rather than being an issue of reporting bias, it may be the case that the vast majority of species in these families produce alkaloids. There is some evidence for this from studies testing large numbers of species for alkaloids where both presences and absences are reported e.g. (Soto-Sobenis et al., 2001). In either case, current data on the presence of alkaloids are relatively uninformative and so is not included in the following analysis. Instead, we use the collected data on alkaloids to catalogue which species have been tested for alkaloids. We use these data to create a binary trait (Tested for Alkaloids) which we use to analyse the relationship between phytochemical knowledge and knowledge of antiplasmodial activity.

An important plant trait indicating potent bioactivity is the degree to which a plant is toxic. As we aim to capture bioactivity in a broad sense, we compiled data on toxicity to any vertebrate and invertebrate animals. We have included this as a binary trait (Poisonous). Poison data were compiled from numerous sources detailing plants considered to be poisonous, outlined in the Supplementary Material , with the majority of the data coming from the LitTox resource (Royal Botanic Gardens, Kew, 2021).

As a major putative role of certain phytochemicals is to protect plants from herbivores (Maldonado et al., 2017), it is plausible that other defence mechanisms have a relation to bioactivity. Furthermore, certain biologically active compounds (e.g. some diterpene alkaloids) are biosynthesised in particular morphological structures (e.g. plant trichomes/hairs) of certain plants (Tomlinson et al., 2022). Here we assess the presence of emergences (hairs or spines) which we include as a binary trait (Emergence). Emergence data have been collated by Gentianales specialists, supplemented by the TRY plant trait database (Kattge et al., 2020) and POWO.

Another morphological trait we consider is plant life-form, which may correspond to occurrences of specific phytochemicals (de Almeida et al., 2005). To facilitate collection and coverage of morphological data, and as life-forms and presence of emergences are often well conserved within genera in these families, we include these traits by using the predominant state at the genus level. As multiple life-forms may appear within a single genus, the life-form data are one-hot encoded giving a set of binary traits (herb, liana, succulent, shrub, subshrub, tree). Life form data were initially retrieved from the WCVP and Flora do Brasil (Jardim Botânico do Rio de Janeiro, 2022), then reviewed and modified by Gentianales specialists.

To examine the relationship between prevalence of malaria in a given geographic area and the number of tested species, we collected data indicating which species are found in regions where malaria transmission occurs. We identified those regions from various sources, including the World Health Organization Database (WHO, 2022a) and the World Bank Development Indicators (The World Bank, 2022) (see Supplementary Material for full details). Regions indicated in these sources were then mapped onto the World Geographical Scheme for Recording Plant Distributions (Level 3) (Brummitt et al., 2001). We then used the WCVP distribution data to identify which species occur in these malarial regions (either native or introduced), assigned to each species as a binary trait In Malarial Region.

There is some evidence of environmental impacts on bioactive metabolite concentrations and diversity, for example, (Defossez et al., 2021). To characterise the environmental niche of species, we followed the methodology of Zu et al. (2021). We first extracted geographic occurrence records from the Global Biodiversity Information Facility (GBIF) ² for each species using the rgbif package (Chamberlain et al., 2022) in R. Occurrence data from GBIF contain many inconsistencies (Meyer et al., 2016). Initially we cleaned the data by removing: records collected before 1945, records with no given coordinates or impossible coordinates, records with coordinate uncertainty over 20km, records with rounded coordinates and records where the quantity of species occurrences (individual counts) is zero. Next, using the CoordinateCleaner package in R (Zizka et al., 2019) we removed: records with zero longitude or latitude, records with equal longitude and latitude, records outside reported country, records within country or province centroids, records in country capitals, records with institutional coordinates and records with GBIF Head Quarters coordinates. Finally, we discarded occurrences where species were reported to be outside of their native or introduced botanical regions according to the WCVP.

We quantified species’ environmental conditions using a set of 17 soil, climate, and topographic variables essential to plant survival, growth and reproduction. We extracted five soil traits (nitrogen content, pH, organic carbon stock (ocs), water capacity) from the SoilGrids database (Hengl et al., 2017; Poggio et al., 2021), which were averaged over a 30cm depth, as well as soil depth to bedrock. The eight bioclimatic traits we used were (bio1, bio4, bio10, bio11, bio12, bio15, bio16, bio17); representing temperature (mean annual, seasonality, daily mean of the warmest quarter, daily mean of the coldest quarter), precipitation (annual amount, seasonality, mean monthly amount of the wettest quarter, mean monthly amount of the driest quarter). These were extracted from the CHELSA database V2.1 (Karger et al., 2017; Karger et al., 2021). We also extracted the Köppen-Geiger climate classification (kg mode) from GloH2O (Beck et al., 2018). Elevation and breakline elevation were extracted from GMTED2010 (Danielson and Gesch, 2011) and slope was calculated from the elevation data using the terra package in R (Hijmans, 2022).

To match the resolution of the occurrence records, all environmental rasters were upscaled to 10 arc-minutes (c. 20 km) using the aggregate function of the terra package and environmental traits were extracted for each species occurrence using the extract function. For the continuous traits, median values were then calculated across all occurrences of each species and for the categorical variable kg mode the mode of all occurrences of each species was used. To capture coarse spatial information we also included median latitude and longitude for each species, calculated from the occurrence records.

To generate a comprehensive dataset of antiplasmodial activity, we conducted a thorough literature review for details of antiplasmodial tests in Apocynaceae, Loganiaceae and Rubiaceae and assigned activity labels to species based on the available reports of in vitro and in vivo studies. As with many biological datasets (Bender and Cortes-Ciriano, 2021), providing class labels is a nontrivial problem as there are many variations on the experiments and methods used for reporting activity. A detailed summary of the designated classification scheme we chose is given in the Supplementary Material . In general, for in vitro studies testing activity against Plasmodium parasites, the potency of IC50 values for crude extracts follows the definitions given in (Rasoanaivo et al., 2004a) i.e. < 10μg/ml is active and ≥ 10μg/ml is inactive. For tests of isolated compounds, according to the Medicines for Malaria Venture ³ compounds with IC50 values < 1 µM are designated as active and of interest for further investigation, thus we use this threshold in our data. For fractions, we use a threshold of 5 μg/ml, which in general corresponds with published author decisions of activity categories. For in vivo studies, we use the published author decisions regarding activity.

For some traits and datasets, presences are commonly reported but absences are not. For example, there are various datasets listing poisonous plants but published data on ‘safe’ plants are sparse. In many cases, this is likely a result of reporting bias, however there are multiple possible reasons for this. For certain traits there are presence biases e.g. in the case of poisons, once a plant has been found to be poisonous it can be reported as such; however if a plant is assessed for its toxicity, there are various caveats which limit the ability to confidently say the plant is safe. Examples of such caveats include the effect of extraction or preparation method on toxicity, the specific plant part tested, and which organisms the plant is toxic to. These variables exist in addition to methodological differences in assessing toxicity and also that in vitro studies may not correlate with effects in vivo (Houghton et al., 2007).

Where missing data give a strong indication of a genuine absence, i.e. for Common Name, Poisonous, Medicinal, Wiki Page, Antimalarial Use, Emergence, we take these pseudoabsences to be absences and fill missing values with 0. Missing values for other traits are left as NA and, where necessary, will be imputed.

In the preprocessing step, the categorical traits Genus, Family and kg mode are target encoded. The traits are then scaled by removing the mean and scaling to unit variance. We then use the scikit-learn (Pedregosa et al., 2011) k-Nearest Neighbor imputer, trained using the training data and the unlabelled data, to impute any missing values. Finally, we use Principal Component Analysis (PCA), implemented in scikit-learn, to reduce the dimensionality of the highly colinear continuous environmental traits. The PCA is trained using the training data and unlabelled data and the number of components used in the PCA is selected such that at least 80% of the variance is explained by the components. The traits In Malarial Region and Tested for Alkaloids were collected for the analysis of sampling bias rather than as predictive traits and so are not included in the machine learning models.

The Logit, SVC and XGB classifiers are trained as follows. Given a set of training folds and a test fold, hyperparameters of the models are tuned via cross validation on the training data using GridSearchCV (Pedregosa et al., 2011). In this step, F _0.5 is used as the evaluation metric and we tune a basic list of hyperparameters in order to minimise under/overfitting and to maximise F _0.5. For the Logit and SVC classifiers, we tune the regularisation parameter C, as well as the class_weight parameter. For the XGB classifier, we tune the max_depth parameter. Once the best hyperparameters for the models have been generated the models are retrained on all the given training data (with/without sample weights depending on the evaluation setting).

For the BNN classifier, we use two layers of 10 and 5 nodes, respectively, and tanh activation function. We train the model through 100,000 Markov chain Monte Carlo iterations, as implemented in npBNN (Silvestro and Andermann, 2020), with/without sample weights depending on the evaluation setting. We use 1,000 posterior samples of the parameters when generating predictions.

Following the scheme for classifying activity described in Section 2.1.7, we designated 132 species as active and 150 species as inactive, providing 282 labelled species from the 21,111 species in Apocynaceae, Loganiaceae and Rubiaceae. In these labelled data, all species are given trait values for each of the traits except for those trait values which rely on GBIF occurrence records where data are missing for five species.

Figure 2 provides a brief overview of the collected data, summarising relationships between some of the traits from all the collected data. The heatmap gives a visualisation of the co-occurrences of the binary traits, and the given values correspond to the mean values of traits in the y axis when traits in the x axis are present, while ‘All Species’ provides a comparison with mean values of traits in the underlying population. For example, 1% of all species are used traditionally for malaria while 18% of poisonous species are used traditionally for malaria. Similarly 10% of all species are used as traditional medicines while 77% of poisonous species are used as medicines. With regards to activity, the first column provides mean trait values for active species which indicate stark differences with the underlying population (e.g. 52% of active species are poisonous compared to 3% in the underlying population). However, these differences are more a reflection of the sampling biases rather than any strong relationships between the traits and antiplasmodial activity.

The most common feature motivating the selection of plants to test for antiplasmodial activity is traditional knowledge of use for malaria, for example (Andrade-Neto et al., 2003; Bourdy et al., 2004; Bertania et al., 2005; Ramalhete et al., 2008; Ezike et al., 2016; Taek et al., 2021). We found that 48% of labelled species are traditionally used for malaria while only 1% of species in the underlying population are traditionally used for malaria. Similarly, plants are frequently tested based on more general traditional medicinal usage (not specific to malaria), e.g. (Kaushik et al., 2013; Mothana et al., 2014; Singh et al., 2015; Satish et al., 2017). 77% of labelled species are traditionally used as medicines while 10% of species in the underlying population are traditionally used as medicines.

As previous successes in finding plants with antiplasmodial activity have linked their alkaloid content to the antiplasmodial activity, tests of antiplasmodial activity are often conducted on plants known/expected to contain alkaloids. For example (Wright et al., 1992; Solis et al., 1995; Likhitwitayawuid et al., 1999; Weniger et al., 2001; Mitaine-Offer et al., 2002; Federici et al., 2009). Moreover, in many reports where plants are tested for antiplasmodial activity, those studies also include tests for (and find) alkaloids e.g. (Likhitwitayawuid et al., 1999; Muhammad et al., 2003; Suksamrarn et al., 2003; Wong et al., 2011). As a result, 69% of labelled species and 82% of active species have been tested for presence of alkaloids, while only 5% of species in the underlying population have been tested for presence of alkaloids.

Another potential factor influencing sampling is the geographic location of species, i.e. plants occurring in regions with malaria are commonly selected to test for antiplasmodial activity, for example (Rasoanaivo et al., 2004b; Bertania et al., 2005; Al-Musayeib et al., 2012; Kantamreddi and Wright, 2012; Taek et al., 2021). As a result, 99% of labelled species are found in malarial regions compared to 89% in the underlying population. In fact, there is only one tested species which is not found in a malarial region (Gardenia urvillei Montrouz. (Rubiaceae) which is native to New Caledonia) and three Ochrosia Juss. (Apocynaceae) species (native to Fiji, Tonga and New Caledonia) whose activity is known through the presence of antiplasmodial compounds (not themselves explicitly tested) which are not found in malarial regions.

It is also common to test plants taxonomically related to known antiplasmodial plants (Weenen et al., 1990; Frédérich et al., 2002; Philippe et al., 2005; dos Santos Torres et al., 2013; Brandão et al., 2020). For example, some genera known to contain active species are frequently tested e.g. Aspidosperma (Apocynaceae: Gentianales) (18 labelled species) and Strychnos (Loganiaceae: Gentianales) (36 labelled species).

For almost all the quantitative traits, the difference between the labelled data and underlying population (as measured by Chi-squared test for the discrete traits and the Kolmogorov–Smirnov 2-Sample test for the continuous traits) is significant (corrected p values < 0.05) with the exception of life-forms (lianas and succulents). The most diverging traits are Antimalarial Use and Tested for Alkaloids (corrected p values = 0, Chi-squared statistic 3137 and 2218 respectively). We can therefore conclude that the labelled data significantly differ from the underlying population. Overall, it is apparent that the approaches used to select plants for antiplasmodial tests have biased the available data on antiplasmodial activity.

When testing the Correction Model in 10 iterations of 10-fold stratified cross validation, the mean Brier score was 0.0097 (SD = 0.001), indicating an accurate fit to the data and so, a reliable prediction of the selection probability. A visual comparison of the bias-corrected data and the underlying population is given in Figures 3 , 4 . For readability, the mean values of the continuous traits are rescaled between 0 and 1 using the MinMaxScaler from scikit-learn (Pedregosa et al., 2011). We can see that for the majority of the traits, the mean values of the corrected data closely resemble the underlying population compared to the values in the labelled data.

For those species that we deem highly unlikely to be tested, ( P ( s | x ) < 0.0014 ) only 2 such species are in the labelled data where one is known to be active, while 9997 are in the unlabelled data. When the SVC model is trained on all the available data and used to predict the activity of these species in the unlabelled data, 2358 are estimated to be active. This gives a 95% confidence interval of 1300 – 1522 active species when the model precision is accounted for. Note that this is a conservative approximation as we are only considering species that the model predicts to be active and correcting for the estimated false positives. However, as visible in the Precision-Recall curves, recall of the models is not perfect and it is highly likely that there are also a significant number of species that the model predicts to be inactive species which are in fact active.