Abstract

Our aim was to assess step by step the evolution of the herbaceous plant community in reclaimed landfill areas by network analysis techniques. We examined the network of interactions among the members of the plant community at two levels; global and local. The experiment was conducted at one landfill that contained three reclaimed sites of different ages after closure (2 years [R2], 9 years [R9] and 13 years [R13]), plus a neighbor seminatural grazed grassland (SM-NAT). In R2 and R9 sites, the plant community consisted of almost 20 species; among them legumes like Melitotus albus and strong competitors like the exotic Solanum elaeagnifolium. The high heterogeneity of the interspecific relationships, the limited number of influential species and the moderate network centrality, recorded in R2 and R9 sites, were signs of competition for habitat exploitation. In contrast, the higher number of species (around 30), the relatively low heterogeneity of interspecific relationships and the even distribution of influence among many species in R13 and SM-NAT sites indicated moderation of competition. The global network metrics in sites R13 and SM-NAT appeared to converge although the composition of the community did not. The local scale analysis revealed the coexistence of two separate and well organized faction/communities within the same site (Mediterranean and widespread species), while within most factions the analysis distinguished further sub-categories like those formed by species of poor habitats and species of disturbed habitats. The identification of critical scales provides powerful theoretical and practical means for conservation and restoration practices.

1. Introduction

To study ecological communities where no species act in isolation, it is vital to understand the network of interactions among species. On this issue, Proulx et al. ¹ suggested that the structure of assemblages of coexisting species can be effectively modeled as a network of nodes representing species and ties standing for interactions among species. Within this context, community ecology has recently paid special attention to network theory in order to analyze patterns of relationships among nodes of various sorts ^{2, 3} such as competitors, predators or collaborators. Considering the type of interactions, Deng et al. ⁴ distinguished three network types; antagonistic, competitive and mutualistic.

According to Stamou and Papatheodorou ⁵, network analysis effectively depicted changes in the topology of ties among soil microbial variables that occur during the advancement of ecological succession processes. Since the establishment of the field of succession ecology, succession has been modelled as the progressive and directional change in the structure and the dynamics of biological communities ⁶. Lately, it has been hypothesized that abiotic factors act as decisive modulators of the vegetation succession process, particularly during the early phases ⁷, while the invasion of alien species and surrounding biotic configurations would strongly influence the way that succession unfolds ⁸. The unfolding of the ecological succession process is controlled by two types of interactions among plant species; facilitation and inhibition. Some already established species facilitate or at least tolerate the establishment of the incoming species, whereas some other beginners inhibit the establishment of the next colonizers. Accordingly the community networks appear as complex hybrids depicting both mutualistic and competitive interactions. As a result, recent mainstream theoretical models tackle succession as a complex procedure leading to multiple equilibrimium states or unpredictable outcomes ^{9, 10}. The unpredictability of the outcomes is further support by the idea of species functional redundancy i.e. the fact that more than one species share similar roles in each community.

The primary objective of this study is to evaluate whether network analysis could be an effective tool to assess changes in the organization of the plant herbaceous community developing during the ecological succession course launched with the reclamation of landfill sites. Actually, we monitored the process of reclamation in terms of spatial coexistence of plant species. The experiment was conducted in a single landfill that contained three reclaimed sections of different ages after closure (2, 9 and 13 years), in addition to a neighbouring natural grassland. Because of the management history (grazing) of the wider area, the grassland was deemed representative of the final stage of the succession process and was used as a control area.

2. Materials and Methods

This study was conducted in 2011 at Tagarades landfill (40° 27’ 30’’ to 40° 28’ 30’’ N latitude, 23° 2’ 30’’ E longitude, 120-180 m altitude above sea level), which is located about 21 km to the south-southeast of the city centre of Thessaloniki, Greece (straight-line distance). The average annual temperature for the last 40 years was 15.7 °C and the average annual precipitation was 448.7 mm. The prevailing wind direction was north-west almost year-round at this site. The average monthly air humidity was 67%. All weather data were collected from the nearest weather station, which is located 5 km away from the landfill site (http://www.hnms.gr/hnms/greek/climatology/ climatology_region_diagrams_html?dr_city=Thessaloniki_Mikra). The area surrounding the Tagarades landfill site is mainly agricultural and urban, with semi-natural vegetation of dry grasslands and thickets of Quercus coccifera shrubs existing along the boundary.

The Tagarades landfill was active from 1981 to 2009, and received the major proportion of the waste produced by the city of Thessaloniki. The landfill may be divided into three sections, which were filled at different time intervals; specifically, the first section was saturated 13 years ago, the second section before nine years and the third section before two years. These sections are reported hereafter as R13, R9 and R2, respectively.

After each section of the landfill was saturated, it was covered by a layer of soil about 0.5-1.5 m depth. Then, geotextile and aluminate layers were superimposed in order to encapsulate wastes, and the closure of each landfill section was completed by adding a suitable surface soil layer (of about 1 m depth) for vegetation growth. For all three landfill sections, the soil was collected from the same area, thus large differences in the availability of native species in the seed bank were not expected. Furthermore, systems for irrigation and the collection of leachates and biogases were installed in all three sections.

In Greece analogous restoration sequences are not easily found. Consequently, the particular sequence of Tagarades offered the unique opportunity to study the unfolding of an ecological succession course in restored landfill sites.

In total, 32 plots (1 m²) were sampled in the Tagarades landfill in June. In each of the sections R2, R9 and R13, eight 1 m² plots were sampled for the assessment of the herbaceous community composition. To avoid pseudoreplication given the constraint of sampling within single restored sites, the plots were non-randomly dispersed within each site to intentionally capture the variety in the number of species and the % cover. Intentional sampling techniques are non-probabilistic methods that select a group of individuals with the purpose of meeting specific prescribed criteria ¹¹. Actually, sampling in distinct microsites in relation to the characteristics of vegetation ensures taking genuine replicates ¹². Additionally, Bowker et al. ¹³ suggested that intentional sampling has the advantage to take plots rich in the most abundant species, and at the same time containing a variety of species mixtures with different abundances. Thus, it makes easier to isolate the actual contribution of the different species to the structure and performance of the biocommunity.

Eight additional plots were intentionally sampled in the seminatural dry grassland (SM-NAT), which was located on land adjacent (50, 300 and 500 m away from R13, R9 and R2 respectively) to the landfill. Therefore, four study sites were delineated in the area; three restored sites within the landfill and one grazed grassland site. All vascular plants were recorded in each plot, and their cover was estimated using the 9-grade refined Braun-Blanquet scale ¹⁴. The nomenclature of taxa follows Euro+Med ¹⁵ and Greuter et al. (1984-1989; ¹⁶).

Network analysis aims to represent and quantify structures and borrows from graph theory, the field of mathematics analyzing patterns of ties (edges) among entities (nodes). In this paper, nodes represent species and edges stand for the strength of co-existence of two species, while the network depicts how the associative spatial relationships are organized.

Network analysis allows for the simultaneous consideration and even depiction of both global and local aspects of the community organization (glocal approach which encapsulates the importance of both global and local structural attributes). The analysis began with the construction of the community matrix which herewith was based on species joint occurrence. Specifically, the probability of joint occurrence was calculated according to the overlap index of Pianka ¹⁷:

where: p_ik and p_jk are the proportional abundances of species i and j in the k sample respectively.

It should be noted that measures of species overlap reflect the intensity of their interactions: the more frequently the two taxa coming together, the higher their overlap and, hence, the more intense their relation. If the overlap of two taxa was below 10% of the recorded maximum, their co-occurrence was considered accidental, and the corresponding entry in the matrix was set to zero. To account for the heterogeneity in the intensity of pairwise interactions, we estimated the variability index CV% (i.e. the ratio of the standard deviation to the mean value expressed as a percentage). Moreover the binary version of the community matrix was also considered. The elements Q_ij of the latter were set either to 0 if the nodes i and j were not connected or otherwise to 1.

Community matrices were constructed separately for each study site, and were further analysed by using the network analysis software UCINET 6 ¹⁸. In the resulting graphs, the nodes standing for plant species and ties representing interactions among species were depicted employing a circular layout ¹⁹.

To describe the global aspects of the networks, we estimated parameters concerning a) the network cohesion, b) the tendency of nodes to form clusters and c) the network centrality (Table 1). Cohesion assesses the extent of connectedness of a network ^{20, 21}, clustering refers to the tendency of nodes towards bundling ²², while centrality measures the structural importance of a given node by evaluating the extent to which it occupies a more influential position than another ²³. In practice, centrality assesses how a given network is consolidated on certain nodes.

Cohesion is estimated by means of density, clustering coefficient and shortest path length. In this paper we estimated the density for weighted networks, i.e. the total of all ties’ intensities divided by the maximum number of possible ties, as well as the density of the binary version of the network ^{24, 25}. Finally, the shortest path length that denotes the minimum possible number of steps to move from one node to another was also assessed ²⁶.The overall clustering coefficient was estimated by averaging the local clustering coefficients of all nodes (see later) and it is used to assess the degree of modularity of the network architecture ⁴.

Further, we used the algorithm of UCINET 6 to partition nodes into mutually exclusive cohesive subgroups of nodes within the network (factions; ^{26, 27}). Also, the value of the coefficient E-I was calculated as the number of ties external to factions minus the number of ties that are internal to factions divided by the total number of ties ²⁸. Negative values of this parameter indicate numerical dominance of the internal interconnections and well organized structures.

We further calculated the number of the unordered transitive triplets and estimated transitivity. An unordered triplet of valued nodes a, b, c, form a transitive path if a is linked to b, b is linked to c and c is linked to a. Then, the unordered transitivity is the number of the transitive triplets divided by the sum of the transitive and non- transitive triplets.

Another partition of nodes to clusters comes from the consideration of their roles and positions (Table 1). The role or position of a node is established by means of its relational profile with other participants. In this paper we identified classes of structurally equivalent nodes. Two nodes of a network are considered structurally equivalent if the profiles of their relationships with others display enhanced similarity. It is expected that constraints operating on one member of a class operate also on the others ²⁶.

We further estimated eigenvector centrality to assess the structural importance of each node ^{29, 30, 31}. Given the community matrix, higher scores for this metric were assigned to nodes connected with well-connected neighbours rather than to nodes connected with poorly-connected neighbours ³². The average centrality of the whole network was expressed as a percentage of the potentially maximum centrality. Potentially maximum is the centrality of the “star network”, i.e. a network with the same number of nodes, where only one focal node is connected to all others, whereas the latter are only connected to the focal node. We further estimated the percentage contribution of the global relative to the local configurations to the structure of the network, which was accounted for by the ratio of the largest to the next largest eigenvalue ²⁶.

To describe the local aspects of the networks structure we focused on estimations relating to ego networks, i.e. the smallest possible network (Table 1). It includes a focal node (ego node) which is embedded in its relational setting consisted of the nodes with which the ego has direct connection (alters), the ties of the ego with alters, as well as the ties among all alters ²⁶. Given the grouping of nodes into factions, the consideration of homophily is of particular interest. Homophily refers to the extent to which egos tend to have ties with alters of the same faction as themselves, or in the words of Kossinets and Watts ³³ ‘the tendency of like to associate with like’. This tendency is often expressed as the ‘homophily principle’ and has been shown to be a key structural property in social networks. In this paper we measured ego’s homophily by employing the Yules’ Q metric. Its values range from -1 in case of perfect heterophily to +1 in case of perfect homophily ²⁶.

Finally, in order to get an insight on the dynamic of changing of networks structures during the unfolding of the succession process, we used the longitudinal algorithm of Ucinet 6 and measured new, lost and kept ties during the transition R2-R9, R9-R13 ³⁴. In addition, ties presented in SM-NAT but not in R13, those presented in R13 but not in SM-NAT, as well as ties presented in both R13 and SM-NAT were also enumerated.

3. Results

The species composition in each site is given in Annex 1 (nomenclature is according to Dimopoulos et al. ³⁶; to the same source is due also information about life history strategies, habitat, life form, and distribution).

Values of the global descriptors are included in Table 2. As shown, values of network density vary in a narrow range from 15.2 to 22.8% (weighted version) or between 34.7 and 56.1 (binary version). Enhanced CV% values were recorded in the younger sites R2 and R9 and they were accompanied by lower values estimated for the clustering coefficient. Estimations for transitivity vary from 0.349 to 0.476. Along all sites the values calculated for the mean length of the shortest path are roughly proportional to the Log of the number of nodes. Highest average eigenvector centrality was estimated for the R2 network, a precipitous decline was recorded in the R9 network, while moderate values were estimated in R13 and SM-NAT plots.

Despite the same number of nodes/species, radical compositional changes occurred during the transition from the R2 to the R9 structure. Actually, 158 new, 126 lost and only 6 kept ties were recorded. Even greater changes in the number of node/species, and the number of newly created ties were recorded also during the transition from the R9 to the R13 structure. By contrast, the number of the lost ties was low, while 52 ties were conserved. In comparison with the R13, the SM-NAT structure displays 348 different ties, 380 ties that were recorded in R13 but were not present in the SM-NAT site, while 84 ties were present in both.

The nodes of the R2 network were partitioned into three factions (Figure 1). Foremost, the nodes of the 1^st faction (blue symbols) stand for species inhabiting rather gentle habitats such as the Mediterranean grasslands (M. albus is exempted), displayed the highest density (43%), while all of them were perfectly homophiles. Among them the species M. albus is the more influential. Finally, with the exception of this latter, negative E-I values for the members of the 1^st faction were estimated (Table 3).

For most, the second faction consists of widespread early colonizers (red symbols). They are strong competitors and residents of ruderal habitats. The corresponding nodes displayed moderate density (33%), five out of nine species exhibited perfect homophily, while negative E-I values were estimated for all participants. The species H. incana, S. oleraceus and C. dactylon are the more influential.

The four species included in the 3^rd faction (green symbols) inhabited drier sites, the corresponding nodes displayed high density (43%), only one species showed perfect homophily while only for three out of four species negative E-I values were estimated. Moreover, four classes of structural equivalence were revealed. In the first (red colored labels) were included the taller and most aggressive competitors, in the second class (labels in black) participated widespread tolerant and weed species, the third class (labels in blue) is consisted essentially of Mediterranean species inhabiting poor and/or disturbed soils, whereas in the fourth class (labels in green) were included widespread species tolerant to extreme water conditions, either water excess or water deficit.

The features of this network show enhanced centrality, enhanced homophily and superiority of the within faction ties in combination with the highest CV% value (Table 3), illustrate this network as the most complex.

The nodes of this network were partitioned into three factions and few species display high values of eigenvector centrality (Figure 2). Foremost, the nodes of the first faction, (red symbols) stand for European, mostly Mediterranean species, displayed moderate density (35%), only four out of 10 exhibited perfect homophily, while negative E-I values were estimated for all of them.

The second faction include nodes (blue symbols) representing widespread cosmopolitan or paleotemperate species. In this faction the density equals to 30%, only two out of the nine species showed perfect homophily and only for M. minima a positive E-I value was estimated. Finally, in the third faction only the species V. hybrid was classified.

The nodes of the network were further partitioned into three classes of structural equivalence. The nodes of the first class (labels in red) stand for species of the 1^st (5 species) the 2^nd (5 species) and 3^rd (1 species) factions which were preferably established at well leached soils. In the 2^nd class (labels in green) were classified three (2^nd faction) and one (1^st faction) species able to colonize a large range of sites including disturbed Mediterranean habitats. The 3^rd class (labels in blue) represent one species distributed in worldwide xeric grasslands.

Finally, moderate centrality, existence of small world structures and superiority of the within faction ties, in combination with a relatively high value for the CV% indicate a rather complex network structure.

Most nodes of the R13 network (Figure 3) were classified into the first and principal faction (red symbols) and displayed modest density (0.43). Only six out of 24 species showed perfect homophily, seven more species showed enhanced homophily (Yule’s Q>0.80), while negative E-I values were estimated for all species. Moreover, four and three species constitute the 2^nd (blue symbols) and 3^rd (green symbols) minor factions, respectively. For all species of these factions positive E-I values were estimated. The above features in conjunction with the lowest heterogeneity recorded for the ties indicate a less complex network.

The nodes of the 1^st faction were further partitioned into four classes of structural equivalence. Constituting the 1^st class, five species of the 1^st faction, along with three species of the 2^nd faction (green labels) inhabit harsh Mediterranean habitats. Moreover nine species of the 1^st faction, foremost inhabiting dry and or disturbed habitats formed the second class. Eight Mediterranean species (blue labels), foremost inhabiting mainly ruderal sandy soils, constituted the 3^rd class. Finally, five species widely dispersed in grasslands formed the 4^th class.

The nodes of this network were partitioned into three factions (Figure 4). For most, the nodes of the 1^st faction, (blue symbols), stand for Mediterranean native plants, inhabiting ruderal habitats. These nodes displayed the highest density (48%), all of them showed perfect homophily (Yule’s Q=1), while, with the exception of C. teuniflorum and S. elaeagnifolium, for the rest members of the faction the E-I values were on the negative. Finally most influential were the species C. teuniflorum, and S. eleagnifolium.

In the second faction are included widespread and winter Mediterranean plants, dwelling xeric grasslands. The corresponding nodes (red symbols) exhibited low density equal to 38%, whereas only four out of 14 exhibited perfect homophily and only for T. angustifolium a positive E-I value was estimated. Finally, few peripheral nodes constituted the 3^rd fraction (green symbols).

The nodes of the factions were further subdivided into four groups of structurally equivalent species. The 1^st class (blue labels; density 69.6%) included taller and strong competing species preferentially inhabiting anthropogenic rough soils. In the 2^nd class (green labels; density 61.4%) were classified one species of the 1^st faction and five species of the 2^nd faction; these species are characterized as widely distributed species. The 3^rd class (red label, density 37.1%) was comprised of species preferring sandy pasture habitats, whereas the 4^th class (black label, density 42.6%) included species preferentially inhabiting rocky Mediterranean sites.

4. Discussion

The rather high number of nodes in the R2 network reflects a reasonably rich herbaceous community dominated by pioneer and alien species. Because the time since closure was not enough for other mechanisms of natural transportation to develop, the initial colonisation of the abandoned landfill should be associated with the richness of the seed bank in the cover soil transported from a neighbouring natural site. It is remarkable that up to 9 years after reclamation when a rather high complex network of interactions occurred, the richness of the community did not advance, probably due to the inhibition exerted by the species C. dactylon and S. elaeagnifolium since both have been reported to be efficient competitors ³⁷ and may affect the resistance of the community to change ³⁸. The influential role of these strong competitor species, in conjunction with the enhanced network density, the relatively modest centrality and the enhanced heterogeneity of ties’ intensity, indicate strong competition among plants for resources’ partition.

On the contrary, the net increase in species richness recorded in the R13 network, along with the less heterogeneity in the intensity of interactions (CV %) and the more or less even distribution of influence among nodes, indicate that the competition was calming down. This is further supported by the enhanced soil enzymatic activity in R13, compared to the rest reclaimed areas, which reflects an enhanced rate of nutrients cycling in this site ³⁹. Moreover, the high number of newly created ties and the low number of lost ties recorded during the transition from R9 to R13 are indicative of precipitous structural changes. However, these changes were not accompanied by increasing network connectedness, as expected according to suggestions by Parker et al. ⁴⁰.

The inspection of the eigenvector centrality values showed that R2 was the most centralized network. Actually, the structure of this network relies on enhanced dependence on few focal species. In such configurations the focal nodes operate as highly influential hubs that ensure the connectedness of the whole network ⁴¹. More to the point, Lusseau et al. ⁴² identified focal nodes as bottlenecks influencing decisively the whole network stability and connectivity and finally its robustness ^{25, 43}. In contrast, the R13 and the SM-NAT networks were the most decentralized. Networks of this type are considered more resilient; all nodes are more or less equally influential and the links among them are evenly distributed ⁴⁴.

Above remarks are in favor of the view that in our reclamation area the changes in the global structural characteristics of the networks (Table 2) converged progressively to the global characteristics of the SM-NAT network, probably reflecting analogous mechanisms of network formation ^{2, 45}. Actually, the successional route in Tagarades involved many species with a small number of interactions. Few facilitator legumes and competitor ruderals occupy focal position within the networks, maintaining multiple links to alters and thus determining the unfolding of the process in all stages.

However, the convergence process was not linear. To explain non-linear convergence in plant communities several hypotheses have been proposed suggesting either that the stochasticity of the succession process is determined by the site-specific environmental conditions ⁴⁶, or that it is induced by the idiosyncratic character of the colonisation sequence ^{47, 48} and the identity and density of the neighbouring plant species ⁴⁵. In this line was the study of Stadler et al. ⁴⁹ who found that the successional trajectories across sites were idiosyncratic. These statements are within the general theoretical assumption postulating that the succession course may follow a diversified array of potential restoration paths, possibly leading to vegetation compositions different from those in adjacent natural formations ⁵⁰. In this paper we proposed an additional reason for the non-linear convergence. Essentially, we believe that due to the idiosyncratic effect of legumes which create favorable conditions in the immediate vicinity and alien species which are reported as extremely forceful competitors that inhibit other species from thriving or even establishing ⁸, the occasional coexistence of different focal species in successive phases reinforced the stochasticity of the course of succession. Actually, the fact that different legumes and aliens held focal position in networks R2 (M. albus and C. dactylon, respectively) and R9 (M. minima and S. eleagnifolium, respectively), allows suggesting that these species may have differentially diverted (i.e. “side track”) the succession route.

In general, the findings of this study seem to fit to succession models proposed for arid and semi-arid environments. These models assume that restoration unfolds through distinct succession stages separated by thresholds ⁷. In Tagarades, the threshold might be set between R9 and R13, and could be identified with interspecific competition. With the caveat that it is risky to infer process from patterns, it seems likely that in the site R13 the interspecific competition decreased and it allowed overcoming the threshold. Then the network complexity dropped and the structural characteristics of the network converged with those of the semi-natural grassland, regardless of the identity of nodes.

From the consideration of the global density of the binary networks which exhibited medium values (34.7 to 56.1%), it is revealed slight superiority of direct relationships at the expense of indirect ones ²⁶, and presence of small world configurations. Actually, in all cases emerged characteristics such as the shortest path which is analogous to the Log of the number of nodes, the relatively high values of the global clustering coefficient, the values of the unordered transitivity ranging from 0.349 to 0.476 (values for transitivity >0.4 indicate modular structure; ²³) and the enhanced homophily, are indicative of cohesive small world structures ⁵¹. According to Telesford et al. ⁵¹ and Kurvers et al. ⁵² such small world configurations show that the participant species share habitat and nutrient resources readily. Moreover, following Telesford et al. ⁵¹ acknowledging that modularity is an essential attribute of the network robustness, we concluded that the modular biocommunities in reclaimed Tagarades sites are relatively resistant to change in the face of environmental perturbations ¹.

The identification of the network modularity offers interesting ecological insights. In this study, two scales of modularity were distinguished. In relation to the first scale with the exception of R13 network, the nodes of the rest networks were distributed into two main and one minor faction. The two main factions included Mediterranean and widespread species respectively. Obviously, this pattern reflects phylogenetic clustering of closely related species ⁵³. In relation to the second scale, most main factions were further subdivided into two ecologically equivalent classes of species reflecting ecological relatedness relative either to the exploitation of abiotic resources associated with habitat characteristics, (colonizers of poor sandy soils vs dwellers of disturbed soils) ^{53, 54} or their capability to cope with common environmental constraints (tolerants vs competitors species; ²⁶).

5. Conclusions

In this paper we sought to assess step by step the output of the successional route at the whole community level. In particular, we achieved to simultaneously perform analysis at different levels (gobal and local; glocal approach), as befits to modular topologies. Actually, second level analysis revealed the coexistence within the same site of two separate and well organized faction/communities (Mediterranean and widespread species), while within most factions a third level analysis distinguished further sub-categories like those formed by species of poor habitats and species of disturbed habitats. Obviously, the identification of such critical meso-scales, provides powerful theoretical and practical means for the validation of conservation or restoration practices.

References