Methodology for the identification and analysis of species associations in biological communities using statistical distributions

Methodology for the identification and analysis of species associations in biological communities using statistical distributions | 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 Methodology for the identification and analysis of species associations in biological communities using statistical distributions Dmitriy G. Seleznev , Alexander A. Prokin , Ekaterina M. Kurina doi: https://doi.org/10.1101/2024.08.02.606312 Dmitriy G. Seleznev 1 Papanin Institute for Biology of Inland Waters, Russian Academy of Sciences , Borok, Yaroslavl Oblsat, 152742, Russia 2 A. N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences , Leninsky Prospekt 33, Moscow, 119071, Russia Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: dmitriy{at}seleznev.name Alexander A. Prokin 1 Papanin Institute for Biology of Inland Waters, Russian Academy of Sciences , Borok, Yaroslavl Oblsat, 152742, Russia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ekaterina M. Kurina 2 A. N. Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences , Leninsky Prospekt 33, Moscow, 119071, Russia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Preview PDF Abstract A methodology for identifying pairs of mutually associated and mutually incompatible species pairs based on qualitative ecological data in the ‘species-samples’ format is described. Statistical criteria for assessing the direction and strength of species associations using discrete hypergeometric and binomial distributions are provided. Additionally, there are proposed methods for visualising selected associations using undirected graphs, combining associated species into groups, and analysing them. An open access web application implementing the described methodology is also presented. The methodology is illustrated with data on 37 alien macrozoobenthic species from 255 samples collected in the Volga River reservoirs. Introduction In nature, the mutual association or incompatibility of species depends on their interspecific relationships, which are the basis for the origination and existence of biological systems. Species can form assemblages, including communities, based on interspecific interactions, which may be direct (competition, predation, parasitism, etc.) or indirect, i.e., mediated or transmitted by a third species or environment (habit facilitation, keystone predation, exploitation, competition, etc.) ( Odum & Barrett, 1971 ; Wootton, 1994 ; Moon et al., 2010 ; Sotomayor & Lortie, 2015 ; etc.). To study such interactions using species presence-absence matrix data, it is necessary to define pairs of species that are mutually associated or incompatible. At present, ecological studies rarely address this problem. This may be due to the lack of obvious mathematical methods for this purpose. For this task, in addition to the correlation coefficients and similarity/dissimilarity measures traditionally used in ecology, such as the Jaccard and Ochiai indices ( Jaccard, 1912 ; Ochiai, 1957 ) or variations of the Sörensen measure (Sörensen, 1948; Bray & Curtis, 1957 ), some specialised matrix indices have been proposed. For instance, there are the checkerboard index CHECKER ( Diamond, 1975 ), the C-score index ( Stone & Roberts, 1990 ), Schluter’s variance test ( Schluter, 1984 ), and the nestedness temperature ( Atmar & Patterson, 1993 ). In addition to the necessary normalisation, all of these methods are based on the arithmetic calculation of species co-occurrence in the qualitative matrix. Therefore, the researcher does not have comprehensible criteria for identifying mutually associated species, even such unreliable one as the significance level. In the case of an observed co-occurrence of two species, the significance of the index, as well as the decision on the association of species, directly depends on the total number of samples, which matrix indices often do not take into account. Indeed, having obtained a normalised C-score index of 0.96, the researcher has no formal reason to assert that these species are mutually incompatible and that their low co-occurrence is not due to chance. This issue has been proposed to be resolved by creating a randomly generated null model or a theoretical distribution of the matrix index with the observed species occurrence and the total number of samples ( Gotelli, 2000 ). It is crucial to carefully choose a randomisation algorithm because different permutation schemes can yield vastly different results. The decision about the association or incompatibility of a pair of species is made by comparing the observed index and given confidence intervals of the theoretical distribution. This approach is proposed to be used by the most common ecological R package “vegan” ( Oksanen et al., 2019 ), as well as by the specialised package “EcoSimR” ( Gotelli et al., 2015 ). In practice, the study of species associations using matrix indices and a simulated null model can encounter certain difficulties. An analysis of 255 samples of alien macrozoobenthic species from Volga River reservoirs with the C-score index and a permutational fixed-equiprobable (FE) null model ( Kurina & Seleznev, 2019 ) revealed 82 positively associated and 5 negatively associated pairs of species. However, 27 positively associated pairs had a normalised C-score index greater than 0.4, reaching a value of 0.55 for the pair of species Katamysis warpachowskyi (Sars, 1893) (n 1 =16) and Pterocuma rostrata (Sars, 1894) (n 2 =25), which co-occurred in 5 samples. Therefore, interpreting the normalised C-score index is a challenging task. Furthermore, the tendency of permutational null models to overdiagnose, producing the type I error, has been noted ( Sanderson, 2004 ; Gotelli et al., 2010 ), as illustrated by the example above. An estimation of the significance level of indices using a randomisation procedure has a probabilistic nature. This approach may lead to problems with the reproducibility of the results obtained, as well as their instability. In a set of samples with several pairs of species having the same individual and joint co-occurrence, only some of them may show a significant association. Moreover, the set of these pairs may change during repeated modelling. Identification of associated pairs of species using statistical distributions A probabilistic model has recently been proposed for the precise calculation of the probability of a certain event: the co-occurrence of two species in x out of N samples ( Veech, 2013 ). In that model, the classical definition of probability as a quotient of the number of successful outcomes and the total number of outcomes was used. This formula was subsequently reduced to the probability mass function (PMF) of the hypergeometric distribution ( Griffith et al., 2016 ), which describes the probability of obtaining a specific number of successes in draws without replacement from a finite population. When applied to the co-occurrence of two species, the hypergeometric distribution allows calculations of the probability of finding particular species in samples that already contain other species with a total number of samples N and given that species are distributed independently within the community. Let m denote the occurrence of a rarer species and n denote the occurrence of a more common species. The probability that the species will be observed together precisely in x samples of a total of N samples is determined by the probability mass function of the hypergeometric distribution: Therefore, the task of determining association or, conversely, incompatibility of two species at a given significance level can be reduced to checking whether the observed joint co-occurrence value falls out of the confidence intervals of the hypergeometric cumulative distribution function (CDF): In other words, the hypergeometric CDF allows one to test the null hypothesis that species are independent of each other and that their co-occurrence in samples is random. An alternative is the assumption that the observed positive or negative association of species is not accidental and is due to their biological characteristics or ecological reasons. It should be noted that the Fisher’s exact test on the contingency table 2×2 (species A, species B; occurrence, co-occurrence) leads to the same result. This is because, as Ronald Fisher points out, such a test for frequency independence leads to a hypergeometric distribution. A discrete hypergeometric distribution tends to be a binomial distribution with an infinite increase in the general population size ( Fig. 1 ). In this case, a “draws without replacement” model becomes “draws with replacement,” also known as the Bernoulli scheme. When using the hypergeometric distribution, we considered all samples taken as the general population, while samples in which one of the species was found were considered as the draws. Download figure Open in new tab Figure 1. Simulated probability mass function of the hypergeometric (black curve) and binomial (red curve) distributions with 10 6 simulations of 50 trials and a success probability of 1/3. The vertical lines indicate the boundaries of the 90% confidence intervals. Download figure Open in new tab Figure 2. The graph of positively associated pairs of macrozoobenthic species in rivers of the East European Plain. A binomial distribution with a significance level of 0.05 and a multi-level modularity optimisation algorithm for clustering were employed. A – Kamada–Kawai algorithm node layout, B – Multidimensional scaling (MDS) node layout. For species abbreviations, see Golovatyuk & Seleznev (in press). Let us reformulate the aforementioned conditions. We will consider the entire potentially possible (infinite) number of samples in a given habitat as the general population. The taken samples will be regarded as draws, and the joint co-occurrence of two species in the samples will be defined as draw success. Since we cannot know the a priori probability of success a priori , we will use it’s a posteriori estimation ( Bowers & Brown, 1982 ; Veech, 2006 ): The probability of the joint co-occurrence of two species in x samples under these conditions is defined by the probability mass function of the binomial distribution: The estimation of the association rate at the given level of significance is also performed by comparing the observed co-occurrence value with the confidence intervals of the binomial cumulative distribution function. The possibility of making decisions about a probabilistic process based solely on the significance level with a conventional threshold of 0.05 is currently a subject of active debate within the scientific community (see studies published in: The American Statistician, 2019 , Vol. 73, Issue Supp 1; Wasserstein & Lazar, 2016 ). Given that we have defined the distribution function of the probabilistic process a priori, it is possible to determine pairs of mutually associated or incompatible species using an approach based on Bayesian inference rather than the traditional method of testing statistical hypotheses. According to this approach, one or more uniformly distributed characteristics are taken as the evaluated parameter. These may be the occurrence of species or the total number of samples. The probability mass function is then used as a generative function. The posterior distribution, filtered by the observed values of the evaluated parameters, allows us to estimate the probability of the joint co-occurrence of two species. On this basis, we can make one of three decisions: both species occur together by chance, both species are associated, or these two species are incompatible ( Fig. 4 ). It should be noted that the application of Bayesian inference may lead to issues with the reproducibility and stability of the results due to the use of a random variable as a prior distribution, as is the case with the use of randomisation procedures. The analysis can also use a transposed matrix with initial qualitative data in the “samples-species” format. As a result of the analysis, we obtain a list of pairs of samples whose similarity or difference in species composition significantly differed from random. This result can be used to determine the faunal similarity of habitats as well as species habitat preferences. Visualisation and grouping of mutually associated species It is not uncommon for a single species to be associated with multiple others. It is convenient to visualise such cases using an undirected graph, in which the nodes represent the species from the list of associated ones and the edges represent the identified associations ( Fig. 3 ). The size of the node marker will be proportional to the occurrence or the total abundance (or biomass) of the species, and the edge thickness represents the degree of association. This degree can be defined by various methods, including the use of one of the previously described normalised matrix indices or the probability of obtaining the observed joint co-occurrence, which is calculated using the selected distribution probability mass function. If desired, pairs of incompatible species can be added to the graph, indicating positive or negative associations with different edge colors. Download figure Open in new tab Figure 3. Graph visualisation of positively associated species: A – FE permutation model of the C-score index; B – hypergeometric distribution. The significance level was 0.01. For species abbreviations see ( Kurina & Seleznev 2019 ). Download figure Open in new tab Figure 4. The posterior distribution of the joint co-occurrence of the two species. The blue lines show the limits of the double-sided 90% confidence intervals, and the red lines show the limits of the double-sided 98% confidence intervals. The layout of graph nodes can be defined using a variety of algorithms, among which we can highlight the Kamada–Kawai algorithm ( Kamada & Kawai, 1989 ), which is based on the physical model of springs ( Fig. 2A ), and the force-directed Fruchterman–Reingold algorithm ( Fruchterman & Reingold, 1991 ). A further approach is possible, whereby additional meaning is added to the layout of the graph nodes. In particular, one may employ one of the ordination methods based on the frequency of occurrence, total abundance, or biomass of the associated species. In this case, the layout of the nodes will be defined by species scores resulting from the ordination procedure ( Fig. 2B ). It is convenient to combine groups of mutually associated species into clusters, which in graph theory terms are referred to as “communities”. The majority of graph node clustering algorithms are based on modularity maximisation, which estimates the extent to which the density of intragroup edges in a given partition is greater than that in a random partition of the same number of edges. There are numerous algorithms of this kind, and they do not require the number of clusters to be prespecified, unlike k-means clustering. In each cluster, the species with the largest number of the strongest connections can be identified as the “core” of the species assemblage, which determines its properties. The species clustering results can be used as input data for the analysis of species assemblages and interactions between species within them. Analysis of mutually associated species groups It is obvious that the selected groups will differ in both the density of connections within them and in the average connection strength. The density of connections (or cluster connectivity) can be calculated as the quotient of the number of realised connections and the number of all possible connections in a cluster: Cn = 2k/(n·(n−1)) , where k is the number of realised connections and n is the number of species in the cluster. The average connection strength in the cluster can be estimated as the average of any edge thickness values proposed above or as the average normalised joint co-occurrence Cf = 2·x/(n+m) , where n and m are the number of samples containing the first and second species, respectively, and x is the number of samples containing both species together. Let us reduce the value of the cumulative distribution function P x to fractions of one. For negatively associated pairs, the value is given by 1− P x /p-value . For positively associated pairs, the value is given by 1−(1−P x )/p-value . The total connectivity strength, which characterises both the density and connection strength in the cluster, is calculated as the average value of the distribution function in fractions of one of all significant connections, including zeros for insignificant connections. Consider a cluster of three species in which two pairs of positively associated species have been identified. The significance level for these associations was set at 0.05. The inverse values of the cumulative distribution function of these associations (1−P x ) were 0.04 and 0.01, respectively. The total connectivity strength of this cluster is therefore (0.2 + 0.8 + 0)/3 = 0.33. This characteristic is normalised in the range [0,1], thus enabling comparisons between different clusters. Furthermore, it is balanced, in that if there are few statistically significant associations in the cluster, a large number of zeros will affect the total connectivity strength, and conversely, if there are many significant associations, the connections’ strength will affect it more. Practical application The proposed statistical approach was employed for the repeated analysis of data on the occurrence of 37 alien macrozoobenthic species in 255 samples from Volga River reservoirs. Previously, these data were analysed using the C-score index, 1000 FE permutations and a double-sided 98% confidence interval ( Kurina & Seleznev, 2019 ). The results of the comparison of the models are presented in Table 1 . View this table: View inline View popup Download powerpoint Table 1. Comparison of models of associated species: A – C-score index model, B – hypergeometric distribution model with confidence intervals, C – hypergeometric distribution model with Bayesian credible intervals, D – binomial distribution model with confidence intervals, E – binomial distribution model with Bayesian credible intervals. A double-sided 98% confidence interval was used; * – averaged values for 5 simulations. Models based on the C-score index and hypergeometric distribution yielded comparable results, although there was a slight discrepancy in the set of species pairs. However, there was a significant difference in the computational resources required for each method. The confidence interval method produces a slightly more conservative result than Bayesian inference in both statistical distribution models. The result of the binomial distribution model is more conservative due to more “heavy tails” in the region of low and high values of x ( Fig. 1 ). Based on the results of the analysis using the permutation model of the C-score index and the confidence intervals of the hypergeometric distribution, we constructed graphs ( Fig. 3 ) and compared them. The thickness of the edges was calculated using the value of the cumulative distribution function P x reduced to fractions of one, as described above. The node marker size was proportional to the species occurrence frequency. The node layout is generated using the Kamada–Kawai algorithm, while clustering is performed by the multi-level modularity optimisation algorithm ( Blondel et al., 2008 ). A comparison of the graphs of associated species obtained using the C-score index ( Fig. 3A ) and the hypergeometric distribution ( Fig. 3B ) revealed that the latter had greater modularity (0.41 vs. 0.32). This is due to the isolation of some species into separate clusters, while the number of intergroup connections is approximately equal (16 vs. 18). At the same time, two main groups of species remain: deep-water (I) and shore (II) species assemblages. Additionally, three species are separated into an intermediate, depth-independent assemblage associated with sandy sediments. Thus, the new grouping of associated species allows us to identify an additional abiotic factor that affects the spatial distribution of species. The two groups obtained by the two different methods are analysed below ( Table 2 ). View this table: View inline View popup Table 2. Characteristics of deep-water (I) and shore (II) clusters of associated macrozoobenthic species selected using the C-score index (C-score) and hypergeometric distribution function (HyperH). The characteristics of the clusters demonstrated that when a hypergeometric distribution was used, there were three pairs of species fewer in both groups, which resulted in an increase in connectivity. Despite a slight decrease in the average connection strength, the total connectivity strength increased. The aforementioned pair of species, Katamysis warpachowskyi (observed in 16 samples) and Pterocuma rostrata (observed in 25 samples) occurred together in 5 out of 255 samples. The pair was analysed using different methods. The C-score index for this pair is 220 and does not always fall outside the upper 99% confidence interval of the simulated index distribution. Consequently, the mutual association of these species is on the border of statistical significance. The value of the hypergeometric distribution function, i.e., the probability of obtaining the actual (5) or greater co-occurrence of species for the given parameters (individual occurrence), is 0.0125. This value does not fall outside the upper 99% limit of the confidence interval. In order to solve the same task, we use Bayesian inference. As a prior distribution, we simulate 10 6 values for the total number of samples N , which are uniformly distributed over the range [25, 255]. We will use the probability mass function of the hypergeometric distribution as a generative function, with the observed individual occurrence of species (25, 16) and the simulated total number of samples. Subsequently, the posterior distribution, which was filtered by the actual number of samples (255), will yielded the distribution of the number of joint co-occurrences of the two species ( Fig. 4 ). Fig. 4 demonstrates that the observed value (5) does not exceed the upper 99% limit of the confidence interval. Consequently, the application of Bayesian inference with a 2% probability of type I error does not allow us to conclude that the joint co-occurrence of Katamysis warpachowskyi and Pterocuma rostrata is statistically different from random. The publication of Veech (2013) represents the results of applying that author’s probabilistic model to 10 different datasets with a significance level of 0.05. Comparisons were made between the published results and those obtained by the methods proposed above on data taken from Galapagos Island finches ( Sanderson, 2000 ), Great Basin Desert rodents ( Bowers & Brown, 1982 ), carabid beetles on islands in lakes of Poland ( Ulrich & Zalewski, 2006 ; Gotelli & Ulrich, 2010 ), as well as bog and forest ants in New England ( Gotelli & Ellison, 2002 ) ( Table 3 ). View this table: View inline View popup Download powerpoint Table 3. Comparison of different models of associated species pairs. The number of positive and negative associations is given with a slash. * – rounded results of 5 simulation series; ** – from ( Veech, 2013 ). The results of the probabilistic model and the hypergeometric distribution model are found to be highly correlated, with the FE resampling null model exhibiting a tendency to overestimate the number of positively associated pairs and underestimate the number of negatively associated pairs. The only exceptions to this are the bog and forest ant datasets, for which both the probabilistic model and the FE resampling null model were found to be more conservative than the hypergeometric distribution model. It is also worth noting that the binomial distribution model is the most conservative of all those considered. Ecological applications In order to ensure the successful using of the proposed methodology in future investigations, it is essential to have a clear understanding of the ecological meaning of species falling into an associated pair or group. To date, the methodology has been employed in six studies, three of which were based on macrozoobenthos data (Kurina et al., 2019; Prokin et al., 2021 ; Kurina et al., 2022 ) and three on zooplankton data ( Kotov et al., 2022 ; Chertoprud et al., 2023; Andreeva et al., 2023 ). To identify associated pairs of species, the initial data was filtered in various ways in all studies except ( Prokin et al., 2021 ), thus allowing us to discuss analyses of taxocoenoses rather than complete biological communities. In all our studies, the principal criteria for categorising species into associated groups were environmental factors, such as depth, current velocity, salinity and other physico-chemical parameters of the water, the presence of water plants, the type of bottom sediment or the type and level of anthropogenic influence. It is likely that the common requirements of environmental factors also explain the occurrence of ecosystem engineers such as the quagga mussel Dreissena bugensis (Andrusov, 1897) and the zebra mussel D. polymorpha (Pallas, 1771) with connected species in the associated groups (Kurina et al., 2019; Kurina et al., 2022 ). A preliminary analysis of the results obtained thus far suggests that the observed associations between species indicate a similar selectivity to a set of environmental factors that limit their distribution at the scale of a particular study. However, changes in the area or scale of the study may result in considerable alterations to the hierarchy and importance of factors and consequently, the groups of associated species. To test and refine this assumption, further studies on data from different groups of organisms, as well as meta-analyses of these studies, are required. Sofware The senior author’s open access web application ( http://apps.ibiw.ru/coobs ) implements the described methodology and allows users to obtain results in text and graphic formats. As an analysis method, hypergeometric or binomial distributions can be selected. The choice of pairs of associated species is performed using confidence intervals or Bayesian credible intervals of selected distribution at the given significance level. In the case of Bayesian credible intervals, 10 6 uniformly distributed total number of samples is used as a prior distribution, with the individual species occurrences fixed. An undirected graph is constructed by several methods, using positive or negative associations only or using all selected associations. To group graph nodes, various algorithms are used that take into account the strength of the connections between graph nodes as a positive edge weight factor. In addition, it is possible to highlight the “Spin glass” algorithm ( Traag & Bruggeman, 2008 ), which is capable of handling negative associations. This algorithm is able to work correctly with negative edge weights, in contrast to other methods that consider them to be null. The “Optimal community structure” algorithm ( Brandes et al., 2008 ) is also worthy of note because it solves the graph node grouping task as an integer programming problem and builds a graph with maximum modularity. However, this method also fails to consider negative associations, such as species incompatibility, and has exponential complexity. This is because the task of maximising modularity appears to belong to the NP class ( Brandes et al., 2006 ). The application was developed in the R 4.2 statistical analysis environment using the package igraph ( Csardi & Nepusz, 2006 ). The web interface was designed using the Shiny package ( Chang et al., 2021 ). Conclusion The proposed methodology differs from those that use permutation null models in its simplicity of implementation and in its low requirements for computing resources. This approach allows to obtain reproducible and easily interpreted results with a given significance level in the form of a list of pairs of mutually associated or incompatible species. This list can later be visualised and analysed both mathematically and in terms of studying the ecological characteristics of individual species and the features of the formation of species assemblages and communities. One disadvantage of the methodology is that it uses only qualitative data on the presence or absence of species in the samples, without considering their quantitative characteristics, such as total abundance or biomass. It is possible that two species may occur in samples together rarely but in large numbers and separately often but in small numbers (or vice versa). However, this phenomenon is not detected when discrete distributions are used. To account for quantitative data, it seems that the use of continuous counterparts of discrete distributions or some other methods will be required. Funding The work was carried out under the framework of the Ministry of Science and Higher Education of the Russian Federation, project No. 124032500016-4 (DGS and AAP) and project No. 0089-2021-0006 (EMK). Acknowledgements The authors are grateful to D.Y. Medvedeva and P.N. Petrov (Moscow) for their help with improving the English of the manuscript and V.K. Shitikov (Tolyatti) for a constructive discussion of this study. References ↵ Andreeva L.V. , Seleznev D.G. , Garibian P.G. , Kotov A.A. Local Faunistic Associations of Water Fleas (Crustacea: Cladocera) in the Water Bodies of the Left Bank of the Lena River (Yakutia, Russia) . Inland Water Biology . 2023 . doi: 10.1134/S1995082923050024 OpenUrl CrossRef ↵ Atmar , W. , Patterson , B. D. ( 1993 ) The measure of order and disorder in the distribution of species in fragmented habitat . Oecologia . Vol. 96 . P. 373 – 382 OpenUrl CrossRef Web of Science ↵ Blondel , V. D. , Guillaume , J.-L. , Lambiotte , R. , Lefebvre , E. ( 2008 ) Fast unfolding of communities in large networks . Journal of Statistical Mechanics: Theory and Experiment . P10008 . doi: 10.1088/1742-5468/2008/10/P10008 OpenUrl CrossRef ↵ Bowers , M. A. , Brown , J. H. ( 1982 ) Body size and coexistence in desert rodents: chance or community structure? Ecology . Vol. 63 , No. 2 . P. 391 – 400 . doi: 10.2307/1938957 OpenUrl CrossRef Web of Science ↵ Brandes , U. , Delling , D. , Gaertler , M. , Görke , R. , Hoefer , M. , Nikolski , Z. & Wagner , D. ( 2006 ) On Modularity – NP-Completeness and Beyond . Fakultät für Informatik, Universität Karlsruhe . Vol. 19 . 33 pp. doi: 10.5445/IR/1000005777 OpenUrl CrossRef ↵ Brandes , U. , Delling , D. , Gaertler , M. , Gorke R. , Hoefer , M. , Nikoloski , Z. & Wagner , D. ( 2008 ) On modularity clustering . IEEE. Transactions on knowledge and data engineering . Vol. 20 , Issue 2 . P. 172 – 188 . doi: 10.1109/TKDE.2007.190689 OpenUrl CrossRef Web of Science ↵ Bray , J. R. , Curtis , J.T. ( 1957 ) An ordination of upland forest communities of southern Wisconsin . Ecological Monographs . Vol. 27 , Issue 4 : P. 325 – 349 . doi: 10.2307/1942268 OpenUrl CrossRef ↵ Chang , W. , Cheng , J. , Allaire , J. , Sievert , C. , Schloerke , B. , Xie , Y. , Allen , J. , McPherson , J. , Dipert , A. & Borges , B. ( 2021 ). Shiny: web application framework for R . R package version 1.7.1 . https://CRAN.R-project.org/package=shiny Chertoptud E.S. , Seleznev D.G. , Garibian P.G. , Kotov A.A. Microcrustaceans (Cladocera and Copepoda) of the Boreal/Tropical Transition Zone in the Russian Far East: A Case Study of Species Associations in Three Large Lakes . Diversity . 2023 . doi: 10.3390/d15030338 OpenUrl CrossRef ↵ Csardi , G. , Nepusz , T. ( 2006 ) The igraph software package for complex network research . International journal of complex systems . Vol. 1695 ( 5 ). P. 1 – 9 . OpenUrl ↵ M.L. Cody and J.M. Diamond Diamond , J.M. ( 1975 ) Assembly of species communities . In: Ecology and Evolution of communities ( M.L. Cody and J.M. Diamond , eds.). Belknap Press , Cambridge. (MA): Harvard Univ. Press . P. 342 – 444 . ↵ Fruchterman , T.M.J. , Reingold , E.M. ( 1991 ). Graph drawing by force-directed placement . Software – Practice and Experience . Vol. 21 , Issue 11 . P. 1129 – 1164 . doi: 10.1002/spe.4380211102 OpenUrl CrossRef Web of Science Golovatyuk , L.V. , Seleznev , D.G. (in prep.) Analysis of river macrozoobenthos species associations on the territory of the East European Plain under the conditions of changing natural zones . ↵ Gotelli , N.J. ( 2000 ) Null model analysis of species co-occurrence patterns . Ecology . Vol. 81 , Issue 9 . P. 2606 – 2621 . doi: 10.1890/0012-9658 (2000)081[2606:NMAOSC]2.0.CO;2 OpenUrl CrossRef PubMed Web of Science ↵ Gotelli , N.J. , Ellison , A.M. ( 2002 ) Biogeography at a regional scale: determinants of ant species density in New England bogs and forests . Ecology . Vol. 83 , Issue 6 , P. 1604 – 1609 . doi: 10.1890/0012-9658 (2002)083[1604:BAARSD]2.0.CO;2 OpenUrl CrossRef Web of Science ↵ Gotelli , N.J. , Ulrich , W. ( 2010 ) The empirical Bayes approaches a tool to identify non-random species associations . Oecologia . Vol. 162 . P. 463 – 477 . doi: 10.1007/s00442-009-1474-y OpenUrl CrossRef PubMed ↵ Gotelli , N.J. , Hart , E.M. , Ellison , A.M. ( 2015 ) EcoSimR: null model analysis for ecological data . R package version 0.1.0 . http://github.com/gotellilab/EcoSimR ↵ Griffith , D.M. , Veech , J.A. , Marsh , C.J. ( 2016 ) Cooccur: probabilistic species cooccurrence analysis in R . Journal of Statistical Software . Vol. 69 , Code Snippet 2. P. 1 – 17 . doi: 10.18637/jss.v069.c02 OpenUrl CrossRef ↵ Jaccard , P. , ( 1912 ). The distribution of the flora in the Alpine zone . New Phytologist . Vol. 11 , Issue 2 . P. 37 – 50 . doi: 10.1111/j.1469-8137.1912.tb05611.x OpenUrl CrossRef ↵ Kamada , T. , Kawai , S. ( 1989 ) An algorithm for drawing general undirected graphs . Information Processing Letters . Vol. 31 , Issue 1 . P. 7 – 15 . doi: 10.1016/0020-0190(89)90102 OpenUrl CrossRef Web of Science ↵ Kotov A.A. , Seleznev D.G. , Garibian P.G. , Korovchnsky N.M. Neretina A.N. , Sinev A.Y. , Jeong H.-G. , Yang H.-M. Lee W. History of Colonization of Jeju Island (Republic of Korea) by the Water Fleas (Crustacea: Cladocera) Is Reflected by the Seasonal Changes in Their Fauna and Species Associations . Water . 2022 . doi: 10.3390/w14213394 OpenUrl CrossRef ↵ Kurina , E.M. , Seleznev , D.G. ( 2019 ) Analysis of the patterns of organization of species complexes of Ponto-Caspian and Ponto-Azovian macrozoobenthosin the Middle and Lower Volga Reservoirs . Russian Journal of Ecology . Vol. 50 . P. 65 – 74 . doi: 10.1134/S1067413619010053 OpenUrl CrossRef ↵ Kurina E.M. , Seleznev D.G. , Sherysheva N.G. Distribution of Alien Species of Macrozoobenthos and the Species Cenotic Complexes in the Kama Reservoirs . Russian Journal of Biological Invasions . 2022 . doi: 10.1134/S2075111722010088 OpenUrl CrossRef ↵ Moon , D. C. , Moon , J. & Keagy , A. ( 2010 ). Direct and Indirect Interactions . Nature Education Knowledge . Vol. 3 , No 10 . 50 . https://www.nature.com/scitable/knowledge/library/direct-and-indirect-interactions-15650000 OpenUrl ↵ Odum , E.P. , Barrett , G.W. ( 1971 ) Fundamentals of Ecology . Third edition. Philadelphia : Saunders . 574 pp. ↵ Ochiai , A. ( 1957 ) Zoogeographical studies on the soleoid fishes found in Japan and its neigbouring regions . Bulletin of the Japanese Society of Scientific Fisheries . Vol. 22 , No. 9 . P. 522 – 525 . doi: 10.2331/suisan.22.522 OpenUrl CrossRef ↵ Oksanen J. , Blanchet F. G. , Friendly M. , Kindt R. , Legendre P. , McGlinn D. , Minchin P. R. , O’Hara R. B. , Simpson G. L. , Solymos P. , Stevens M. H. H. , Szoecs E. and Wagner H. ( 2019 ). vegan: Community Ecology Package . R package version 2.5-6 . https://CRAN.R-project.org/package=vegan (accessed on 10 August 2023 ) ↵ Prokin A.A. , Seleznev D.G. , Tsvetkov A.I. Influence of environmental factors on the interannual variability of macrozoobenthos of the floodplain lakes . Ecosystem Transformation . 2021 . doi: 10.23859/estr-210306 OpenUrl CrossRef ↵ Sanderson , J.G. ( 2000 ) Testing ecological patterns: a well-known algorithm from computer science aids the evaluation of species distributions . American Scientist . Vol. 88 , No. 14 . P. 332 – 339 . doi: 10.1511/2000.29.332 OpenUrl CrossRef GeoRef ↵ Sanderson , J.G. ( 2004 ) Null model analysis of communities on gradients . Journal of Biogeography . Vol. 31 , No. 6 . P. 879 – 883 . doi: 10.1111/j.1365-2699.2004.01054.x OpenUrl CrossRef ↵ Schluter , D. ( 1984 ) A variance test for detecting species associations, with some example applications . Ecology . Vol. 65 , No. 3 . P. 998 – 1005 . doi: 10.2307/1938071 OpenUrl CrossRef Web of Science Sorensen , T. ( 1948 ) A Method of establishing groups of equal amplitude in plant sociology based on similarity of species content and its application to analyses of the vegetation on Danish Commons . Biologiske Skrifter. Bind V , No. 4 . Ikommissionhos E. Munksgaard, København . P. 1 – 34 . ↵ Sotomayor , D. A. , and C. J. Lortie . ( 2015 ). Indirect interactions in terrestrial plant communities: emerging patterns and research gaps . Ecosphere . Vo. 6 , No 6 . 103 . doi: 10.1890/ES14-00117.1 OpenUrl CrossRef ↵ Stone , L. , Roberts , A. ( 1990 ) The checkerboard score and species distributions . Oecologia . Vol. 85 . P. 74 – 79 . doi: 10.1007/BF00317345 OpenUrl CrossRef Web of Science ↵ The American Statistician ( 2019 ) Volume 73 , Issue Supp 1 : Statistical inference in the 21st century: A World beyond p < 0.05..401 pp . ↵ Traag , V.A. , Bruggeman , J. ( 2008 ) Community detection in networks with positive and negative links . Physical Review E . Vol. 80 , Issue 3 . P036115 . doi: 10.1103/PhysRevE.80.036115 OpenUrl CrossRef ↵ Ulrich , W. , Zalewski , M. ( 2006 ) Abundance and co-occurrence patterns of core and satellite species of ground beetles on small lake islands . Oikos . Vol. 114 , No. 2 . P. 338 – 348 . doi: 10.1111/j.2006.0030-1299.14773.x OpenUrl CrossRef Web of Science ↵ Veech , J.A. ( 2006 ) A probability-based analysis of temporal and spatial co-occurrence in grassland birds . Journal of Biogeography . Vol. 33 , Iissue 12. P. 2145 – 2153 , doi: 10.1111/j.1365-2699.2006.01571.x OpenUrl CrossRef Web of Science ↵ Veech , J.A. ( 2013 ) A probabilistic model for analysing species co-occurrence . Global Ecology and Biogeography . Vol. 22 , Issue 2 . P. 252 – 260 . doi: 10.1111/j.1466-8238.2012.00789.x OpenUrl CrossRef ↵ Wassersteina , R.L. , Lazar , N.A. ( 2016 ) The ASA’s statement on p-values: context, process, and purpose . The American Statistician . Vol. 70 , Issue 2 . P. 129 – 133 . doi: 10.1080/00031305.2016.1154108 OpenUrl CrossRef ↵ Wootton , J. T. ( 1994 ). The nature and consequences of indirect effects on ecological communities . Annual Review of Ecology and Systematics . Vol. 25 . P. 443 – 466 . doi: 10.1146/annurev.es.25.110194.002303 OpenUrl CrossRef Web of Science View the discussion thread. Back to top Previous Next Posted August 06, 2024. Download PDF 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 Methodology for the identification and analysis of species associations in biological communities using statistical distributions 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 Methodology for the identification and analysis of species associations in biological communities using statistical distributions Dmitriy G. Seleznev , Alexander A. Prokin , Ekaterina M. Kurina bioRxiv 2024.08.02.606312; doi: https://doi.org/10.1101/2024.08.02.606312 Share This Article: Copy Citation Tools Methodology for the identification and analysis of species associations in biological communities using statistical distributions Dmitriy G. Seleznev , Alexander A. Prokin , Ekaterina M. Kurina bioRxiv 2024.08.02.606312; doi: https://doi.org/10.1101/2024.08.02.606312 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 Ecology Subject Areas All Articles Animal Behavior and Cognition (7646) Biochemistry (17728) Bioengineering (13917) Bioinformatics (42038) Biophysics (21489) Cancer Biology (18637) Cell Biology (25554) Clinical Trials (138) Developmental Biology (13403) Ecology (19941) Epidemiology (2067) Evolutionary Biology (24368) Genetics (15624) Genomics (22547) Immunology (17764) Microbiology (40475) Molecular Biology (17208) Neuroscience (88756) Paleontology (667) Pathology (2842) Pharmacology and Toxicology (4834) Physiology (7659) Plant Biology (15175) Scientific Communication and Education (2047) Synthetic Biology (4304) Systems Biology (9835) Zoology (2272)