More stories

  • in

    The impact of industrial agglomeration on urban green land use efficiency in the Yangtze River Economic Belt

    Research areaThe YREB covers Shanghai, Jiangsu, Zhejiang, Anhui, Jiangxi, Hubei, Hunan, Chongqing, Sichuan, Guizhou, and Yunnan. It includes the Yangtze River Delta urban agglomerations (YRDUA), Yangtze River midstream urban agglomeration (YRMUA), and Chengdu-Chongqing urban agglomeration (CCUA). With a regional area of 2.05 million km2, the YREB runs through the eastern, central and western regions in China32. In 2019, the total GDP of YREB is 45.8 trillion yuan, accounting for 46.2% of the national GDP. The YREB plays a pivotal strategic support and leading role in the overall situation of stable economic growth in China33. At the same time, the contradiction between the shortage of land resources and economic growth in the YREB is very prominent. Therefore, this paper selects 107 cities in YREB as the research sample. The specific geographic locations are shown in Fig. 2. This article uses ARCGIS 10.2 version to draw the map. The URL link is http://demo.domain.com:6080/arcgis/services.Figure 2The geographic location of the YREB in China.Full size imageResearch methodsGlobal Malmquist–Luenberger indexUGLUE refers to the effective utilization degree of land elements under certain input of other elements. The green utilization of urban land mainly comes from three aspects: first, improve the utilization intensity of the existing actual input land, that is, increase the input intensity of other elements of the unit land area. Second, reduce the input of land in the production process to avoid excessive waste of land. Third, promote the optimal allocation of land elements among production units. Technical efficiency refers to the maximum degree that all factor inputs need to expand or shrink in equal proportion when all production units reach the production frontier. However, for production units with high technical efficiency, the factor allocation structure may not be reasonable. The land factors may still have the problem of under-input or over-input, resulting in the reduction of UGLUE.Pastor and Lovell34 proposed a global index, which uses all the inspection periods of each decision-making unit as a benchmark to construct the production frontier. According to the current benchmark construction period t, the production possibility set reference set is defined as follows:$$P_{C}^{t} (x^{t} ) = left{ {left. {(y^{t} ,b^{t} )} right|x^{t} {kern 1pt} can{kern 1pt} , produce{kern 1pt} , b^{t} ,y^{t} } right}$$
    (1)
    The global benchmark is defined as: (P_{G} = P_{C}^{1} , cup ,P_{C}^{2} , cup , cdots ,P_{C}^{t}), The subscripts C and G represent the current benchmark and the global benchmark respectively. The ML index of decision-making unit i is calculated based on the current reference benchmark:$$ML^{S} (x^{t} ,y^{t} ,b^{t} ,x^{t + 1} ,y^{t + 1} ,b^{t + 1} ) = frac{{1 + D_{C}^{S} (x^{t} ,y^{t} ,b^{t} )}}{{1 + D_{C}^{S} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} )}}$$
    (2)
    Among them, the superscript S indicates two adjacent periods, t period and t + 1 period. The subscript C indicates the current benchmark, which is a simplified directional distance function. (ML^{s} > 1), indicates that the productivity increases. (ML^{s} < 1), indicates that the productivity decreases.According to Hofmann et al.35, the GMLI is defined as follows:$$GMLI^{t,t + 1} (x^{t} ,y^{t} ,b^{t} ,x^{t + 1} ,y^{t + 1} ,b^{t + 1} ) = frac{{1 + D_{G}^{T} (x^{t} ,y^{t} ,b^{t} )}}{{1 + D_{G}^{T} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} )}}$$ (3) Among them, (D_{G}^{T} (x,y,b) = max left{ {alpha |(y - alpha y,b - alpha b) in P_{G} (x)} right}). (GMLI^{t,t + 1} > 1) indicates that the productivity has increased. (GMLI^{t,t + 1} < 1) indicates that the productivity decreases. The GMLI is further broken down as follows:$$begin{aligned} & GMLI^{t,t + 1} (x^{t} ,y^{t} ,b^{t} ,x^{t + 1} ,y^{t + 1} ,b^{t + 1} ) = frac{{1 + D_{G}^{T} (x^{t} ,y^{t} ,b^{t} )}}{{1 + D_{G}^{T} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} )}} \ & quad = frac{{1 + D_{G}^{t} (x^{t} ,y^{t} ,b^{t} )}}{{1 + D_{G}^{t + 1} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} )}} times left[ {frac{{(1 + D_{G}^{T} (x^{t} ,y^{t} ,b^{t} ))/(1 + D_{C}^{T} (x^{t} ,y^{t} ,b^{t} ))}}{{(1 + D_{G}^{T} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} ))/(1 + D_{C}^{T} (x^{t + 1} ,y^{t + 1} ,b^{t + 1} ))}}} right] \ & quad = frac{{TE^{t + 1} }}{{TE^{t} }} times left( {frac{{BPG_{t + 1}^{t + 1} }}{{BPG_{t}^{t + 1} }}} right) = EC_{t}^{t + 1} times BPC_{t}^{t + 1} \ end{aligned}$$ (4) Among them, TE is the change of technological progress. EC is the change of technological efficiency. The change of technological progress reflects the change of the highest technical level. The improvement of the highest technical level often requires the introduction and innovation of advanced technology, which often requires a large amount of investment. The change of technical efficiency reflects the gap with the highest technical level. Narrowing the gap with the highest technical level often requires improvements in internal management and governance structures. (BPG_{t}^{t + 1}) is the “best practitioner gap” between the current period and overall technological frontier. (BPC_{t}^{t + 1}) measures the changes in the “best practitioner gap” between two periods (technological changes). (BPC_{t}^{t + 1} , > , 1 ,) indicates technological progress. (BPC_{t}^{t + 1} < 1) indicates technology regress.Econometric techniques of industrial agglomeration on UGLUEIn recent years, many scholars used the traditional SPM for empirical analysis, which is a basic measurement model suitable for panel data. Therefore, this article firstly uses the traditional SPM to analyze the impact of industrial agglomeration on UGLUE. The formula is:$$begin{aligned} ln UGLUE_{it} & = alpha_{0} + alpha_{1} ln RZI_{it} + alpha_{2} ln RZI_{it} *ln RZI_{it} + alpha_{3} ln RDI_{it} + alpha_{4} ln EC_{it} \ & quad + alpha_{5} ln GDP_{it} + alpha_{6} ln TEC_{it} + alpha_{7} ln ROAD_{it} + alpha_{8} ln GOV_{it} + varepsilon_{it} \ end{aligned}$$ (5) Among them, ε is the disturbance term. i represents the city, i in this paper involves 107 cities in YREB. t represents the time, and the range of t in this paper is from 2007 to 2016. UGLUE is the explained variable, which represents the UGLUE. RZI and RDI are explanatory variables, representing industrial specialization agglomeration and industrial diversification agglomeration. EC is the industrial structure. GDP is the level of economic development. TEC is the level of technology. ROAD is the level of infrastructure. GOV is the degree of government intervention. (alpha_{1}) to (alpha_{8}) is the coefficient of each variable.Formula (5) assumes that the UGLUE changes with the changes of various influencing factors in the current period. That is, there is no time lag effect. But in reality, land use often has a time lag effect. The previous level has a non-negligible impact on the current results. Therefore, this paper selects the dynamic panel model for empirical analysis. However, there is often a two-way causal relationship between industrial agglomeration and UGLUE, which may cause endogenous bias. For example, cities with higher UGLUE levels tend to have higher levels of economic development, which promotes industrial agglomeration in this city. Therefore, this paper adopts the method of system GMM for regression analysis of dynamic panel model. Compared with mixed OLS, system GMM can make full use of sample information, select appropriate lag terms as instrumental variables36. It can effectively solve the endogeneity problem between industrial agglomeration and UGLUE. Based on the above analysis, this paper introduces the first-order lag term of UGLUE on the basis of formula (5). The DPM is as follows:$$begin{aligned} ln UGLUE_{it} & = beta_{0} + tau ln UGLUE_{i(t - 1)} + beta_{1} ln RZI_{it} + beta_{2} ln RZI_{it} times ln RZI_{it} + beta_{3} ln RDI_{it} \ & quad + beta_{4} ln EC_{it} + beta_{5} ln GDP_{it} + beta_{6} ln TEC_{it} + beta_{7} ln ROAD_{it} + beta_{8} ln GOV_{it} + varepsilon_{{{text{it}}}} \ end{aligned}$$ (6) Among them, (tau) is the first-order lag coefficient of UGLUE, reflecting the time lag effect of UGLUE.Variable descriptionExplained variableThe GMLI is used to measure the UGLUE of 107 cities in YREB. According to existing research37, the following core evaluation index of UGLUE are selected (see Table 1). Regarding input indicators, we mainly choose land element input M, labor element input L, and capital element input K as input indicators. Regarding output indicators, we choose the added value of the secondary and tertiary industries in the municipal area as the expected output, and use the GDP deflator to convert it into a comparable value. At the same time, pollution indicators such as industrial wastewater emissions, industrial sulfur dioxide emissions, and industrial smoke (dust) emissions are selected as undesired output. Since the GMLI reflects the growth rate of UGLUE, this paper assumes that the GMLI in 2006 is 1, and then multiplies the calculated GMLI year by year to obtain the development level of UGLUE in each city from 2007 to 2016.Table 1 Input and output index.Full size tableExplanatory variablesIndustrial specialization index ZI is usually used to measure the specialization level of urban industries. The specialization index is represented by the share of the employment of the industry in the total employment of the city:$$ZI_{i} = Max_{j} (S_{ij} )$$ (7) Nextly, we use the relative specialization index to make a horizontal comparison of the specialization level between different cities:$$RZI_{i} = Max(S_{ij} /S_{j} )$$ (8) The most common measure of the level of industrial diversification is the Herfindahl–Hirschman Index (HHI). For city i, the HHI is the sum of the square sum of employment shares of all industries in the city. The diversification index is the reciprocal of the HHI:$$DZ_{i} = frac{1}{{sumlimits_{j} {S_{ij}^{2} } }}$$ (9) The expression of relative diversification index is as follows:$$RDI_{i} = {1 mathord{left/ {vphantom {1 {sumlimits_{j} {left| {S_{ij} - S_{j} } right|} }}} right. kern-0pt} {sumlimits_{j} {left| {S_{ij} - S_{j} } right|} }}$$ (10) Among them, Sij is the employment proportion of j industry in city i, and Sj is the proportion of the total employment of the national j industry. The greater value of RZI and RDI, the higher level of industrial specialization and diversification.Control variablesRegarding control variables, we choose the following variables as control variables.Industrial structure (EC): The continuous optimization of industrial structure promotes the improvement of UGLUE through three aspects: saving land, increasing land income and promoting the optimal allocation of land resources. This paper selects the added value of the tertiary industry as a percentage of GDP (take the logarithm) to express.Technological level (TEC): The higher the technological innovation level of a city is, the more it promotes the use of input elements and the transformation of innovation results, thereby improving the UGLUE. This paper selects the proportion of science and technology expenditure to fiscal expenditure (take the logarithm) to represent.Economic development level (GDP): The continuous economic development promote the rational allocation of various production factors and increase the level of urban land output, thereby improving the UGLUE. This paper selects GDP per capita (take the logarithm) to express.Road infrastructure level (ROAD): The continuous improvement of infrastructure reduces transportation costs and transaction costs, and promotes communication externalities between producers, consumers, and between producers and consumers. This paper selects the average road area per capita (take the logarithm) to express.Government behavior (GOV): Fiscal expenditure is an important means for the government to carry out macro-control. Appropriate fiscal expenditure makes up for market shortages, improves factor flow and resource allocation efficiency, and realizes positive economic externalities. This paper selects the proportion of fiscal expenditure to GDP (take the logarithm) to express. We can see the meaning of specific variables from Table 2.Table 2 The descriptive statistics of variables.Full size tableData sourceThe object of this thesis is the 107 cities in YREB from 2007 to 2016. The urban construction land area data comes from the "China Urban Construction Statistical Yearbook", and the rest of the index data all come from the "China City Statistical Yearbook". The URL link is https://www.cnki.net/. In order to maintain the integrity of the data, this article uses the average method to fill in the missing values. In addition, because Chaohu City began to be under the jurisdiction of Hefei City in 2011, Bijie City and Tongren City in Guizhou Province only became prefecture-level cities in 2011. The three cities and Pu'er City are taken from the sample to maintain the continuity of data. More

  • in

    Artificial intelligence for automated detection of large mammals creates path to upscale drone surveys

    Chapman, A. It’s okay to call them drones. J. Unmanned Veh. Syst. 2, iii–v (2014).Article 

    Google Scholar 
    Chabot, D., Hodgson, A. J., Hodgson, J. C. & Anderson, K. ‘Drone’: Technically correct, popularly accepted, socially acceptable. Drone Syst. Appl. 10, 399–405 (2022).Article 

    Google Scholar 
    Chabot, D. & Bird, D. M. Wildlife research and management methods in the 21st century: Where do unmanned aircraft fit in?. J. Unmanned Veh. Syst. 3, 137–155 (2015).Article 

    Google Scholar 
    Christie, K. S., Gilbert, S. L., Brown, C. L., Hatfield, M. & Hanson, L. Unmanned aircraft systems in wildlife research: Current and future applications of a transformative technology. Front. Ecol. Environ. 14, 241–251 (2016).Article 

    Google Scholar 
    Whitehead, K. & Hugenholtz, C. H. Remote sensing of the environment with small unmanned aircraft systems (UASs), part 1: A review of progress and challenges. J. Unmanned Veh. Syst. 2, 69–85 (2014).Article 

    Google Scholar 
    Barnas, A. et al. Evaluating behavioral responses of nesting lesser snow geese to unmanned aircraft surveys. Ecol. Evol. 8, 1328–1338 (2018).Article 

    Google Scholar 
    Mulero-Pázmány, M. et al. Unmanned aircraft systems as a new source of disturbance for wildlife: A systematic review. PLoS ONE 12, e0178448 (2017).Article 

    Google Scholar 
    Linchant, J., Lisein, J., Semeki, J., Lejeune, P. & Vermeulen, C. Are unmanned aircraft systems (UAS s) the future of wildlife monitoring? A review of accomplishments and challenges. Mammal Rev. 45, 239–252 (2015).Article 

    Google Scholar 
    Whitehead, K. et al. Remote sensing of the environment with small unmanned aircraft systems (UASs), part 2: Scientific and commercial applications. J. Unmanned Veh. Syst. 2, 86–102 (2014).Article 

    Google Scholar 
    Barasona, J. A. et al. Unmanned aircraft systems for studying spatial abundance of ungulates: Relevance to spatial epidemiology. PLoS ONE 9, e115608 (2014).Article 
    ADS 

    Google Scholar 
    Chrétien, L. P., Théau, J. & Ménard, P. Wildlife multispecies remote sensing using visible and thermal infrared imagery acquired from an unmanned aerial vehicle (UAV). Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 40, 241 (2015).Article 

    Google Scholar 
    Guo, X. et al. Application of UAV remote sensing for a population census of large wild herbivores—Taking the headwater region of the yellow river as an example. Remote Sens. 10, 1041 (2018).Article 
    ADS 

    Google Scholar 
    Hu, J., Wu, X. & Dai, M. Estimating the population size of migrating Tibetan antelopes Pantholops hodgsonii with unmanned aerial vehicles. Oryx 54, 101–109 (2020).Article 

    Google Scholar 
    Mulero-Pázmány, M., Stolper, R., Van Essen, L. D., Negro, J. J. & Sassen, T. Remotely piloted aircraft systems as a rhinoceros anti-poaching tool in Africa. PLoS ONE 9, e83873 (2014).Article 
    ADS 

    Google Scholar 
    Rey, N., Volpi, M., Joost, S. & Tuia, D. Detecting animals in African Savanna with UAVs and the crowds. Remote Sens. Environ. 200, 341–351 (2017).Article 
    ADS 

    Google Scholar 
    Schroeder, N. M., Panebianco, A., Gonzalez Musso, R. & Carmanchahi, P. An experimental approach to evaluate the potential of drones in terrestrial mammal research: A gregarious ungulate as a study model. R. Soc. Open Sci. 7, 191482 (2020).Article 
    ADS 

    Google Scholar 
    Su, X. et al. Using an unmanned aerial vehicle (UAV) to study wild yak in the highest desert in the world. Int. J. Remote Sens. 39, 5490–5503 (2018).Article 

    Google Scholar 
    Vermeulen, C., Lejeune, P., Lisein, J., Sawadogo, P. & Bouché, P. Unmanned aerial survey of elephants. PLoS ONE 8, e54700 (2013).Article 
    ADS 
    CAS 

    Google Scholar 
    Mallory, M. L. et al. Financial costs of conducting science in the Arctic: Examples from seabird research. Arct. Sci. 4, 624–633 (2018).Article 

    Google Scholar 
    Sasse, D. B. Job-related mortality of wildlife workers in the United States, 1937–2000. Wildl. Soc. Bull. 31, 1015–1020 (2003).
    Google Scholar 
    Loarie, S. R., Joppa, L. N. & Pimm, S. L. Satellites miss environmental priorities. Trends Ecol. Evol. 22, 630–632 (2007).Article 

    Google Scholar 
    IUCN. The IUCN Red List of Threatened Species. IUCN Red List of Threatened Species https://www.iucnredlist.org/en (2021).Mech, L. D. & Barber, S. M. A critique of wildlife radio-tracking and its use in National Parks: a report to the National Park Service. (2002).Patterson, C., Koski, W., Pace, P., McLuckie, B. & Bird, D. M. Evaluation of an unmanned aircraft system for detecting surrogate caribou targets in Labrador. J. Unmanned Veh. Syst. 4, 53–69 (2015).Article 

    Google Scholar 
    Hodgson, J. C. et al. Drones count wildlife more accurately and precisely than humans. Methods Ecol. Evol. 9, 1160–1167 (2018).Article 

    Google Scholar 
    Seymour, A. C., Dale, J., Hammill, M., Halpin, P. N. & Johnston, D. W. Automated detection and enumeration of marine wildlife using unmanned aircraft systems (UAS) and thermal imagery. Sci. Rep. 7, 1–10 (2017).Article 

    Google Scholar 
    COSEWIC. COSEWIC assessment and status report on the caribou (Rangifer tarandus) eastern migratory population, Torngat mountain population in Canada. (COSEWIC, Committee on the Status of Endangered Wildlife in Canada, 2017).Albawi, S., Mohammed, T. A. & Al-Zawi, S. Understanding of a convolutional neural network. in 2017 international conference on engineering and technology (ICET) 1–6 (IEEE, 2017).Gu, J. et al. Recent advances in convolutional neural networks. Pattern Recognit. 77, 354–377 (2018).Article 
    ADS 

    Google Scholar 
    Teuwen, J. & Moriakov, N. Convolutional neural networks. in Handbook of medical image computing and computer assisted intervention 481–501 (Elsevier, 2020).Corcoran, E., Winsen, M., Sudholz, A. & Hamilton, G. Automated detection of wildlife using drones: Synthesis, opportunities and constraints. Methods Ecol. Evol. 12, 1103–1114 (2021).Article 

    Google Scholar 
    Corcoran, E., Denman, S., Hanger, J., Wilson, B. & Hamilton, G. Automated detection of koalas using low-level aerial surveillance and machine learning. Sci. Rep. 9, 3208 (2019).Article 
    ADS 

    Google Scholar 
    Gray, P. C. et al. Drones and convolutional neural networks facilitate automated and accurate cetacean species identification and photogrammetry. Methods Ecol. Evol. 10, 1490–1500 (2019).Article 

    Google Scholar 
    Gray, P. C. et al. A convolutional neural network for detecting sea turtles in drone imagery. Methods Ecol. Evol. 10, 345–355 (2019).Article 

    Google Scholar 
    Peng, J. et al. Wild animal survey using UAS imagery and deep learning: modified Faster R-CNN for kiang detection in Tibetan Plateau. ISPRS J. Photogramm. Remote Sens. 169, 364–376 (2020).Article 
    ADS 

    Google Scholar 
    Borowicz, A. et al. Multi-modal survey of Adélie penguin mega-colonies reveals the Danger Islands as a seabird hotspot. Sci. Rep. 8, 3926 (2018).Article 
    ADS 

    Google Scholar 
    Francis, R. J., Lyons, M. B., Kingsford, R. T. & Brandis, K. J. Counting mixed breeding aggregations of animal species using drones: Lessons from waterbirds on semi-automation. Remote Sens. 12, 1185 (2020).Article 
    ADS 

    Google Scholar 
    Santangeli, A. et al. Integrating drone-borne thermal imaging with artificial intelligence to locate bird nests on agricultural land. Sci. Rep. 10, 1–8 (2020).Article 

    Google Scholar 
    Bowley, C., Mattingly, M., Barnas, A., Ellis-Felege, S. & Desell, T. An analysis of altitude, citizen science and a convolutional neural network feedback loop on object detection in unmanned aerial systems. J. Comput. Sci. 34, 102–116 (2019).Article 

    Google Scholar 
    Bowley, C., Mattingly, M., Barnas, A., Ellis-Felege, S. & Desell, T. Detecting wildlife in unmanned aerial systems imagery using convolutional neural networks trained with an automated feedback loop. in International Conference on Computational Science 69–82 (Springer, 2018).Delplanque, A., Foucher, S., Lejeune, P., Linchant, J. & Théau, J. Multispecies detection and identification of African mammals in aerial imagery using convolutional neural networks. Remote Sens. Ecol. Conserv. 8, 166–179 (2021).Article 

    Google Scholar 
    Eikelboom, J. A. J. et al. Improving the precision and accuracy of animal population estimates with aerial image object detection. Methods Ecol. Evol. 10, 1875–1887 (2019).Article 

    Google Scholar 
    Kellenberger, B., Marcos, D. & Tuia, D. Detecting mammals in UAV images: Best practices to address a substantially imbalanced dataset with deep learning. Remote Sens. Environ. 216, 139–153 (2018).Article 
    ADS 

    Google Scholar 
    Hooge, I. T. C., Niehorster, D. C., Nyström, M., Andersson, R. & Hessels, R. S. Is human classification by experienced untrained observers a gold standard in fixation detection?. Behav. Res. Methods 50, 1864–1881 (2018).Article 

    Google Scholar 
    Barnas, A. F., Darby, B. J., Vandeberg, G. S., Rockwell, R. F. & Ellis-Felege, S. N. A comparison of drone imagery and ground-based methods for estimating the extent of habitat destruction by lesser snow geese (Anser caerulescens caerulescens) in La Pérouse Bay. PLoS ONE 14, e0217049 (2019).Article 
    CAS 

    Google Scholar 
    Brook, R. K. & Kenkel, N. C. A multivariate approach to vegetation mapping of Manitoba’s Hudson Bay Lowlands. Int. J. Remote Sens. 23, 4761–4776 (2002).Article 

    Google Scholar 
    Shilts, W. W., Aylsworth, J. M., Kaszycki, C. A., Klassen, R. A. & Graf, W. L. Canadian shield. in Geomorphic Systems of North America vol. 2 119–161 (Geological Society of America Boulder, Colorado, 1987).Barnas, A. F., Felege, C. J., Rockwell, R. F. & Ellis-Felege, S. N. A pilot (less) study on the use of an unmanned aircraft system for studying polar bears (Ursus maritimus). Polar Biol. 41, 1055–1062 (2018).Article 

    Google Scholar 
    Ellis-Felege, S. N. et al. Nesting common eiders (Somateria mollissima) show little behavioral response to fixed-wing drone surveys. J. Unmanned Veh. Syst. 10, 1–4 (2021).
    Google Scholar 
    Barnas, A. F. et al. A standardized protocol for reporting methods when using drones for wildlife research. J. Unmanned Veh. Syst. 8, 89–98 (2020).Article 

    Google Scholar 
    Ren, S., He, K., Girshick, R. & Sun, J. Faster R-CNN: Towards real-time object detection with region proposal networks. Adv. Neural Inf. Process. Syst. 28, 91–99 (2016).
    Google Scholar 
    Chen, T., Xu, B., Zhang, C. & Guestrin, C. Training Deep Nets with Sublinear Memory Cost. ArXiv160406174 Cs (2016).Pinckaers, H. & Litjens, G. Training convolutional neural networks with megapixel images. ArXiv180405712 Cs (2018).Abadi, M. et al. TensorFlow: Large-scale machine learning on heterogeneous systems. (2015).Janocha, K. & Czarnecki, W. M. On loss functions for deep neural networks in classification. ArXiv Prepr. ArXiv170205659. (2017).Murata, N., Yoshizawa, S. & Amari, S. Learning curves, model selection and complexity of neural networks. Adv. Neural Inf. Process. Syst. 5, 607–614 (1992).
    Google Scholar 
    Hänsch, R. & Hellwich, O. The truth about ground truth: Label noise in human-generated reference data. in IGARSS 2019–2019 IEEE International Geoscience and Remote Sensing Symposium 5594–5597 (IEEE, 2019).Bowler, E., Fretwell, P. T., French, G. & Mackiewicz, M. Using deep learning to count albatrosses from space: Assessing results in light of ground truth uncertainty. Remote Sens. 12, 2026 (2020).Article 
    ADS 

    Google Scholar 
    Brack, I. V., Kindel, A. & Oliveira, L. F. B. Detection errors in wildlife abundance estimates from Unmanned Aerial Systems (UAS) surveys: Synthesis, solutions, and challenges. Methods Ecol. Evol. 9, 1864–1873 (2018).Article 

    Google Scholar 
    Jagielski, P. M. et al. The utility of drones for studying polar bear behaviour in the Canadian Arctic: Opportunities and recommendations. Drone Syst. Appl. 10, 97–110 (2022).Article 

    Google Scholar 
    Williams, P. J., Hooten, M. B., Womble, J. N. & Bower, M. R. Estimating occupancy and abundance using aerial images with imperfect detection. Methods Ecol. Evol. 8, 1679–1689 (2017).Article 

    Google Scholar 
    Link, W. A., Schofield, M. R., Barker, R. J. & Sauer, J. R. On the robustness of N-mixture models. Ecology 99, 1547–1551 (2018).Article 

    Google Scholar 
    Horvitz, D. G. & Thompson, D. J. A generalization of sampling without replacement from a finite universe. J. Am. Stat. Assoc. 47, 663–685 (1952).Article 
    MATH 

    Google Scholar 
    Corcoran, E., Denman, S. & Hamilton, G. New technologies in the mix: Assessing N-mixture models for abundance estimation using automated detection data from drone surveys. Ecol. Evol. 10, 8176–8185 (2020).Article 

    Google Scholar 
    Lunga, D., Arndt, J., Gerrand, J. & Stewart, R. ReSFlow: A remote sensing imagery data-flow for improved model generalization. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 14, 10468–10483 (2021).Article 
    ADS 

    Google Scholar 
    Fromm, M., Schubert, M., Castilla, G., Linke, J. & McDermid, G. Automated detection of conifer seedlings in drone imagery using convolutional neural networks. Remote Sens. 11, 2585 (2019).Article 
    ADS 

    Google Scholar 
    Velumani, K. et al. Estimates of maize plant density from UAV RGB images using Faster-RCNN detection model: Impact of the spatial resolution. Plant Phenomics 2021, 9824843 (2021).Article 
    CAS 

    Google Scholar 
    Hodgson, A., Peel, D. & Kelly, N. Unmanned aerial vehicles for surveying marine fauna: Assessing detection probability. Ecol. Appl. 27, 1253–1267 (2017).Article 

    Google Scholar 
    Ferguson, M. C. et al. Performance of manned and unmanned aerial surveys to collect visual data and imagery for estimating arctic cetacean density and associated uncertainty. J. Unmanned Veh. Syst. 6, 128–154 (2018).Article 

    Google Scholar 
    Zmarz, A. et al. Application of UAV BVLOS remote sensing data for multi-faceted analysis of Antarctic ecosystem. Remote Sens. Environ. 217, 375–388 (2018).Article 
    ADS 

    Google Scholar 
    Lyons, M. B. et al. Monitoring large and complex wildlife aggregations with drones. Methods Ecol. Evol. 10, 1024–1035 (2019).Article 

    Google Scholar  More

  • in

    Genetic structure and relatedness of juvenile sicklefin lemon shark (Negaprion acutidens) at Dongsha Island

    Dulvy, N. K., Sadovy, Y. & Reynolds, J. D. Extinction vulnerability in marine populations. Fish Fish. 4, 25–64 (2003).Article 

    Google Scholar 
    Fowler S. L. et al. Sharks, Rays and Chimaeras: The Status of the Chondrichthyan Fishes. IUCN/SSC Shark Specialist Group, Gland, Switzerland and Cambridge, UK (2005).Dulvy, N. K. et al. Extinction risk and conservation of the world’s sharks and rays. Elife 3, e00590 (2014).Article 

    Google Scholar 
    Lack M. & Sant G. Illegal, Unreported and Unregulated Shark Catch: A review of current knowledge and action. Department of the Environment, Water, Heritage and the Arts and TRAFFIC, Canberra http://www.traffic.org/fish/ (2008).Rose D.A. An Overview of World Trade in Sharks and Other Cartilaginous Fishes. TRAFFIC International, Cambridge, UK (1996).Lam, V. Y. & Sadovy, M. Y. The sharks of South East Asia–unknown, unmonitored and unmanaged. Fish Fish 12, 51–74 (2011).Article 

    Google Scholar 
    Kessel S.T. Investigation into the behaviour and population dynamics of the lemon shark (Negaprion brevirostris). Cardiff University (United Kingdom) (2010).Morrissey, J. F. & Gruber, S. H. Habitat selection by juvenile lemon sharks Negaprion brevirostris. Environ. Biol. Fishes 38, 311–319 (1993).Article 

    Google Scholar 
    Filmalter, J. D., Dagorn, L. & Cowley, P. D. Spatial behaviour and site fidelity of the sicklefin lemon shark Negaprion acutidens in a remote Indian Ocean atoll. Mari. Biol. 160, 2425–2436 (2013).Article 

    Google Scholar 
    DiBattista, J. D. et al. A genetic assessment of polyandry and breeding site fidelity in lemon sharks. Mol. Ecol. 17, 3337–3351 (2008).Article 

    Google Scholar 
    Wetherbee, B. M., Gruber, S. H. & Rosa, R. S. Movement patterns of juvenile lemon sharks Negaprion brevirostris within Atol das Rocas, Brazil: A nursery characterized by tidal extremes. Mar. Ecol. Prog. Seri. 343, 283–293 (2007).Article 
    ADS 

    Google Scholar 
    Feldheim, K. A. et al. Two decades of genetic profiling yields first evidence of natal philopatry and long-term fidelity to parturition sites in sharks. Mol. Ecol. 23, 110–117 (2014).Article 

    Google Scholar 
    Stevens J. D. et al. Diversity, abundance and habitat utilisation of sharks and rays: Final report to West Australian Marine Science Institute. CSIRO, editor. Hobart (2009).Schultz, J. K. et al. Global phylogeography and seascape genetics of the lemon sharks (genus Negaprion). Mol. Ecol. 17, 5336–5348 (2008).Article 
    CAS 

    Google Scholar 
    Mourier, J., Buray, N., Schultz, J. K., Clua, E. & Planes, S. Genetic network and breeding patterns of a sicklefin lemon shark (Negaprion acutidens) population in the Society Islands, French Polynesia. PLoS ONE 8, e73899 (2013).Article 
    ADS 
    CAS 

    Google Scholar 
    Speed, C. W. et al. Reef shark movements relative to a coastal marine protected area. Reg. Stud. Mar. Sci. 3, 58–66 (2016).
    Google Scholar 
    Huang, Z. Marine Species and Their Distribution in China’s Seas (Krieger Publishing Company, 2001).
    Google Scholar 
    Chang, C. W., Huang, C. S. & Wang, S. I. Species composition and sizes of fish in the lagoon of dongsha island (Pratas Island), Dongsha Atoll of the South China sea. Platax 2012, 25–32 (2012).
    Google Scholar 
    Pillans, R. D. et al. Long-term acoustic monitoring reveals site fidelity, reproductive migrations, and sex specific differences in habitat use and migratory timing in a large coastal shark (Negaprion acutidens). Front. Mar. Sci. 8, 616633 (2021).Article 

    Google Scholar 
    Daly-Engel, T. S. et al. Global phylogeography with mixed-marker analysis reveals male-mediated dispersal in the endangered scalloped hammerhead shark (Sphyrna lewini). PLoS ONE 7, e29986 (2012).Article 
    ADS 
    CAS 

    Google Scholar 
    Félix-López, D. G. et al. Possible female philopatry of the smooth hammerhead shark Sphyrna zygaena revealed by genetic structure patterns. J. Fish Biol. 94, 671–679 (2019).Article 

    Google Scholar 
    Nosal, A. P., Caillat, A., Kisfaludy, E. K., Royer, M. A. & Wegner, N. C. Aggregation behavior and seasonal philopatry in male and female leopard sharks Triakis semifasciata along the open coast of southern California, USA. Mar. Ecol. Prog. Ser. 499, 157–175 (2014).Article 
    ADS 

    Google Scholar 
    Jirik, K. E. & Lowe, C. G. An elasmobranch maternity ward: Female round stingrays Urobatis halleri use warm, restored estuarine habitat during gestation. J. Fish. Biol. 80(5), 1227–1245 (2012).Article 
    CAS 

    Google Scholar 
    Jacoby, D. M., Croft, D. P. & Sims, D. W. Social behaviour in sharks and rays: Analysis, patterns and implications for conservation. Fish Fish 13(4), 399–417 (2012).Article 

    Google Scholar 
    Su, S. H., Liu, S. Y. V., Liu, K. M. & Tsai, W. P. Development and characterization of novel microsatellite loci for an endangered hammerhead shark Sphyrna lewini by using shotgun sequencing. Taiwania 65(2), 261–263 (2020).
    Google Scholar 
    Dieringer, D. & Schlötterer, C. Microsatellite analyser (MSA): A platform independent analysis tool for large microsatellite data sets. Mol. Ecol. Notes 3, 167–169 (2003).Article 
    CAS 

    Google Scholar 
    Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. & Shipley, P. Micro-checker: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538 (2004).Article 

    Google Scholar 
    Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).Article 
    CAS 

    Google Scholar 
    Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Dominant markers and null alleles. Mol. Ecol. Notes 7, 574–578 (2007).Article 
    CAS 

    Google Scholar 
    Earl, D. A. & VonHoldt, B. M. Structure harvester: A website and program for visualizing structure output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361 (2012).Article 

    Google Scholar 
    Jakobsson, M. & Rosenberg, N. A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23(14), 1801–1806 (2007).Article 
    CAS 

    Google Scholar 
    Peakall, R. & Smouse, P. E. GenAlEx 6.5: Genetic analysis in excel population genetic software for teaching and research–an update. Bioinformatics 28, 2537–2539 (2012).Article 
    CAS 

    Google Scholar 
    Kamvar, Z. N., Tabima, J. F. & Grünwald, N. J. POPPR: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ 2, e281 (2014).Article 

    Google Scholar 
    Kalinowski, S. T., Wagner, A. P. & Taper, M. L. ML-Relate: A computer program for maximum likelihood estimation of relatedness and relationship. Mol. Ecol. Resour. 6, 576–579 (2006).Article 
    CAS 

    Google Scholar 
    Do, C. et al. NeEstimator v2: Re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol. Ecol. Resour. 14, 209–214 (2014).Article 
    CAS 

    Google Scholar 
    Oh, B. Z. et al. Contrasting patterns of residency and space use of coastal sharks within a communal shark nursery. Mar. Freshw. Res. 68, 1501–1517 (2017).Article 

    Google Scholar 
    McClelland J. Genetic Assessment of Breeding Patterns and Population Size of the Sicklefin Lemon Shark Negaprion acutidens in a Tropical Marine Protected Area: Implications for Conservation and Management (Doctoral dissertation, University of York) (2020).Compagno L. J .V. FAO species catalogue Sharks of the world: An annotated and illustrated catalogue of shark species known to date. FAO Fish. Synop. No. 125 Rome 4, 1–655 (1984).Stevens, J. D. Life-history and ecology of sharks at aldabra Atoll. Indian Ocean. Proc R Soc. B 222, 79–106 (1984).ADS 

    Google Scholar 
    Kool, J. T., Moilanen, A. & Treml, E. A. Population connectivity: Recent advances and new perspectives. Landsc. Ecol. 28, 165–185 (2013).Article 

    Google Scholar 
    Ruzzante, D. E. et al. Effective number of breeders, effective population size and their relationship with census size in an iteroparous species Salvelinus fontinalis. Proc. R Soc. B 283, 20152601 (2016).Article 

    Google Scholar 
    Van Wyngaarden, M. et al. Identifying patterns of dispersal, connectivity and selection in the sea scallop, Placopecten magellanicus, using RADseq-derived SNPs. Evol. Appl. 10, 102–117 (2017).Article 

    Google Scholar 
    Frankham, R., Bradshaw, C. J. A. & Brook, B. W. Genetics in conservation management: Revised recommendations for the 50/500 rules, Red list criteria and population viability analyses. Biol. Conserv. 170, 56–63 (2014).Article 

    Google Scholar 
    Pazmiño, D. A., Maes, G. E., Simpfendorfer, C. A., Salinas-de-León, P. & van Herwerden, L. Genome-wide SNPs reveal low effective population size within confined management units of the highly vagile Galapagos shark (Carcharhinus galapagensis). Conserv. Genet. 18, 1151–1163 (2017).Article 

    Google Scholar 
    Waples, R. S. & Do, C. Linkage disequilibrium estimates of contemporary Ne using highly variable genetic markers: A largely untapped resource for applied conservation and evolution. Evol. Appl. 3, 244–262 (2010).Article 

    Google Scholar 
    Dudgeon, C. L. & Ovenden, J. R. The relationship between abundance and genetic effective population size in elasmobranchs: An example from the globally threatened zebra shark Stegostoma fasciatum within its protected range. Conserv. Genet. 16, 1443–1454 (2015).Article 

    Google Scholar 
    Feldheim, K. A., Gruber, S. H. & Ashley, M. V. Population genetic structure of the lemon shark (Negaprion brevirostris) in the western Atlantic: DNA microsatellite variation. Mol. Ecol. 10, 295–303 (2001).Article 
    CAS 

    Google Scholar 
    Feldheim, K. A., Gruber, S. H. & Ashley, M. V. The breeding biology of lemon sharks at a tropical nursery lagoon. Proc. R. Soc. Lond. B 269, 1471–2954 (2002).Article 

    Google Scholar 
    Portnoy, D., McDowell, J. R., Thompson, K., Musick, J. A. & Graves, J. E. Isolation and characterization of five dinucleotide microsatellite loci in the sandbar shark, Carcharhinus plumbeus. Mol. Ecol. Notes 6, 431–433 (2006).Article 
    CAS 

    Google Scholar  More

  • in

    Acclimatization of a coral-dinoflagellate mutualism at a CO2 vent

    Steffen, W. Introducing the Anthropocene: The human epoch. Ambio 50, 1784–1787 (2021).Article 

    Google Scholar 
    Keys, P. W. et al. Anthropocene risk. Nat. Sustain. 2, 667–673 (2019).Article 

    Google Scholar 
    Bell, G. Evolutionary rescue and the limits of adaptation. Philos. Trans. R. Soc. B Biol. Sci. 368, 20120080 (2013).Article 

    Google Scholar 
    Byrne, M. & Przeslawski, R. Multistressor impacts of warming and acidification of the ocean on marine invertebrates’ life histories. Integr. Comp. Biol. 53, 582–596 (2013).Article 
    CAS 

    Google Scholar 
    Feely, R. A. et al. Impact of anthropogenic CO2 on the CaCO3 system in the oceans. Science 305, 362–366 (2004).Article 
    CAS 

    Google Scholar 
    Doney, S. C., Fabry, V. J., Feely, R. A. & Kleypas, J. A. Ocean acidification: the other CO2 problem. Ann. Rev. Mar. Sci. 1, 169–192 (2009).Article 

    Google Scholar 
    Hill, T. S. & Hoogenboom, M. O. The indirect effects of ocean acidification on corals and coral communities. Coral Reefs https://doi.org/10.1007/s00338-022-02286-z (2022).Biagi, E. et al. Patterns in microbiome composition differ with ocean acidification in anatomic compartments of the Mediterranean coral Astroides calycularis living at CO2 vents. Sci. Total Environ. 724, 138048 (2020).Article 
    CAS 

    Google Scholar 
    Chen, B. et al. Microbiome community and complexity indicate environmental gradient acclimatisation and potential microbial interaction of endemic coral holobionts in the South China Sea. Sci. Total Environ. 765, 142690 (2021).Article 
    CAS 

    Google Scholar 
    Palumbi, S. R., Barshis, D. J., Traylor-Knowles, N. & Bay, R. A. Mechanisms of reef coral resistance to future climate change. Science 344, 895–898 (2014).Article 
    CAS 

    Google Scholar 
    Wood, R. The ecological evolution of reefs. Annu. Rev. Ecol. Syst. 29, 179–206 (1998).Article 

    Google Scholar 
    Drake, J. L. et al. How corals made rocks through the ages. Glob. Chang. Biol. 26, 31–53 (2020).Article 

    Google Scholar 
    Stanley, G. D. Photosymbiosis and the evolution of modern coral reefs. Science 312, 857–858 (2006).Article 
    CAS 

    Google Scholar 
    Kitahara, M. V., Cairns, S. D., Stolarski, J., Blair, D. & Miller, D. J. A comprehensive phylogenetic analysis of the scleractinia (Cnidaria, Anthozoa) based on mitochondrial CO1 sequence data. PLoS One. 5, e11490 (2010).Article 

    Google Scholar 
    Dubinsky, Z. & Jokiel, P. Ratio of energy and nutrient fluxes regulates symbiosis between zooxanthellae and corals. Pac. Sci. 48, 313–324 (1994).
    Google Scholar 
    Falkowski, P. G., Dubinsky, Z., Muscatine, L. & Porter, J. W. Light and the bioenergetics of a symbiotic coral. Bioscience 34, 705–709 (1984).Article 
    CAS 

    Google Scholar 
    Frankowiak, K., Roniewicz, E. & Stolarski, J. Photosymbiosis in Late Triassic scleractinian corals from the Italian Dolomites. PeerJ 9, e11062 (2021).Article 

    Google Scholar 
    Davy, S. K., Allemand, D. & Weis, V. M. Cell biology of cnidarian-dinoflagellate symbiosis. Microbiol. Mol. Biol. Rev. 76, 229–261 (2012).Article 
    CAS 

    Google Scholar 
    Kremer, P. Ingestion and elemental budgets for Linuche unguiculata, a scyphomedusa with zooxanthellae. J. Mar. Biol. Assoc. UK. 85, 613–625 (2005).Article 

    Google Scholar 
    Welsh, D. T., Dunn, R. J. K. & Meziane, T. Oxygen and nutrient dynamics of the upside down jellyfish (Cassiopea sp.) and its influence on benthic nutrient exchanges and primary production. Hydrobiologia 635, 351–362 (2009).Article 
    CAS 

    Google Scholar 
    Muscatine, L., McCloskey, L. R. & Marian, R. E. Estimating the daily contribution of carbon from zooxanthellae to coral animal respiration. Limnol. Oceanogr. 26, 601–611 (1981).Article 
    CAS 

    Google Scholar 
    Ferrier‐Pagès, C. & Leal, M. C. Stable isotopes as tracers of trophic interactions in marine mutualistic symbioses. Ecol. Evol. 9, 723–740 (2019).Article 

    Google Scholar 
    Teixidó, N. et al. Ocean acidification causes variable trait shifts in a coral species. Glob. Chang. Biol. 26, 6813–6830 (2020).Article 

    Google Scholar 
    Fantazzini, P. et al. Gains and losses of coral skeletal porosity changes with ocean acidification acclimation. Nat. Commun. 6, 7785 (2015).Article 
    CAS 

    Google Scholar 
    Prada, F. et al. Coral micro- and macro-morphological skeletal properties in response to life-long acclimatization at CO2 vents in Papua New Guinea. Sci. Rep. 11, 19927 (2021).Article 
    CAS 

    Google Scholar 
    Kerrison, P., Hall-Spencer, J. M., Suggett, D. J., Hepburn, L. J. & Steinke, M. Assessment of pH variability at a coastal CO2 vent for ocean acidification studies. Estuar. Coast. Shelf Sci. 94, 129–137 (2011).Article 
    CAS 

    Google Scholar 
    Johnson, V. R., Russell, B. D., Fabricius, K. E., Brownlee, C. & Hall-Spencer, J. M. Temperate and tropical brown macroalgae thrive, despite decalcification, along natural CO2 gradients. Glob. Chang. Biol. 18, 2792–2803 (2012).Article 

    Google Scholar 
    Caroselli, E. et al. Low and variable pH decreases recruitment efficiency in populations of a temperate coral naturally present at a CO2 vent. Limnol. Oceanogr. 64, 1059–1069 (2019).Article 
    CAS 

    Google Scholar 
    González-Delgado, S. & Hernández, J. C. The importance of natural acidified systems in the study of ocean acidification: what have we learned? Adv. Mar. Biol. 80, 57–99 (2018).Article 

    Google Scholar 
    Capaccioni, B., Tassi, F., Vaselli, O., Tedesco, D. & Poreda, R. Submarine gas burst at Panarea Island (southern Italy) on 3 November 2002: A magmatic versus hydrothermal episode. J. Geophys. Res. 112, B05201 (2007).
    Google Scholar 
    Reggi, M. et al. Biomineralization in mediterranean corals: The role of the intraskeletal organic matrix. Cryst. Growth Des. 14, 4310–4320 (2014).Article 
    CAS 

    Google Scholar 
    Prada, F. et al. Ocean warming and acidification synergistically increase coral mortality. Sci. Rep. 7, 1–10 (2017).Article 

    Google Scholar 
    Goffredo, S. et al. Biomineralization control related to population density under ocean acidification. Nat. Clim. Chang. 4, 593–597 (2014).Article 
    CAS 

    Google Scholar 
    Wall, M. et al. Linking internal carbonate chemistry regulation and calcification in corals growing at a Mediterranean CO2 vent. Front. Mar. Sci. 6, 699 (2019).Article 

    Google Scholar 
    Zohary, T., Erez, J., Gophen, M., Berman-Frank, I. & Stiller, M. Seasonality of stable carbon isotopes within the pelagic food web of Lake Kinneret. Limnol. Oceanogr. 39, 1030–1043 (1994).Article 
    CAS 

    Google Scholar 
    Xu, S. et al. Spatial variations in the trophic status of Favia palauensis corals in the South China Sea: Insights into their different adaptabilities under contrasting environmental conditions. Sci. China Earth Sci. 64, 839–852 (2021).Article 

    Google Scholar 
    Horwitz, R., Borell, E. M., Yam, R., Shemesh, A. & Fine, M. Natural high pCO2 increases autotrophy in Anemonia viridis (Anthozoa) as revealed from stable isotope (C, N) analysis. Sci. Rep. 5, 1–9 (2015).Article 

    Google Scholar 
    Chen, B., Zou, D., Zhu, M. & Yang, Y. Effects of CO2 levels and light intensities on growth and amino acid contents in red seaweed Gracilaria lemaneiformis. Aquac. Res. 48, 2683–2690 (2017).Article 
    CAS 

    Google Scholar 
    Winters, G., Beer, S., Zvi, B., Brickner, I. & Loya, Y. Spatial and temporal photoacclimation of Stylophora pistillata: zooxanthella size, pigmentation, location and clade. Mar. Ecol. Prog. Ser. 384, 107–119 (2009).Article 

    Google Scholar 
    Fitt, W. K., McFarland, F. K., Warner, M. E. & Chilcoat, G. C. Seasonal patterns of tissue biomass and densities of symbiotic dinoflagellates in reef corals and relation to coral bleaching. Limnol. Oceanogr. 45, 677–685 (2000).Article 
    CAS 

    Google Scholar 
    Wangpraseurt, D., Larkum, A. W. D., Ralph, P. J. & Kühl, M. Light gradients and optical microniches in coral tissues. Front. Microbiol. 3, 1–9 (2012).Article 

    Google Scholar 
    Krief, S. et al. Physiological and isotopic responses of scleractinian corals to ocean acidification. Geochim. Cosmochim. Acta. 74, 4988–5001 (2010).Article 
    CAS 

    Google Scholar 
    Scucchia, F., Malik, A., Zaslansky, P., Putnam, H. M. & Mass, T. Combined responses of primary coral polyps and their algal endosymbionts to decreasing seawater pH. Proc. R. Soc. B Biol. Sci. 288, 20210328 (2021).Article 
    CAS 

    Google Scholar 
    Anthony, K. R. N., Connolly, S. R. & Willis, B. L. Comparative analysis of energy allocation to tissue and skeletal growth in corals. Limnol. Oceanogr. 47, 1417–1429 (2002).Article 

    Google Scholar 
    LaJeunesse, T. C. et al. Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts. Curr. Biol. 28, 2570–2580.e6 (2018).Article 
    CAS 

    Google Scholar 
    Howells, E. J. et al. Coral thermal tolerance shaped by local adaptation of photosymbionts. Nat. Clim. Chang. 2, 116–120 (2012).Article 

    Google Scholar 
    Brading, P. et al. Differential effects of ocean acidification on growth and photosynthesis among phylotypes of Symbiodinium (Dinophyceae). Limnol. Oceanogr. 56, 927–938 (2011).Article 
    CAS 

    Google Scholar 
    Takahashi, T., Broecker, W. S. & Langer, S. Redfield ratio based on chemical data from isopycnal surfaces. J. Geophys. Res. 90, 6907 (1985).Article 
    CAS 

    Google Scholar 
    Xu, Z. et al. Changes of carbon to nitrogen ratio in particulate organic matter in the marine mesopelagic zone: A case from the South China Sea. Mar. Chem. 231, 103930 (2021).Article 
    CAS 

    Google Scholar 
    Crawford, D. W. et al. Low particulate carbon to nitrogen ratios in marine surface waters of the Arctic. Glob. Biogeochem. Cycles. 29, 2021–2033 (2015).Article 
    CAS 

    Google Scholar 
    Kikumoto, R. et al. Nitrogen isotope chemostratigraphy of the Ediacaran and Early Cambrian platform sequence at Three Gorges, South China. Gondwana Res. 25, 1057–1069 (2014).Article 
    CAS 

    Google Scholar 
    DeNiro, M. J. & Epstein, S. Influence of diet on the distribution of carbon isotopes in animals. Geochim. Cosmochim. Acta 42, 495–506 (1978).Article 
    CAS 

    Google Scholar 
    Benavides, M., Bednarz, V. N. & Ferrier-Pagès, C. Diazotrophs: Overlooked key players within the coral symbiosis and tropical reef ecosystems? Front. Mar. Sci. 4, 10 (2017).Article 

    Google Scholar 
    Wannicke, N., Frey, C., Law, C. S. & Voss, M. The response of the marine nitrogen cycle to ocean acidification. Glob. Chang. Biol. 24, 5031–5043 (2018).Article 

    Google Scholar 
    Bourne, D. G., Morrow, K. M. & Webster, N. S. Insights into the coral microbiome: underpinning the health and resilience of reef ecosystems. Annu. Rev. Microbiol. 70, 317–340 (2016).Article 
    CAS 

    Google Scholar 
    Palladino, G. et al. Metagenomic shifts in mucus, tissue and skeleton of the coral Balanophyllia europaea living along a natural CO2 gradient. ISME Commun. 2, 65 (2022).Article 

    Google Scholar 
    Muscatine, L. et al. Stable isotopes (δ13C and δ15N) of organic matrix from coral skeleton. Proc. Natl Acad. Sci. 102, 1525–1530 (2005).Article 
    CAS 

    Google Scholar 
    Lesser, M. P. et al. Nitrogen fixation by symbiotic cyanobacteria provides a source of nitrogen for the scleractinian coral Montastraea cavernosa. Mar. Ecol. Prog. Ser. 346, 143–152 (2007).Article 
    CAS 

    Google Scholar 
    Alamaru, A., Loya, Y., Brokovich, E., Yam, R. & Shemesh, A. Carbon and nitrogen utilization in two species of Red Sea corals along a depth gradient: Insights from stable isotope analysis of total organic material and lipids. Geochim. Cosmochim. Acta. 73, 5333–5342 (2009).Article 
    CAS 

    Google Scholar 
    Lesser, M. P., Mazel, C. H., Gorbunov, M. Y. & Falkowski, P. G. Discovery of symbiotic nitrogen-fixing cyanobacteria in corals. Science 305, 997–1000 (2004).Article 
    CAS 

    Google Scholar 
    Lesser, M. P., Morrow, K. M., Pankey, S. M. & Noonan, S. H. C. Diazotroph diversity and nitrogen fixation in the coral Stylophora pistillata from the Great Barrier Reef. ISME J. 12, 813–824 (2018).Article 
    CAS 

    Google Scholar 
    Marcelino, V. R., Morrow, K. M., Oppen, M. J. H., Bourne, D. G. & Verbruggen, H. Diversity and stability of coral endolithic microbial communities at a naturally high pCO2 reef. Mol. Ecol. 26, 5344–5357 (2017).Article 
    CAS 

    Google Scholar 
    Rädecker, N., Pogoreutz, C., Voolstra, C. R., Wiedenmann, J. & Wild, C. Nitrogen cycling in corals: the key to understanding holobiont functioning? Trends Microbiol. 23, 490–497 (2015).Article 

    Google Scholar 
    Santos, H. F. et al. Climate change affects key nitrogen-fixing bacterial populations on coral reefs. ISME J. 8, 2272–2279 (2014).Article 

    Google Scholar 
    Olson, N. D., Ainsworth, T. D., Gates, R. D. & Takabayashi, M. Diazotrophic bacteria associated with Hawaiian Montipora corals: Diversity and abundance in correlation with symbiotic dinoflagellates. J. Exp. Mar. Bio. Ecol. 371, 140–146 (2009).Article 
    CAS 

    Google Scholar 
    Zheng, X. et al. Effects of ocean acidification on carbon and nitrogen fixation in the hermatypic coral Galaxea fascicularis. Front. Mar. Sci. 8, 644965 (2021).Article 

    Google Scholar 
    Lewis, E. & Wallace, D. Program developed for CO2 system calculations. Ornl/Cdiac-105 1–21 (1998).Dickson, A. G. & Millero, F. J. A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media. Deep Sea Res. Part A. Oceanogr. Res. Pap. 34, 1733–1743 (1987).Article 
    CAS 

    Google Scholar 
    Dickson, A. G. Thermodynamics of the dissociation of boric acid in potassium chloride solutions from 273.15 to 318.15 K. J. Chem. Eng. Data. 35, 253–257 (1990).Article 
    CAS 

    Google Scholar 
    Mehrbach, C., Culberson, C. H., Hawley, J. E. & Pytkowicx, R. M. Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure. Limnol. Oceanogr. 18, 897–907 (1973).Article 
    CAS 

    Google Scholar 
    Ivancic, I. & Degobbis, D. An optimal manual procedure for ammonia analysis in natural waters by the indophenol blue method. Water Res. 18, 1143–1147 (1984).Article 
    CAS 

    Google Scholar 
    Parson, T. R., Maita, Y. & Llli, C. M. A manual of chemical & biological methods for seawater analysis. (Elsevier, 1984). https://doi.org/10.1016/C2009-0-07774-5Schreiber, U. Pulse-Amplitude-Modulation (PAM) fluorometry and saturation pulse method: an overview. in Chlorophyll a Fluorescence 1367, 279–319 (Springer Netherlands, 2004).Grover, R., Maguer, J. F., Reynaud-Vaganay, S. & Ferrier-Pagès, C. Uptake of ammonium by the scleractinian coral Stylophora pistillata: Effect of feeding, light, and ammonium concentrations. Limnol. Oceanogr. 47, 782–790 (2002).Article 

    Google Scholar 
    Tremblay, P., Grover, R., Maguer, J. F., Hoogenboom, M. & Ferrier-Pagès, C. Carbon translocation from symbiont to host depends on irradiance and food availability in the tropical coral Stylophora pistillata. Coral Reefs. 33, 1–13 (2014).Article 

    Google Scholar 
    Pupier, C. A. et al. Productivity and carbon fluxes depend on species and symbiont density in soft coral symbioses. Sci. Rep. 9, 17819 (2019).Article 

    Google Scholar 
    Ritchie, R. J. Universal chlorophyll equations for estimating chlorophylls a, b, c, and d and total chlorophylls in natural assemblages of photosynthetic organisms using acetone, methanol, or ethanol solvents. Photosynthetica 46, 115–126 (2008).Article 
    CAS 

    Google Scholar 
    Goffredo, S., Arnone, S. & Zaccanti, F. Sexual reproduction in the Mediterranean solitary coral Balanophyllia europaea (Scleractinia, Dendrophylliidae). Mar. Ecol. Prog. Ser. 229, 83–94 (2002).Article 

    Google Scholar 
    Barshis, D. J. et al. Genomic basis for coral resilience to climate change. Proc. Natl Acad. Sci. 110, 1387–1392 (2013).Article 
    CAS 

    Google Scholar 
    Moore, R. B. Highly organized structure in the non-coding region of the psbA minicircle from clade C Symbiodinium. Int. J. Syst. Evol. Microbiol. 53, 1725–1734 (2003).Article 
    CAS 

    Google Scholar 
    LaJeunesse, T. C. & Thornhill, D. J. Improved resolution of reef-coral endosymbiont (Symbiodinium) species diversity, ecology, and evolution through psbA non-coding region genotyping. PLoS One. 6, e29013 (2011).Article 
    CAS 

    Google Scholar 
    LaJeunesse, T. C. et al. Revival of Philozoon Geddes for host-specialized dinoflagellates, ‘zooxanthellae’, in animals from coastal temperate zones of northern and southern hemispheres. Eur. J. Phycol. 57, 166–180 (2022).Article 

    Google Scholar 
    Anderson, M. J. PERMANOVA: a FORTRAN computer program for permutational multivariate analysis of variance. Wiley StatsRef: Statistics Reference Online (2005). More

  • in

    Evolution of snow algae, from cosmopolitans to endemics, revealed by DNA analysis of ancient ice

    Classification of snow algae in the ice core based on ITS2 sequencesWe used high-throughput sequencing to obtain DNA sequences of algae from 19 layers of an ice core drilled on a glacier in central Asia, dated from present time to 8000 years ago (Fig. 1 and Table S1). In total, 17,016 unique sequences (phylotypes) for the fast-evolving algal nuclear rDNA internal transcribed spacer 2 (ITS2) region were determined in the ice core, from which 290 OTUs were defined with ≥98% nt sequence identity among all OTUs.The ITS2 sequences were classified at the species level according to the genetic species concept based on secondary structural differences in the ITS2 region, which correlate with the boundaries of most biological species [38]. The ITS2 sequences from ice core samples were classified into 24 subgroups consisting of 17 chlorophycean, 5 trebouxiophycean, and 2 ulvophycean groups based on their secondary structures and BLASTn results (Fig. S1 and Tables S3–S4). The 17 subgroups of Chlorophyceae were subdivided into 10 subgroups of the Chloromonadinia clade, 1 subgroup of the Monadinia clade (recently treated as the genus Microglena [54]), 3 subgroups of the Reinhardtinia clade, 2 subgroups of the Stephanosphaerinia clade, and 1 subgroup corresponding to an unnamed group (which is related to Ploeotila sp. CCCryo 086-99) (for the clade names, see [55]). Although the Chloromonadinia clade contains several snow species belonging to Chloromonas or Chlainomonas, the 10 subgroups of the Chloromonadinia clade were considered to be Chloromonas. The 5 trebouxiophycean subgroups were composed of 2 subgroups of the Chlorella group, 1 subgroup of the Raphidonema group, 1 subgroup of the Trebouxia group, and 1 subgroup of the Neocystis group. The 2 subgroups of Ulvophyceae were closely related to the genus Chamaetrichon and Planophila, respectively. It is noted that Sanguina (‘Chlamydomonas’-snow group B [6]), Ancylonema, and Mesotaenium, which are snow algal genera found throughout the world [56, 57], were not detected in the ice core samples (Tables S3–S4).Global distribution of the Raphidonema groupTo understand the process by which snow algae form geographically specific population structures and how they migrate globally across the glaciers and snow fields, it is necessary to focus on the microbial species that inhabit the global cryosphere. Previous work elucidated that the Raphidonema group and ‘Chlamydomonas’-snow group B (Sanguina) are the cosmopolitans at both poles [6], but the latter was not detected in ice core samples examined in this study. Therefore, to elucidate the evolutionary history of the Raphidonema group, we further analyzed the ITS2 sequences obtained from the ice core sample as well as the glacier-surface samples from both poles [6] and from the mid-latitudes (samples from 10 sites, obtained in this study) (hereafter, surface samples; Table S2). Members of the Raphidonema group were detected in the older (deep core) layers of the ice core and at the glacier surface of central Asia (Fig. S1 and Tables S3–S4), as well as in the red snow samples from both poles [6]. In central Asia, the Raphidonema group was found in the Russian, Chinese, and Kyrgyz samples but was not detected in the Japanese and Tajik samples (Fig. S1 and Tables S3–S4). Combining these sequences yielded 893,649 reads and 22,389 unique sequences for subsequent detailed analysis (Tables S5–S6). The taxonomic composition of the Raphidonema communities differed among the mid-latitude, ice core, Arctic, and Antarctic samples as determined by PERMANOVA (Table S7). Most of the unique sequences in the Raphidonema group were consistent with an endemic distribution (Tables S8–S10). An average of 77% of the unique sequences of the Raphidonema group were endemic to a specific region (mid-latitude, 96%; Antarctic, 66%; Arctic, 79%), accounting for 40% of the total sequencing reads (mid-latitude, 77%; Antarctic, 74%; Arctic, 22%) (Fig. 2a, b and Tables S9–S10). This result suggested that most of the unique sequences are endemic, indicating that their dispersal has been limited to their respective regions [58,59,60,61].Fig. 2: Distribution types of the Raphidonema group obtained from each region and the ice core based on ITS2 unique sequences.Proportions of unique sequence and number of sequencing reads are shown. a Unique sequences from surface snow and ice-core samples. b Number of sequencing reads from surface snow and ice core samples. c Unique sequences from the indicated locations within the ice core. d Number of sequencing reads of the unique sequences from the indicated locations within the ice core.Full size imageNext, we analyzed the global distribution of the cosmopolitan phylotypes of the Raphidonema group, because a previous study analyzed their distribution only at the poles [6]. Only a limited number of unique sequences were distributed in all regions (mid-latitude, 1.4%; Antarctic, 5.6%; Arctic, 3.1%), accounting for a large proportion of the sequencing reads in polar regions but for only a small proportion in the mid-latitudes (mid-latitude, 2.8%; Antarctic, 20%; Arctic, 55%) (Figs. 2a, b, S2–S3, and Tables S9–S10). The distribution types of the Raphidonema group obtained from each region and the ice core were similar between the USEARCH and DADA2 analyses (Figs. 2, S4). In addition, we note that in ancient samples, post-mortem nt substitutions, such as cytosine to thymine, accumulate over many years of deposition [62], and these are not included in the DADA2 error model, which leads to the elimination of minor sequences in the DADA2 analysis. Therefore, we based our analysis on the results of the USEARCH unique sequences. These results suggested that only a few snow algae in the Raphidonema group were detected in samples from the mid-latitude regions.Snow algae of the Raphidonema group were detected in different ice core layers, corresponding to different time periods. The ice core records revealed that the distribution types of the Raphidonema group have not changed significantly for the last 8000 years, with p = 0.1924 based on a PERMANOVA between the newer (1800–2001 AD) and the older (6000–8000 years before present) layers (Fig. 2c, d). In ice core samples, 77% of the unique sequences of the Raphidonema group were detected only in the ice core samples, accounting for 23% of the total sequencing reads (Fig. S5). Although some of these unique sequences may be artifacts of the post-mortem nt substitution or sequencing errors, because we conducted sequence quality filtering and removed the majority of artifact sequences by removing the singleton clusters, most of the unique sequences in the ice core are not likely to be artifacts, but they could represent endemic phylotypes (Figs. 2a, b, S5).The cosmopolitan phylotypes were detected over a broad period as represented by ice core samples. They were present in approximately similar ratios in the newer and older layers (Fig. 2c, d). The cosmopolitan phylotypes were relatively abundant in the ice core samples (average, 4.0%; range, 0.2–13%), accounting for 13% (0.9–81% in the samples) of the total sequencing reads (Figs. 2c, d and S5).Microevolution of cosmopolitan and endemic phylotypesWe analyzed the evolutionary relationship between cosmopolitan and endemic phylotypes of the Raphidonema group among all snow surface and ice core samples. In total, 22,389 unique sequences of the Raphidonema group were clustered into 170 OTUs that were defined with ≥98% nt sequence identity among sequences within OTUs. The OTU sequences were subdivided into five subgroups (Groups A–E) based on phylogenetic analysis (Figs. S6–S11 and Tables S11–S12). Based on a previous study [63], Groups A–C and Group E were assigned to R. sempervirens and R. nivale, respectively, but Group D was not consistent with any species examined in that study (Fig. S6).The phylotypes were categorized into three subsets: the cosmopolitan phylotypes found at both poles and the mid-latitude regions; the multi-region phylotypes found in any two of the Antarctic, Arctic, and mid-latitude regions; and the endemic phylotypes found in only one of the three regions. Cosmopolitan phylotypes were found in Groups A, B, and C and accounted for 64.6% of the unique sequences. We then analyzed the dispersal of the three groups in detail.MJ networks [47] for the ITS2 sequences in each subgroup revealed that the cosmopolitan phylotypes were located at the center of the networks in Groups A and C that contained any types (endemics, multi-regions, and cosmopolitans) of the phylotypes, whereas the endemic phylotypes were considered to be derived from the cosmopolitan phylotypes (Figs. 3 and S12–S13). Moreover, the outgroup phylotypes were directly connected to the cosmopolitan phylotypes. These findings clearly showed that the cosmopolitan phylotypes were ancestral, whereas the endemic phylotypes were derived. In contrast, there were remarkable differences in the shape of the networks between Group B and the others (Groups A and C). In Group B, the Antarctic endemic phylotypes formed a distinct clade, and multi-region phylotypes seemed to be recently derived from this clade. In addition, the Arctic endemic phylotypes formed another distinct clade. These two Group-B clades split directly from a cosmopolitan phylotype (5.3% of the total sequencing reads). For Groups A and C, however, major portions of the total sequencing reads belonged to cosmopolitan phylotypes in Groups A (48.2%) and C (62.4%), and the endemic and multi-region phylotypes were directly connected to these major cosmopolitan phylotypes in a radial manner—the so-called “star-like” pattern [64]. These contrasting network shapes seem to have been formed as a consequence of the unique evolution of each of these groups. We also found that sequences from ice cores did not represent a basal position (Figs. 3 and S12–S13). This is because the haplotypes found in the modern samples have existed from times earlier than the ice core ages, due to the very small mutation numbers expected to have occurred since the ice core ages. Therefore, detected ice core ages were not included in the molecular evolution calculations of our demographic model. However, the phylogenetic networks themselves do not provide information on the evolutionary time scale. Hence, the ice core samples provide further direct evidence that Raphidonema, especially cosmopolitans belonging to this genus, persistently grew on snow and ice at least during the Holocene, and their ITS2 sequences have not changed over the last 8000 years.Fig. 3: Phylogenetic relationships among phylotypes of the Raphidonema groups.Phylotype networks for ITS2 sequences within Groups A (a), B (b), and C (c) of the Raphidonema group that include the cosmopolitan phylotypes in this study. The median-joining method was used. Circles indicate phylotypes; the size of each circle is proportional to the number of unique sequences. Each notch on the edges represents a mutation. Phylotypes are colored according to geographic region. The arrow represents the phylotype in the outgroup (see Fig. S6).Full size imageReferring to “ancestral” phylotypes as those having a longer history than other, more recently derived phylotypes, it is possible that individuals not closely related can share the same ancestral phylotype. In such cases, if genetically far-related individuals from various geographical regions share the same ancestral phylotype, they appear to be “cosmopolitan” (Fig. S14a). In order to distinguish between these “apparent cosmopolitans”, and “true cosmopolitans” that migrate globally, it is necessary to show that the cosmopolitan and endemic phylotypes have distinct demographic histories rather than being part of a continuous population sharing certain demographic dynamics (Fig. S14). Because phylotype networks are not useful for quantifying the rate(s) of microevolution, we used the coalescent model to quantify phylotype demographics [65]. As numerous phylotypes must be analyzed with this approach, we concentrated on statistical inference based on pairwise comparisons of phylotypes, for which the likelihood can be determined in a practical manner (see Materials and Methods). Histograms for the number of mismatched sites between two phylotypes chosen from a set of phylotypes, which will be called the pairwise mismatch distribution, are shown in Figs. 4 and S15. For Groups A and C, the distribution among cosmopolitans, multi-regions, and endemics was unimodal, in which the modes align from left to right with the order cosmopolitans, multi-regions, and endemics. Rogers and Harpending [48] noted that this “wave” propagation results from the expansion in size of a population, which leads to large mismatches, and the mode shifts to the right (see Fig. 2 of [48]). As time passes, the mode shifts to the left and eventually returns to the origin, i.e., representing a population that has not undergone an expansion event. Rogers and Harpending obtained an approximate solution for the wave and fitted the solution to human mitochondrial sequence data. We improved upon their method based on the coalescent model (see Materials and Methods) and applied it to the ITS2 sequence data for snow algae.Fig. 4: Mismatch distribution based on the number of pairwise differences in each distribution type in Raphidonema groups.The lines represent the observed number of pairwise differences in each distribution type (cosmopolitan, multi-region, endemic) within the Raphidonema Groups A (a), B (b) and C (c). Calculations were performed for all distribution types of Raphidonema Groups A and C, for which various cosmopolitan phylotypes were detected. On the other hand, calculations for only multi-region and endemic phylotypes were performed for Raphidonema group B, because no variation was found in cosmopolitan phylotypes.Full size imageFor Group A, when we fit the single demographic model to all phylotypes, the log-likelihood was –414,487. In contrast, when we fit the demographic model to each subset, that is, cosmopolitans, multi-regions, and endemics, separately, the log-likelihood was –341,964. Because the latter is larger than the former, we fit the model to each subset of phylotypes separately. For Group C, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately, the log-likelihood was –142,106, which is larger than the log-likelihood, –218,080, when we fit the single demographic model to all phylotypes. In contrast to Groups A and C then, we fit the single demographic model to all phylotypes of Group B because the log-likelihood, –196,070, was larger than the log-likelihood, –220,145, when we fit the demographic model to the cosmopolitans, multi-regions, and endemics separately. These results suggested that cosmopolitans, multi-regions, and endemics experienced different demographic histories in Groups A and C, whereas they had the same demographic history in Group B (Table S13). These results indicate the cosmopolitans in Group A and C are true cosmopolitans, whereas the those in Group B can be regarded as an apparent cosmopolitan.The ML estimates of (tau = 2ut_0), (theta _0 = 2N_0u), and (theta _1 = 2N_1u) are shown in Table S13 with standard deviation values. The population expanded t years ago, with the size before and after the expansion being represented by N0 and N1, respectively. The mutation rate (u) was assumed to be 7.9 × 10–8/ sequence/generation, and the generation interval was assumed to be 24 days (Materials and Methods). In Group A, for the cosmopolitans, the estimates of t, N0, and N1 were (33.8/(2 times 7.9) times 10^8 times {textstyle{{24} over {365}}} = 1.4 times 10^7) years, ((0.108 – 0.010)/(2 times 7.9) times 10^8 = (6.8 – 0.63) times 10^5), and ((0.217)/(2 times 7.9) times 10^8 = 1.4 times 10^6), respectively. In the same way, we computed estimates of t, N0, and N1 of other phylotypes and other groups (Table S14). For the endemics, the respective values were 9.2 × 106 years, 80, and 2.1 × 107, and the values were 4.6 × 106 years, 139, and 1.5 × 107 for the multi-regions. Taking into account the minimum and maximum ranges of the mutation rates per generation as well as the generation intervals, t for cosmopolitans was 3.6 × 106–4.0 × 107 years ago, and t for endemics was 2.3 × 106–2.6 × 107 years ago (Table S14). These results suggested that the cosmopolitans existed at least 1.4 × 107 years ago, and the endemics were derived from the cosmopolitans 9.2 × 106 years ago. The size of the endemics expanded 2.6 × 105-fold, which may have resulted from extensive dispersal. The multi-regions tended to mimic the endemics. Note that our demographic model was simplified to avoid overparameterization. In reality, considering the branching patterns of the MJ network, it is plausible that the endemic phylotypes have been repetitively and continuously derived from the cosmopolitans in multiple lineages—from 9.2 × 106 years ago to the present. In the same way, as for Group C, our results suggested that the cosmopolitan population expanded 3.9-fold ~3.2 × 106 years ago, and the endemics were derived from the cosmopolitans 1.9 × 105 years ago. The size of the endemics expanded 59-fold. In contrast to the phylotypes of Groups A and C, those of Group B experienced no significant expansion (Supplementary Results). In Groups A and C, the derived endemics (and multi-regions) expanded greatly as compared with the ancestral cosmopolitans (Table S14). These extraordinary expansions constitute evidence for local adaptation by the endemic/multi-region populations. In contrast, there was no evidence of local adaptation in Group B. The mismatch distribution of the entire Group B (multi-regions + endemics) showed a multimodal pattern (Fig. 4), which is present in the populations with stable sizes for a long period. When the populations finally reach equilibrium, the mismatch distributions show the exponential distribution [48]. Based on our ML estimates (Table S14), the historical population of Group B has been stable. More

  • in

    Half a century of rising extinction risk of coral reef sharks and rays

    Plaisance, L., Caley, M. J., Brainard, R. E. & Knowlton, N. The diversity of coral reefs: what are we missing? PLoS One. 6, e25026 (2011).Article 
    ADS 
    CAS 

    Google Scholar 
    Hoegh-Guldberg, O., Poloczanska, E. S., Skirving, W. & Dove, S. Coral reef ecosystems under climate change and ocean acidification. Front. Mar. Sci. 4, 158 (2017).Article 

    Google Scholar 
    Mora, C. et al. Global human footprint on the linkage between biodiversity and ecosystem functioning in reef fishes. PLoS Biol. 9, e1000606 (2011).Article 
    CAS 

    Google Scholar 
    Burke, L., Reytar, K., Spalding, M. & Perry, A. Reefs at Risk Revisited. 130 pp. (World Resources Institute, Washington, D.C., 2011).Hicks, C. C., Graham, N. A. J., Maire, E. & Robinson, J. P. W. Secure local aquatic food systems in the face of declining coral reefs. One Earth. 4, 1214–1216 (2021).Article 
    ADS 

    Google Scholar 
    Cinner, J. E. et al. Gravity of human impacts mediates coral reef conservation gains. PNAS 115, E6116–E6125 (2018).Article 
    CAS 

    Google Scholar 
    Eddy, T. D. et al. Global decline in capacity of coral reefs to provide ecosystem services. One Earth. 4, 1278–1285 (2021).Article 
    ADS 

    Google Scholar 
    Graham, N. A. J. et al. Human disruption of coral reef trophic structure. Curr. Biol. 27, 231–236 (2017).Article 
    CAS 

    Google Scholar 
    Sherman, C. S., Heupel, M. R., Moore, S. K., Chin, A. & Simpfendorfer, C. A. When sharks are away rays will play: effects of top predator removal in coral reef ecosystems. Mar. Ecol. Prog. Ser. 641, 145–157 (2020).Article 
    ADS 

    Google Scholar 
    Ruppert, J. L. W., Travers, M. J., Smith, L. L., Fortin, M.-J. & Meekan, M. G. Caught in the middle: combined Impacts of shark removal and coral loss on the fish communities of coral reefs. PLoS One. 8, e74648 (2013).Article 
    ADS 
    CAS 

    Google Scholar 
    Last, P. R. et al. Rays of the World. (CSIRO Publishing, 2016).Ebert, D. A., Dando, M. & Fowler, S. Sharks of the World. 2nd edn, 608 (Princeton University Press, 2021).Heupel, M. R., Lédée, E. J. I. & Simpfendorer, C. A. Telemetry reveals spatial separation of co-occurring reef sharks. Mar. Ecol. Prog. Ser. 589, 179–192 (2018).Article 
    ADS 

    Google Scholar 
    Heupel, M. R., Papastamatiou, Y. P., Espinoza, M., Green, M. E. & Simpfendorfer, C. A. Reef shark science – key questions and future directions. Front. Mar. Sci. 6, 12 (2019).Article 

    Google Scholar 
    Roff, G., Brown, C. J., Priest, M. A. & Mumby, P. J. Decline of coastal apex shark populations over the past half century. Commun. Biol. 1, 223 (2018).Article 

    Google Scholar 
    Williams, J. J., Papastamatiou, Y. P., Caselle, J. E., Bradley, D. & Jacoby, D. M. P. Mobile marine predators: an understudied source of nutrients to coral reefs in an unfished atoll. Proc. R. Soc. B. 285, 20172456 (2018).Article 

    Google Scholar 
    Heithaus, M. R., Wirsing, A. J. & Dill, L. M. The ecological importance of intact top-predator populations: a synthesis of 15 years of research in a seagrass ecosystem. Mar. Freshw. Res. 63, 1039–1050 (2012).Article 

    Google Scholar 
    Peel, L. R. et al. Stable isotope analyses reveal unique trophic role of reef manta rays (Mobula alfredi) at a remote coral reef. R. Soc. Open Sci. 6, 190599 (2019).Article 
    ADS 
    CAS 

    Google Scholar 
    O’Shea, O. R., Thums, M., van Keulen, M. & Meekan, M. Bioturbation by stingrays at Ningaloo Reef, Western Australia. Mar. Freshw. Res. 63, 189–197 (2012).Article 

    Google Scholar 
    Takeuchi, S. & Tamaki, A. Assessment of benthic disturbance associated with stingray foraging for ghost shrimp by aerial survey over an intertidal sandflat. Continental Shelf Res. 84, 139–157 (2014).Article 
    ADS 

    Google Scholar 
    Burkholder, D. A., Heithaus, M. R., Fourqurean, J. W., Wirsing, A. & Dill, L. M. Patterns of top-down control in a seagrass ecosystem: could a roving apex predator induce a behaviour-mediated trophic cascade? J. Anim. Ecol. 82, 1192–1202 (2013).Article 

    Google Scholar 
    Creel, S. & Christianson, D. Relationships between direct predation and risk effects. TRENDS Ecol. Evolution. 23, 194–201 (2008).Article 

    Google Scholar 
    Ward-Paige, C. A. et al. Large-scale absence of sharks on reefs in the greater-Caribbean: a footprint of human presence. PLoS One. 5, e11968 (2010).Article 
    ADS 

    Google Scholar 
    Espinoza, M., Cappo, M., Heupel, M. R., Tobin, A. J. & Simpfendorfer, C. A. Quantifying shark distribution patterns and species-habitat associations: implications of marine park zoning. PLoS One. 9, e106885 (2014).Article 
    ADS 

    Google Scholar 
    Graham, N. A., Spalding, M. D. & Sheppard, C. R. Reef shark declines in remote atolls highlight the need for multi-faceted conservation action. Aquat. Conserv.: Mar. Freshw. Ecosyst. 20, 543–548 (2010).Article 

    Google Scholar 
    Nadon, M. O. et al. Re-creating missing population baselines for Pacific reef sharks. Conserv. Biol. 26, 493–503 (2012).Article 

    Google Scholar 
    MacNeil, M. A. et al. Global status and conservation potential of reef sharks. Nature 583, 801–806 (2020).Article 
    ADS 
    CAS 

    Google Scholar 
    Dulvy, N. K. et al. Overfishing drives over one-third of all sharks and rays toward a global extinction crisis. Curr. Biol. 31, 1–15 (2021).Article 

    Google Scholar 
    Walls, R. H. L. & Dulvy, N. K. Eliminating the dark matter of data deficiency by predicting the conservation status of Northeast Atlantic and Mediterranean Sea sharks and rays. Biol. Conserv. 246, 108459 (2020).Article 

    Google Scholar 
    Yan, H. F. et al. Overfishing and habitat loss drives range contraction of iconic marine fishes to near extinction. Science Adv. 7, eabb6026, (2021).Butchart, S. H. M. et al. Using Red List Indices to measure progress towards the 2010 target and beyond. Philos. Trans. R. Soc. B 360, 255–268 (2005).Article 
    CAS 

    Google Scholar 
    Sherman, C. S. et al. Taeniura lymma. The IUCN Red List of Threatened Species, eT116850766A116851089 (2021). 10.2305/IUCN.UK.2021-1.RLTS.T116850766A116851089.enPacoureau, N. et al. Half a century of global decline in oceanic sharks and rays. Nature 589, 567–571 (2021).Article 
    ADS 
    CAS 

    Google Scholar 
    Cardeñosa, D. et al. Small fins, large trade: a snapshot of the species composition of low-value shark fins in the Hong Kong markets. Anim. Conserv. 23, 203–211 (2019).Article 

    Google Scholar 
    Haque, A. B. & Spaet, J. L. Y. Trade in threatened elasmobranchs in the Bay of Bengal, Bangladesh. Fish. Res. 243, 106059 (2021).Article 

    Google Scholar 
    Alcala, A. C. & Russ, G. R. A direct test of the effects of protective management on abundance and yield of tropical marine resources. ICES J. Mar. Sci. 47, 40–47 (1990).Article 

    Google Scholar 
    Serrano, A. et al. Effects of anti-trawling artificial reefs on ecological indicators of inner shelf fish and invertebrate communities in the Cantabrian Sea (southern Bay of Biscay). J. Mar. Biol. Assoc. U. Kingd. 91, 623–633 (2011).Article 

    Google Scholar 
    Cortés, E. Perspectives on the intrinsic rate of population growth. Methods Ecol. Evolution. 7, 1136–1145 (2016).Article 

    Google Scholar 
    McClenachan, L., Cooper, A. B. & Dulvy, N. K. Rethinking trade-driven extinction risk in marine and terrestrial megafauna. Curr. Biol. 26, 1–7 (2016).Article 

    Google Scholar 
    Tamburello, N., Cote, I. M. & Dulvy, N. K. Energy and the scaling of animal space use. Am. Naturalist 186, 196–211 (2015).Article 

    Google Scholar 
    Dulvy, N. K. et al. Challenges and priorities in shark and ray conservation. Curr. Biol. 27, R565–R572 (2017).Article 
    CAS 

    Google Scholar 
    Davidson, L. N. K. & Dulvy, N. K. Global marine protected areas to prevent extinctions. Ecol. Evolution. 1, 1–6 (2017).
    Google Scholar 
    Pauly, D., Zeller, D. & Palomares, M. L. D. Sea Around Us Concepts, Design and Data, (2021).Simpfendorfer, C. A. & Dulvy, N. K. Bright spots of sustainable shark fishing. Curr. Biol. 27, R83–R102 (2017).Article 

    Google Scholar 
    Booth, H., Squires, D. & Milner-Gulland, E. J. The mitigation hierarchy for sharks: a risk-based framework for reconciling trade-offs between shark conservation and fisheries objectives. Fish. Fish. 21, 269–289 (2019).Article 

    Google Scholar 
    Grorud-Colvert, K. et al. The MPA Guide: A framework to achieve global goals for the ocean. Science 373, eabf0861 (2021).Article 
    CAS 

    Google Scholar 
    Enright, S. R., Meneses-Orellana, R. & Keith, I. The Eastern Tropical Pacific Marine Corridor (CMAR): The emergence of a voluntary regional cooperation mechanism for the conservation and sustainable use of marine biodiversity within a fragmented regional ocean governance landscape. Front. Mar. Sci. 8, 674825 (2021).Article 

    Google Scholar 
    Chin, A., Kyne, P. M., Walker, T. I. & McAuley, R. B. An integrated risk assessment for climate change: analysing the vulnerability of sharks and rays on Australia’s Great Barrier Reef. Glob. Change Biol. 16, 1936–1953 (2010).Article 
    ADS 

    Google Scholar 
    Dwyer, R. G. et al. Individual and population benefits of marine reserves for reef sharks. Curr. Biol. 30, 480–489 (2020).Article 
    CAS 

    Google Scholar 
    Speed, C. W., Cappo, M. & Meekan, M. G. Evidence for rapid recovery of shark populations within a coral reef marine protected area. Biol. Conserv. 220, 308–319 (2018).Article 

    Google Scholar 
    Mizrahi, M. I., Diedrich, A., Weeks, R. & Pressey, R. L. A systematic review of the socioeconomic factors that influence how marine protected areas impact on ecosystems and livelihoods. Soc. Nat. Resour. 32, 4–20 (2019).Article 

    Google Scholar 
    IPBES. Global assessment report on biodiversity and ecosystem services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. 1148 (Bonn, Germany, 2019).Butchart, S. H. M. et al. Global biodiversity: indicators of recent declines. Science 328, 1164–1168 (2010).Article 
    ADS 
    CAS 

    Google Scholar 
    Hanh, T. T. H. & Boonstra, W. J. What prevents small-scale fishing and aquaculture households from engaging in alternative livelihoods? A case study in the Tam Giang lagoon, Viet Nam. Ocean Coast. Manag. 182, 104943 (2019).Article 

    Google Scholar 
    Ahmed, N., Troell, M., Allison, E. H. & Muir, J. F. Prawn postlarvae fishing in coastal Bangladesh: challenges for sustainable livelihoods. Mar. Policy. 34, 218–227 (2010).Article 

    Google Scholar 
    Prasetyo, A. P. et al. Shark and ray trade in and out of Indonesia: addressing knowledge gaps on the path to sustainability. Mar. Policy. 133, 104714 (2021).Article 

    Google Scholar 
    McClanahan, T., Polunin, N. & Done, T. Ecological states and the resilience of coral reefs. Conserv. Ecol. 6, 18 (2002).
    Google Scholar 
    Bellwood, D. R., Hughes, T. P. & Hoey, A. S. Sleeping functional group drives coral-reef recovery. Curr. Biol. 16, 2434–2439 (2006).Article 
    CAS 

    Google Scholar 
    Cinner, J. E. et al. Vulnerability of coastal communities to key impacts of climate change on coral reef fisheries. Glob. Environ. Change. 22, 12–20 (2012).Article 
    ADS 

    Google Scholar 
    Víe, J.-C., Hilton-Taylor, C. & Stuart, S. N. Wildlife in a Changing World – An analysis of the 2008 IUCN Red List of Threatened Species. 180 (Gland, Switzerland, 2009).Mace, G. M. et al. Quantification of extinction risk: IUCN’s system for classifying threatened species. Conserv. Biol. 22, 1424–1442 (2008).Article 

    Google Scholar 
    Sherley, R. B. et al. Estimating IUCN Red List population reduction: JARA – A decision-support tool applied to pelagic sharks. Conserv. Lett. 13, e12688 (2019).
    Google Scholar 
    IUCN Red List. Threats Classification Scheme (Version 3.2), (2021).Salafsky, N. et al. A standard lexicon for biodiversity conservation: unified classifications of threats and actions. Conserv. Biol. 22, 897–911 (2008).Article 

    Google Scholar 
    Moore, A. Chiloscyllium arabicum. The IUCN Red List of Threatened Species 2017, e.T161426A109902537 (2017). 10.2305/IUCN.UK.2017-2.RLTS.T161426A109902537.enSadovy de Mitcheson, Y. J. et al. Valuable but vulnerable: Over-fishing and under-management continue to threaten groupers so what now? Mar. Policy. 116, 103909 (2020).Article 

    Google Scholar 
    R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2021).Regression Models for Ordinal Data v. 2019.12.10 (CRAN, 2019).Econometric Tools for Performance and Risk Analysis v. 2.0.4 (2020).Dormann, C. F. et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46 (2013).Article 

    Google Scholar 
    Akinwande, M. O., Dikko, H. G. & Samson, A. Variance inflation factor: As a condition for the inclusion of suppressor variable(s) in regression analysis. Open J. Stat. 5, 754–767 (2015).Article 

    Google Scholar 
    Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociological Methods Res. 33, 261–304 (2004).Article 
    MathSciNet 

    Google Scholar 
    Plots Coefficients from Fitted Models v. 1.2.8 (2022).Fisheries and Aquaculture Software. FishStatJ – Software for Fishery and Aquaculture Statistical Time Series., (2020).Reynolds, R. W., Rayner, N. A., Smith, T. M., Stokes, D. C. & Wang, W. An improved in situ and satellite SST analysis for climate. J. Clim. 15, 1609–1625 (2002).Article 
    ADS 

    Google Scholar 
    NASA Ocean Biology (OB.DAAC). Mean annual sea surface chlorophyll-a concentration for the period 2009-2013 (composite dataset created by UNEP-WCMC). Data obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua Ocean Colour website (NASA OB.DAAC, Greenbelt, MD, USA), (2014).General Bathymetric Chart of the Oceans. GEBCO_2014 Grid. version 20150318. www.gebco.net (2015).XGBoost: A Scalable Tree Boosting System v. 1.4.1.1 (In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (pp. 785–794). New York, NY, USA: ACM, 2016).Elith, J., Leathwick, J. R. & Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 77, 802–813 (2008).Article 
    CAS 

    Google Scholar 
    ArcGIS Pro 2.7.0 (Environmental Systems Research Institute) (2020).Ferreira, L. C. & Simpfendorer, C. Galeocerdo cuvier. The IUCN Red List of Threatened Species 2019, e.T39378A2913541 (2019).Beta Regression v. 3.1-4 (2021).Butchart, S. H. et al. Improvements to the Red List Index. PLoS ONE. 2, e140 (2007).Article 
    ADS 

    Google Scholar 
    Sherman, C. S. et al. Half a century of rising extinction risk of coral reef sharks and rays, sammsherman27/CoralReefSharkRayIUCN: Data and Code Used in Sherman et al. Half a century of rising extinction risk of coral reef sharks and rays v1.0.0. https://doi.org/10.5281/zenodo.7267904 (2022). More

  • in

    Predicting cascading extinctions and efficient restoration strategies in plant–pollinator networks via generalized positive feedback loops

    The Campbell et al. model provides an excellent framework to identify species whose extinction leads to community collapse and species whose reintroduction can restore the community (see Fig. 2 for an illustration of these processes). Our first objective, finding the effect of species extinction on the rest of the species in an established community, is achievable using the concept of Logical Domain of Influence (LDOI)41; the LDOI represents the influence of a (set of) fixed node state(s) on the rest of the components in a system. In this section we first present our proposed method to calculate the LDOI for the Boolean threshold functions governing the Campbell et al. model of plant–pollinator community assembly. Then we verify that the simplified logical functions preserve the LDOI and hence can be implemented to further analyze the effect of extinction in plant–pollinator networks. Next, we address one of the main questions that motivated this study: Can stable motif driver set analysis facilitate the identification of keystone species? We discuss the identification of the driver sets of inactive stable motifs and motif groups and present the results of stabilizing these sets to measure the magnitude of the effect of species extinction on the communities. Lastly we discuss possible prevention and mitigation measures based on the knowledge acquired from driver sets of stable motifs and motif groups.Figure 2Illustration of species extinction and restoration in a hypothetical 6-species community. (a) The interaction network (on the left), and the maximal richness community possible for this network (the community with the most established species). Nodes highlighted with green represent established species. (b) The initial extinction of two species, po_1 and po_2 (left) and the community that results after cascading extinctions (right). Nodes highlighted with grey represent extinct species. (c) An intervention to restore pl_2 (left), which induces the restoration of further species, finally leading to a restored community with all the species present (right). The nodes highlighted with teal represent the restored species.Full size imageLDOI in the Boolean threshold modelThe LDOI concept was originally defined on Boolean functions expressed in a disjunctive prime form. Here we extend it to Boolean threshold functions. We implemented it as a breadth first search on the interaction network, as exemplified in Fig. 3. Assume that we want to find the LDOI of a (set of) node(s) (S_0={n_1,dots ,n_N}) and their specific fixed state (Q(S_0)={sigma _{n_1},dots ,sigma _{n_N}}). Starting from the set (S_0), the next set of nodes (S_1) that can acquire a fixed state due to the influence of (Q(S_0)) consists of the nodes that have an incoming edge from the nodes in the set (S_0) in the interaction network. The nodes in set (S_1) are the subject of the first search level. For each node (n_i in S_0) and (n^prime _i in S_1) we assume a “worst case scenario” (i.e., maximal opposition of the effect of (n_i) on (n^prime _i) from other regulators) to find the possible sufficiency relationships between the two. There are five cases:

    1.

    If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume that all other positive regulators of (n^prime _i) that have an unknown state (i.e., are not in (Q(S_0))) are inactive and all negative regulators of (n^prime _i) that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the active state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to LDOI((Q(S_0))).

    2.

    If (n_i) is a positive regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all other positive regulators of (n^prime _i) that have an unknown state are active and all negative regulators of (n^prime _i) that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the inactive state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

    3.

    If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=1) is a candidate for being sufficient for (sigma _{n^prime _i}=0). We assume all positive regulators of (n^prime _i) that have an unknown state are active and all other negative regulators of (n^prime _i) that that have an unknown state are inactive. If (sum _j W_{ij}le 0) under this assumption, then the active state of (n_i) is sufficient to deactivate (n^prime _i). The virtual node (sim n^prime _i) that corresponds to (sigma _{n^prime _i}=0) is added to LDOI((Q(S_0))).

    4.

    If (n_i) is a negative regulator of (n^prime _i), then (sigma _{n_i}=0) is a candidate for being sufficient for (sigma _{n^prime _i}=1). We assume all positive regulators of (n^prime _i) that have an unknown state are inactive and all other negative regulators of (n^prime _i) that that have an unknown state are active. If (sum _j W_{ij} > 0) under this assumption, then the inactive state of (n_i) is sufficient to activate (n^prime _i). The virtual node (n^prime _i) that corresponds to (sigma _{n^prime _i}=1) is added to the LDOI((Q(S_0))).

    5.

    If none of the past four sufficiency checks are satisfied, the node (n^prime _i) will be visited again in the next search levels.

    The second set of nodes that can be influenced, (S_2), are the nodes that have an incoming edge from the nodes in the set (S_1). The algorithm goes over these nodes in the second search level as described above. This search continues to all the levels of the search algorithm until all nodes are visited (possibly multiple times) and either acquire a fixed state and are added to the LDOI or their state will be left undetermined at the end of the algorithm. In Fig. 3, we illustrate this search to find the LDOI((sim )pl_1). The first search level is (S_1={)po_1, po_3(}); (sim )pl_1 is sufficient to deactivate po_3, but not po_1. As a result, (sim )po_3(in ) LDOI((sim )pl_1). This process continues until all levels are visited and at the end of the algorithm LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Figure 3Breadth first search of the interaction network to find the LDOI of a (set of) fixed note state(s) in Boolean threshold functions governing the dynamics of plant–pollinator networks. (a) An interaction network with five plants and 3 pollinators. (b) The breadth first search in the case of starting from the node state (sim )pl_1. The nodes with incoming edges from pl_1 make up (S_1={)po_1, po_3(}). The second sufficiency check is satisfied for node state (sim )po_3, as a result (sim )po_3(in ) LDOI((sim )pl_1). The same process is applied for node po_1, but none of the sufficiency checks are satisfied, so this node will be visited again later. The next level of the search consists of the nodes that have incident edges from (S_1), i.e., (S_2={)pl_2, pl_3, pl_4, pl_5(}). The second sufficiency check is satisfied for all of these nodes and they are all fixed to their inactive state in the LDOI((sim )pl_1). Lastly, we reach (S_3={)po_1, po_2(}). Node po_1 is reached again, and with both its positive regulators fixed to their inactive states the second sufficiency check is satisfied and node po_1 is fixed to its inactive state as well. The same holds for po_2 and hence LDOI((sim )pl_1()={sim )po_3, (sim )pl_2, (sim )pl_3, (sim )pl_4, (sim )pl_5, (sim )po_1, (sim )po_2 (}).Full size imageTo measure the accuracy of the simplification method originally introduced in28, we analyzed logical domains of influence in 6000 networks with 50–70 nodes. These networks are among the largest in our ensembles and have the most complex structures. We randomly selected (sets of) inactive node states, found their LDOIs using the Boolean threshold functions and the simplified Boolean functions, and compared the two resulting LDOIs. We used 8 single node states and 8 combinations of size 2 to 4 for each network. We found that in all cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions.Next, we analyzed (sets of) active node states and their LDOIs in the same ensembles of networks. Similar to the previous analysis, we used 8 single node states and 8 combinations of size 2 to 4 for each network. Our analysis shows that in 77.1% of the cases the LDOI calculated using the simplified Boolean functions matches the LDOI calculated using the Boolean threshold functions. In 22% of the cases the LDOI calculated from the simplified Boolean functions contains the LDOI calculated from the threshold functions, and it also contains extra active node states, overestimating the LDOI by 57.5% on average. These additional members of the LDOI result from the fact that the simplified Boolean functions contain fewer negative regulators than the threshold functions. The guiding principle of the simplification method is that the probability of (H(x)=1) conserves the probability of each node having an active state across all the states it can have. In contrast, the probability of the propagation of the active state is not necessarily preserved and tends to be higher in the simplified Boolean model; thus the LDOI of the active node states is overestimated in some cases.In the rest of the cases (about 1%), the LDOI calculated from the simplified Boolean functions does not fully capture the LDOI calculated from the threshold functions. This again is caused by the sparsification of the negative edges in the simplified Boolean functions. In the threshold functions, the activation of 4 or more negative regulators of a target node combined with one active positive regulator is sufficient to deactivate the target node, i.e., there might be inactive node states in the LDOI of a set of active node states. However, some of these negative regulators drop in the simplified Boolean model and the inactive state of the target node is not necessarily in the LDOI of the set of active node states in the simplified case. This is the rare mechanism by which the simplified model might underestimate the influence of active node states on the rest of the network.In the following section we are interested in analyzing the effect of species extinction on the established community, i.e., we look at the LDOI of (sets of) inactive node states. Observing that the influence of extinction of species is measured correctly in the simplified Boolean models, we conclude that these models can be utilized to further analyze the process of extinction and its ecological implications.Stable motif based identification of species whose loss leads to cascading extinctionsEach stable motif or motif group can have multiple driver sets; stabilization of each driver set leads to the stabilization of the whole motif or motif group. In plant–pollinator interaction networks, the stable motifs either represent a sub-community (when the constituent nodes stabilize in their active states) or the simultaneous extinction of all species in the group (when the constituent nodes stabilize to their inactive states). Stabilization of the nodes in the driver set of an inactive stable motif results in stabilization of all the nodes in the stable motif to their inactive state, i.e., cascading extinction of the constituent species.The knowledge gained from stable motif analysis and the network of functional relationships offers insight into the cascading effect of an extinction that constitutes a driver set of an inactive stable motif. The magnitude of this effect depends on (i) the number of nodes that the inactive stable motif contains and (ii) the number of virtual nodes (including motifs and motif groups) corresponding to inactive species that are logically determined by the stabilization of the inactive stable motif.To investigate the role of stable motifs in the study of species extinction in plant–pollinator networks, we simulated extinctions that drive inactive stable motifs in 6000 networks with the sizes of 50–70 nodes. We considered driver sets of size 1, 2, or 3, and implemented them by fixing the corresponding node(s) to its (their) inactive state. As a point of comparison, we also performed a “control” analysis using the same networks with the same size of initial extinction; however, the candidates of initial extinction are inactive node states that do not drive stable motifs or motif groups. Based on the properties of the drivers of stable motifs, one expects that following the extinction of driver species, cascading extinctions of other species follow, while the same does not necessarily hold for non-driver species. As a result, we expect to observe greater damage to the original community when driver species become extinct.We assume that the “maximal richness community”—the community (attractor) in which the largest number of species managed to establish—is the subject of species extinction. This maximal richness community results from the stabilization of all active stable motifs. All other attractors that have some established species contain a subset of all active stable motifs and thus will contain a subset of the species of the maximal richness community. While for a generic Boolean model with multiple attractors one expects that a perturbed version of the model also has multiple attractors, this specific perturbation of a plant–pollinator model (namely, extinction of species in the maximal richness community) has a single attractor. We prove this by contradiction. Assume there are two separate attractors in the perturbed model, which means that there is at least one node that has opposite states in these two attractors. Note that this bi-stability is the result of the perturbation and not a property of the original system as the maximal richness community (an attractor) is the starting point for the introduced extinction. Specifically, the inactive state of the extinct node has to lead to the stabilization of another node to its active versus inactive states in the two separate attractors. The only case in which the stabilization of an inactive node state can result in the stabilization of an active node state is if there is a negative edge from the former to the latter in the interaction network after simplification. Since the Boolean function in 2 is inhibitor dominant, the negative regulators that remain in the Boolean model must be in their inactive states in the maximal richness attractor. As they are already inactive (extinct), they are not candidates for extinction. The only nodes that are candidates for extinction are the ones that positively regulate other nodes; perturbing the system by fixing these candidates to their inactive states cannot lead to the active state of a target node. In conclusion, bi-stability is not possible.We found the new attractor of the system given the (combination of) inactive node state(s) using the the functions percolate_and_remove_constants() and trap_spaces() from the pyboolnet Python package. We quantify the effects of the initial extinction(s) on the maximal richness attractor by the percentage change in the number of active species, which we call damage percentage. Note that this choice of maximal richness community as the reference and starting point allows us to detect the cascading extinctions following the initial damage.In Fig. 4 the left column plots show the average damage percentage caused by the extinction of 1 (top panel), 2 (middle panel), or 3 (bottom panel) species that represent driver sets of stable motifs and motif groups, while the right column plots illustrate the average damage percentage as a result of the extinction of 1, 2 or 3 species that represent non-driver nodes. Comparing the two columns, one can notice that stabilization of the driver sets of stable motifs and motif groups leads to considerably larger damage to the communities. This is due to the fact that stabilization of driver sets ensures the stabilization of entire inactive stable motifs and motif groups and hence ensures cascading extinctions. Comparing the plots in the left column, we see that the larger the driver sets are, the larger the damage to the community becomes. This is because larger driver sets are more likely to stabilize larger stable motifs and motif groups. This figure illustrates the significance of stable motifs and their driver sets in the study of species extinction in plant–pollinator communities.Figure 4Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 6000 networks with the size of 50–70 nodes. To study the extinction of species we started from the maximal richness community, then we fixed the nodes that correspond to the focal species to the their inactive states. The original extinctions are excluded from the damage percentages. The left column plots show the average damage percentage caused to the maximal richness community by the extinction of a driver set of size 1 (top), 2 (middle), or 3 (bottom) of an inactive stable motif or motif group. For each network, we determined all the relevant driver sets of one stable motif or motif group, we performed the extinction and calculated the resulting damage, then we calculated the average damage percentage over all data points collected for the same network. The right column plots show the average damage percentage caused to the maximal richness community by the extinction of 1 (top), 2 (middle), and 3 (bottom) non-driver, randomly chosen nodes. Each time a randomly selected combination of non-driver nodes were the subject of simultaneous extinction until all combinations are explored and then we calculated the average damage percentage over all data points collected for each network. The number of networks that qualify for each of these 6 categories differ (e.g., some networks have a stable motif with a driver set of size 2 but no stable motif with a driver set of size 3). In the left column 5529, 3212, and 1980 networks and in the right column 5779, 5626, and 5423 networks qualified respectively. The red lines represent the mean value of all the presented data points in each plot.Full size imageIn Fig. 4 left column, the full driver set of one inactive stable motif or motif group was stabilized. However, the species that become extinct might only contain a subset of a driver set of a stable motif or motif group, i.e., they only stabilize a subset of the inactive node states in the stable motif or motif group. We compare the extinction effect caused by the stabilization of a full driver set of four nodes with the effect of the extinction of four nodes that contain a partial driver set in Fig. 5 using the batch of the largest networks in this study, i.e, the batch that contains networks with 30 nodes representing plant species and 40 nodes representing pollinator species. This choice is due to the fact that the existence of stable motifs and motif groups having a driver set of four node states is highly probable in larger networks. As expected, the stabilization of the complete driver set leads to greater damage. Stabilization of the same number of nodes that contain a partial driver set leads to significantly less damage and species loss in the community; the median damage percentage in the case of stabilization of partial driver sets is 22.6% while it is 69.2% in the case of stabilization of the full driver sets. We also note that damage of more than 90% occurs rarely and is only possible when a full driver set is stabilized (see Fig. 5 right plot). This suggests that the motif groups that lead to total extinction tend to have a driver set with more than four nodes; in other words, only the simultaneous extinction of five or more species would lead to total community collapse.Figure 5Histogram plots illustrating the average percentage of the damage caused in an established community after the extinction of species. This analysis is performed over 1000 networks with the size of 70 nodes (30 nodes representing plant species and 40 nodes representing pollinator species). The original extinctions are excluded from the damage percentages. The left plot shows the average damage percentage caused to the maximal richness community by the extinction of 2 species that are a subset of the 4-node driver set of an inactive stable motif or motif group plus 2 randomly selected non-driver species. The right plot shows the damage percentage caused to the maximal richness community by the extinction of 4-node driver sets of the same inactive stable motifs and motif groups. Each time the driver set of one stable motif or motif group was the subject of extinction and we calculated the average damage percentage over all data points collected for each network. 295 networks qualified for this analysis.Full size imageMotif driver set analysis outperforms structural measures in identifying keystone speciesThe literature on ecological networks offers multiple measures that reflect the importance of each species for community stability. One family of such measures is centrality (quantified by the network measures degree centrality and betweenness centrality). Previous studies45,46 have shown that species (nodes) with higher centrality scores are keystone species in ecological communities (i.e., species whose loss would dramatically change or even destroy the community). The nodes with highest in-degree centrality (such as pl_2 in Fig. 6a) represent generalist species that can receive beneficial interactions from multiple sources and survive. The nodes with highest betweenness centrality (such as pl_2 and po_2 in Fig. 6a) represent species that act as connectors and help the community survive. We find that high centrality corresponds to specific patterns in the expanded network: the inactive state of generalist or connector species is often the driver of a cascading extinction. Indeed, stable motif analysis of the expanded network in Fig. 6b confirms that there is an inactive stable motif (highlighted with grey) driven by the minimal set {(sim )pl_2}. The fact that node pl_2 is a stable motif driver means that in the case of the extinction of pl_2 the whole community collapses.To compare the effectiveness of stable motif analysis to the effectiveness of the more studied structural measures to identify keystone species, we performed an analysis similar to the previous section. We compared the magnitude of cascading extinctions in the case of extinction of stable motif driver nodes and of nodes with high values of previously introduced structural importance measures. Specifically, we used node betweenness centrality, node contribution to nestedness47, and mutualistic species rank (MusRank)22 to find crucial species based on their structural properties. For more details on definition and adaptation of these two measures see “Methods”. In this analysis, we used each measure to target species in the simplified Boolean models as follows:

    1.

    Betweenness centrality: The 10% of species with the highest betweenness centrality are chosen to be candidates for extinction.

    2.

    Node contribution to nestedness: The species with the most interactions tend to contribute the least to the community nestedness. Targeting them most likely leads to a faster community collapse48. As a result, 10% of species with the lowest contribution to network nestedness are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

    3.

    Pollinator MusRank: The pollinator species with the highest MusRank importance are more likely to interact with multiple plants, so the 10% of pollinator species with the highest importance are chosen to be candidates for extinction. For more details on this measure, please see “Methods”.

    4.

    Plant MusRank: The plant species with the highest MusRank importance are more likely to interact with multiple pollinators, so the 10% of plant species with the highest importance are chosen to be candidates for extinction.

    Figure 7 illustrates the results of this analysis in 6000 networks with 50–70 nodes. In each network the 1-node, 2-node, and 3-node driver sets of inactive stable motifs are identified and made extinct. In the same networks 10% of nodes based on betweenness centrality, node contribution to nestedness, and node MusRank score were chosen to be candidates for extinction. To match the “driver set” data, all choices of 1, 2, or 3 nodes in these sets were explored and the damage was averaged over each extinction size for each network. We observe the cascading extinction and calculate the damage percentage relative to the maximal richness attractor. The plot represents the collective data over all initial simultaneous extinction sizes of 1, 2, and 3 species.Comparing the four methods, one notices that the histograms acquired using stable motif driver sets, node betweenness centrality, and node contribution to nestedness are very similar, showing a peak for the 10–20% bin of the damage, and a long tail that reaches a damage percentage of 80–100%. The MusRank score performs less well in identifying the crucial species. Also, the frequency of the higher damage percentages shows that node contribution to nestedness is the closest to the “driver set” method in identifying nodes whose extinction causes the collapse of the whole community, making it the best structural measure out of the three. Nevertheless, the driver set method finds keystone species that cannot be identified via structural measures, as the corresponding damage percentage histogram has the most prominent tail at the right edge of the panel. Indeed, stable motif driver sets identified 82%, 80%, and 546% more species whose extinction leads to 60% or higher damage to the community when compared to betweenness centrality, node nestedness, and node MusRank score based methods respectively.The reason for the higher effectiveness of driver set analysis is illustrated in Fig. 8 in which the MusRank score and node contribution to nestedness are calculated for two example networks. One can see how these two measures might incorrectly identify less vital species. In the left column of Fig. 8, MusRank identifies the node po_2, highlighted with green, as the most important species. However, this node does not have any outgoing edges; its extinction does not lead to any cascading extinction. The inability of the MusRank score to consider the direction of edges causes such misidentification. In the right column, the three nodes highlighted with yellow have the lowest contributions to network nestedness. The expanded network shows that these three nodes together are not able to cause full community collapse, while the three-node driver set of the inactive stable motif can. Since the nestedness definition depends on the number of mutual interactions, it might fail to identify some of the keystone nodes that are necessary to the stability of the community (for more details on node nestedness see “Methods”).Previously it was shown that identifying the stable motifs and their driver sets can successfully steer the system toward a desired attractor or away from unwanted ones37,38,43. Stable motif analysis of the Boolean model offers insight into the dynamical trajectories of the system; hence control strategies can be developed accordingly. In the next section we use stable motif driver sets to suggest control methods and analyze their efficiency.Figure 6Generalist species in the interaction network and the expanded network. (a) A simplified network consisting of 3 plant and 3 pollinator species. pl_2 is a generalist species, i.e., it has two incoming edges indicating that it can survive on either of its sources of pollination, po_1 or po_2. The expanded network in (b) illustrates that the stabilization of the grey stable motif stabilizes all the nodes to their inactive states, and hence causes full community collapse. (sim )pl_2 is the minimal driver set of the grey stable motif, consistent with the strong damage induced by the loss of a generalist species.Full size imageFigure 7Histogram plots illustrating the performance of driver set analysis versus structural measures in identifying keystone species. The analysis was done on 6000 networks with sizes of 50–70 nodes. The starting point is the maximal richness community, i.e., the attractor in which the most species establish. For each network 1, 2, and 3 node(s) were selected and simultaneously fixed to their inactive states. After the cascading damage the new attractor is compared to the maximal richness attractor to calculate the damage percentage. The structural measures—betweenness centrality, node nestedness contribution, and node MusRank score—were calculated for all nodes in each network; the top 10% according to the relevant ordering were candidates to being fixed to their inactive states. The network IDs were matched, i.e., only the networks that had candidate nodes according to all four measures for each extinction size are included in this plot. The total number of data points is 6360. The red solid lines represent the mean and the black dashed lines represent the median over all data points in each plot.Full size imageFigure 8Networks illustrating examples of when structural measures fail to identify keystone species. In both columns simplified networks consisting of 3 plant and 3 pollinator species are presented. The MusRank is calculated for all the nodes in the network in the left column and denoted in the node labels. The expanded network corresponding to this network is shown below. Node contribution to network nestedness is calculated for all the nodes in the network in the right column and denoted in the node labels. Similarly the expanded network that correspond to it is shown below. Note that these two networks have different edges. In the left column MusRank score identifies node po_2, highlighted with green, as the most important, while the expanded network shows that the extinction of po_2 does not cause any further damage to the community, as this node has no outgoing edges. This is due to the fact that MusRank calculation process fails to consider the directed network and replaces all the directed edges with undirected ones. The MusRank score does not identify po_3 as a crucial species; however, virtual node (sim )po_3, outlined with black in the expanded network is a driver of a stable motif that has all other nodes in its LDOI; the extinction of po_3 leads to full community collapse. In the right column, the nodes highlighted with yellow (pl_2, pl_3, and po_2) have the lowest node contribution to nestedness, which predicts that these nodes are likely crucial to the stability of the community. Analyzing the expanded network, one can see that these three nodes together are not able to drive the inactive stable motif highlighted with teal. The minimal driver set for this stable motif, outlined with black, consists of {(sim )po_1, (sim )po_2, (sim )po_3}; together these nodes drive the inactive stable motif and cause full community collapse. The nestedness-based measure was not able to capture the significance of nodes po2 and po_3.Full size imageDamage mitigation measures and strategies for endangered communitiesThere are two substantial questions related to managing the damage induced by species extinction: (1) How can one prevent the damage as much as possible? (2) Once the damage happens, the reintroduction of which species can restore the community and to what extent? In this section we aim to answer these questions in the context of the Campbell et al. model, implementing stable motif based network control. This analysis can inform agricultural and ecological strategies employed to prevent and mitigate damage.Damage preventionOne of the most important questions in ecology is what strategies to use so that we can prevent and avert extinction damage to the community. In this section we analyze how the knowledge from stable motif analysis and driver sets can be implemented to minimize the effect of extinction of keystone species in case of limited resources. Each attractor of the original system can have multiple control sets; stabilizing the node states in each control set ensures that the system reaches that specific attractor. The same information from the attractor control sets can be implemented to prevent the system from converging into unwanted attractors. Zañudo et al. illustrated that by blocking (not allowing to stabilize) the stable motifs that lead to the unwanted attractors, one can decrease the probability (sometimes to zero) that the system arrives in those attractors38. In order to block an attractor, the control sets of that attractor are identified and the negations of the node states in the control sets are externally imposed. This approach eliminates the undesired attractor; however, new attractors might form that are similar to the eliminated attractor. Campbell et al. showed that in order to avoid such new attractors one needs to block the parent motif, which in this case is the largest strongly connected subgraph of the expanded network that contains the inactive virtual nodes44. Here, we investigate how stable motif blocking based attractor control can identify the species whose preservation would offer the highest benefit in avoiding catastrophic damage to the community. This information would aid the development of management strategies in plant–pollinator communities.To avoid all attractors that lead to some degree of species extinction, one needs to block all the driver sets of all inactive stable motifs and motif groups in a given network. Implementing this in 100 randomly selected networks with 25 plant and 25 pollinator nodes, we found that 45.6% of the species in the maximal richness community need to be kept (prevented from extinction) to ensure the lack of cascading extinctions. Given that management resources are usually limited, active monitoring and conservation of almost half of the species in a community seems costly and impractical. Hence, we set a more feasible goal of identifying and blocking the driver set(s) of the largest inactive stable motif or motif group in each network. The same 100 networks containing 50 nodes are the subject of analysis in this section. The reason for performing the analysis in a relatively limited ensemble is that it involves the identification of all driver sets of the largest inactive stable motif or motif group, which is computationally expensive. For each network, the driver set of the largest inactive stable motif or motif group (which corresponds to the extinction of all the species in that group) is identified and blocked (that is, the corresponding species are not allowed to go extinct). Then the same number of species as in the driver set of that stable motif or motif group are selected and stabilized to their inactive state. We considered all combinations of node extinctions outside the blocked subset, calculated the damage percentage relative to the maximal richness community, and then averaged over all data points for each network. As a control, we repeated the analysis without blocking; the size of the initial extinction is the same as in the previous analysis for consistency.Figure 9 shows the result of the analysis described above for 100 networks. The left box and whiskers plot illustrates the damage percentage relative to the maximal richness community when the blocking feature is activated, while the right box and whiskers plot shows the damage percentage relative to the maximal richness community when the blocking is disabled. The average and median damage percentages are 14.96% and 13.04% respectively when the largest inactive stable motif or motif group was blocked and 24.73% and 20.38% when it was not. This (sim )10% difference in the average between the two sets of results, as well as the fewer cases of high-damage outliers in the left plot, demonstrates that by preventing the extinction of species identified by stable motif analysis, one can prevent catastrophic community damage considerably.To estimate the fraction of species that would need to be monitored to prevent their extinction, we compared the size of the maximal richness attractor and the size of the driver set of the largest stable motif. The maximal richness community represents an average of 32% of the original species pool, approximately 15 out of 50 species. The driver sets of the largest stable motifs had an average size of 2.5 node states over all 100 networks, i.e., about 16.6% of the maximal richness community. In ecological terms, given limited resources, the information gained from stable motif driver sets can help direct the conservation efforts toward the keystone species that play a key role in maintaining the rest of the community in a cost-effective manner.Figure 9Box plots comparing the damage communities face if the largest inactive stable motif or motif group is completely blocked, i.e., all the drivers of this inactive stable motif or motif group are prevented from stabilizing versus if the same stable motif or motif group is allowed to stabilize. This analysis was performed over 100 randomly selected networks that contain 25 plant and 25 pollinator nodes. All the driver sets of an inactive stable motif or motif group are identified. From left to right the box and whiskers plots show the average damage percentage relative to the maximal richness community if the largest inactive stable motif is blocked and the same quantity if the largest stable motif or motif group is not blocked respectively. For the left box and whiskers plot, all combinations of inactive node states except the driver sets are considered, and for the right box and whiskers plot all combinations are explored. Due to the computational complexity caused by combinatorial explosion, this analysis was performed over 100 randomly selected 50-node networks.Full size imageRestoration of a group of speciesAlthough human preservation efforts have been directed toward community conservation, there are many industrial activities that lead to ecosystem degradation. Ecologists are interested in developing restoration strategies to be deployed after a stable community is hit by catastrophic damage to recover biodiversity and the ecosystem functions it provides49. Here we propose that stable motif analysis and the driver sets identified from the expanded network can give insight into restoration measures. While we examined the inactive stable motifs in the study of species extinction, here we focus on the active stable motifs as our goal is to restore as much biodiversity as possible.Several network measures have been proposed to identify the species that if re-introduced would restore the community considerably. Two of the most studied algorithms include maximising functional complementarity (or diversity) and maximising functional redundancy50. The first strategy targets the restoration of the species that provide as many functions to the ecosystem as possible; this approach results in a community that has a maximal number of functions provided by different groups of species. Alternatively, maximising the functional redundancy yields a community in which several species perform the same function. While this resultant community might have a limited number of functions, it is robust. Both of these community restoration approaches have been studied extensively (e.g. see21).We hypothesize that restoring the species that constitute driver sets of active stable motifs can help maximise the number of species post-restoration. Since there is evidence that functional diversity correlates with the number of species in the community51, we compare the post-restoration communities identified by stable motif driving with the functional diversity maximisation approach. As discussed in section LDOI in the Boolean threshold model, the Boolean simplification of the threshold functions leads to an overestimation of the LDOI of active node states (compared to the original threshold functions) in some networks. We evaluate the negative effects of this overestimation by checking the effectiveness of the restored species in the original threshold model.The same 6000 networks we examined in the last section were the subject of this analysis. To create an unbiased initial community, we create the damaged communities by eliminating the same number of species from the maximal richness community as the number that will be restored. We identify the inactive stable motif or motif group with the driver set size of 1, 2, or 3 node states that causes the most damage to the maximal richness community. We then eliminate the species corresponding to this driver set to reach the most damaged community for the given size of the initial extinction. This community is the starting point for two analyses. In the stable motif driving approach we stabilized an active stable motif that has a driver set of the same size as the initial extinction to reach a post-restoration community and calculated the percentage of the extinct species that were restored. In the functional diversity maximization based approach we re-introduced the same number of species selected from the to 10% of species in terms of their contribution to functional diversity.To calculate the functional diversity of a community one needs to (1) define and construct a trait matrix, (2) determine the distance (trait dissimilarity) of pairs of species, (3) perform hierarchical clustering based on the distances to create a dendrogram, and (4) calculate the total branch length of the dendrogram, i.e., the sum of the length of all paths51,52. Petchey et al. argued that resource-use traits among plant and pollinator species can be used to classify the organisms into separate functional groups53 and Devoto et al. proposed the use of the adjacency matrix based on the interaction network as the trait matrix21. In this study we do the same and implement the bipartite adjacency matrix to construct the distance matrix.Since the networks of the Campbell et al. model are directed, we modify the algorithm in that we have two separate adjacency matrices, one denoting the edges incoming to plant species and the other denoting the edges incoming to pollinator species. The hierarchical clustering algorithm is then run on each of these matrices separately, resulting in a dendrogram for each adjacency matrix. If extinction occurs in a community, the functional diversity of the survived community can be determined by calculating the total branch length of the subset of the dendrogram that includes only the survived species. The restoration strategy using this method is to re-introduce the nodes whose branches add the most to the total branch length of this subset, i.e., maximise the functional diversity of the survived community54. For more details see “Methods”.In each network, the percentage of the extinct species that were restored was calculated and averaged over all data points for each restoration size and each network. Figure 10 illustrates the results of this investigation. Applied to the simplified Boolean model, the median restoration percentage in the case of active stable motif driver set method (blue plot) is 80%. The functional diversity maximization strategy to restoration (yellow plot) yields a lower median restoration percentage, 73%, as well as a large number of low-restoration outliers. Although one might argue that identifying beneficial species using the functional diversity maximization strategy works well, the higher percentage of the cases of 80–100% restoration in case of the active stable motif driver set analysis indicates that the latter identifies some of the most effective restorative species that are not identified via the former method. As in a minority of cases the simplified Boolean model overestimates the positive impact of the sustained presence of a species (see section LDOI in the Boolean threshold model), we sought to verify the effectiveness of the predicted restoration candidates in the original threshold model. The blue (respectively, yellow) box and whiskers plot on the right represents the restoration percentages of the same species as in the left blue (respectively, yellow) plot when these species are restored in the threshold model. The median of the right blue plot is 70%, while the median of the right yellow is 63%, preserving the advantage of the stable motif driver sets. We conclude that although the simplified Boolean model overestimates the restoration effectiveness of certain driver sets (visible in the fact that the lower whisker of the blue plot on the right goes well below the lower whisker of the blue plot on the left), stable motif driver sets are more effective in both comparisons.Figure 10Box and whiskers plots illustrating the average percentage of the extinct species that are restored following the stable motif driver set restoration strategy (blue) versus the functional diversity based approach (yellow). This analysis is performed over 6000 networks with sizes of 50–70 nodes. Starting from the maximal richness community, for each network one inactive stable motif with a driver set of 1, 2 or 3 nodes was stabilized to reach a new damaged community. This task was performed until the community with the most extinct species was identified. This is the community we set as the starting point for the restoration process using both methods. The pair on the left represents the two methods applied to the simplified Boolean model. For both methods we identified 1, 2, or 3 influential nodes for community restoration and we calculated the percentage of the extinct species that could be restored. The pair on the right represents restoring the same species identified by each method in the previous analysis in the original threshold model. In all analyses the community restoration percentage was averaged over all combinations of the same size, for each network and each method. The IDs of all networks are matched.Full size imageCommunity restoration via attractor controlAs illustrated in section “Restoration of a group of species”, stable motif analysis identifies promising and cost-effective group restoration strategies. In this section we aim to go further and identify interventions that can maximally restore a community. Previous stable motif based network control methods37,38,55 require a search for the smallest set of node states to control the system once the stable motif stabilization trajectories are identified. This smallest set may not contain a node from each stable motif in the sequence. In this work, however, we know that each stable motif or motif group needs to be controlled individually28 because the stabilization of none of the motifs results in the stabilization of another. As a result, the control set of each attractor is the same as the union of the driver sets of all members in the consistent combination corresponding to that attractor.In this section we examined this attractor control method by setting the communities with 70% or more of the species in the maximal richness community as the target, i.e., the attractors that have 70% of the species in the maximal richness community are assumed to be the desired attractors. We then recorded the size of the minimal control set needed to achieve each of these attractors. Note that stabilizing each of these control sets guarantees that the system reaches the corresponding attractor38.For this section, we analyzed 6000 networks that have 50–70 nodes. Figure 11 represents box-and-whiskers plots of the size of the minimal set of species that need to be restored, where the target community sizes are classified into three groups based on the percentage of the species relative to the maximal richness attractor. One can see that in half of the cases, the restoration of either 1 or 2 species manages to restore more than 70% of the maximal richness community. The largest set has 8 species that need to be restored; however, this data point is an outlier. As illustrated, driver set analysis and stable motif based attractor control can efficiently identify the species that play an influential restorative role and suggest management strategies that are effective at the scale of the whole community. To assess the impact of the LDOI inflation on this result, we used the restoration candidates identified by control sets of the attractors of the Boolean model in the threshold functions of a subset of networks. The results of comparing the restoration percentage is shown in Fig. 14. The first quartile, median and third quartile values are 78.26%, 86.6%, and 100% for the simplified Boolean models and 43.78%, 72.41%, and 85.71% for the threshold model.To further compare the results of restoration obtained from the two models we sorted the species in the order of their contribution to community restoration following a catastrophic damage. We randomly selected 100 of the largest (70-node) networks, which have the highest probability of a discrepancy between the threshold functions and the simplified Boolean model. In 72% of the cases the two rankings matched completely, and in the majority of the remaining cases only one species was misplaced in the simplified Boolean model-based ranking. To conclude, there is a significant advantage to the implementation of the simplified Boolean model and the drawback can be addressed by a follow-up checking on the original threshold functions.Figure 11The number of species that need to be restored to save 70% of more of the species in the maximal richness community. In this analysis 6000 networks with 50–70 nodes were the subject. For each networks all the attractors that have 70% or more of the species in the maximal richness attractor are identified and set to be the target attractors. The control set of these attractors are then classified into three groups based on the percentage as illustrated in the figure. From left to right, the box and whiskers represent the size of the control set of attractors that have 70–80%, 80–90%, and 90–100% of the species in the maximal richness attractor respectively.Full size image More

  • in

    Plant nitrogen retention in alpine grasslands of the Tibetan Plateau under multi-level nitrogen addition

    Study siteThe field experiment was conducted at Namco Station (30°47’N, 90°58’E, altitude 4730 m) of the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS), which is located in the alpine steppes of TP in China. The experiment was permitted by ITPCAS, complied with local and national guidelines and regulations. From 2006 to 2017, the mean annual temperature (MAT) and mean annual precipitation (MAP) was about − 0.6 °C and 406 mm, respectively. Monthly mean temperature varied from − 10.8 °C in January to 9.1 °C in July and most of the precipitation occurred from May to October37,38. During our six-year observations (2010, 2011, 2012, 2013, 2015 and 2017), climate change during the growing season from May to September varied differently, with the annual precipitation ranged from 255.9 mm to 493.8 mm and the MAT from 6.7 to 7.4 °C. Androsace tapete, Kobresia pygmaea, Stipa purpurea and Leontopodium pusillum were the dominant plant species at the alpine steppe.Experimental design and treatmentsThe long-term experiment began in May, 2010. Three homogenous plots were randomly arranged as replicates at the alpine steppe and six subplots (~ 13 m2) were distributed in each plot by a cycle, with a 2 m buffer zone between each adjacent subplot (Appendix S1: Fig. S1). In this experiment, six treatments of N fertilization rate (0, 1, 2, 4, 8, and 16 g N m−2 yr−1) were clockwise applied in each subplot. The subplots of 0 g N m−2 yr−1 were control group. We sprayed NH4NO3 solution on the first day of each month in the growing season (from May to September) each year. After fertilizing, we rinsed plant residual fertilizer with a little deionized water (no more than 2 mm rainfall). For the control groups, we added equivalent amount of water. The experiment was conducted from 2010 to 2017 (it should be pointed out that there was no fertilization in 2014 and 2016).Sampling and measurementsThe samples were collected with the training and permission of ITPCAS and involved plants that are common species and not endangered or protected. The identification of the plants was done by referring to a book of Chen and Yang39. Pictures of the corresponding specimens can be seen on the website of ITPCAS (http://itpcas.cas.cn/kxcb/kxtp/nmc_normal_plant/).Vegetation samples were collected in August in 2011 and repeated at the same time in 2012, 2013, 2015 and 2017. We established one 50 × 50 cm quadrat in each subplot, clipped aboveground biomass (AGB) and sorted species by families. The biomass was used to measure ANPP (g m−2 yr−1). Following aboveground portion collected, we used three soil cores (5 cm diameter) to collect the belowground roots at 0–30 cm depth and mixed into one sample, which were used to assess belowground net primary productivity (BNPP, g m−2 yr−1). The roots were cleaned with running water to remove sand and stones.Both plant and root samples were dried at 75 °C for 48 h and then ground into powder (particle size ~ 5 μm) by a laboratory mixer mill (MM400, Retsch). To determine N and C content of plants, we weighed the samples into tin capsules and measured with the elemental analyzer (MAT253, Finnigan MAT GmbH, Germany).Estimation of the critical N rate (Ncr), N retention fraction (NRF), N retention capacity and N-induced C gainAccording to the N saturation hypothesis, plant productivity increases gradually during N addition, reaches a maximum at the Ncr, and eventually declines16,17. We considered the Ncr to be the rate where ANPP no longer remarkably changed with N addition (Fig. 1).We defined plant N retention fraction (NRF, %; Eq. 1) as the aboveground N storage caused by unit N addition rate, and N retention capacity (g N m−2 yr−1; Eq. 2) was the increment of N storage due to exogenous N addition compared to the control40. The equations are as following:$$N;retention;fraction = frac{{ANPP_{tr} times N;content_{tr} – ANPP_{ck} times N;content_{ck} }}{N;rate}$$
    (1)
    $$N;retention;capacity = ANPP_{tr} times N;content_{tr} – ANPP_{ck} times N;content_{ck}$$
    (2)
    where ANPPtr and N contenttr (%) refer to those in the treatment (tr) groups, and ANPPck and N contentck refer to those in the control (ck) groups. These expressions are also used in the following equations (Eqs. 3–5).The N-induced C gain (g C m−2 yr−1; Eq. 3) was estimated by the increment of C storage owing to exogenous N addition compared to the control40. Maximum N retention capacity (MNRC, Eq. 4) and maximum N-induced C gain (Eq. 5) mean the maximum N and C storage increment in plant caused by exogenous N input at Ncr, respectively. The formulas are as following:$$N{text{-}}induced;C;gain = ANPP_{tr} times C;content_{tr} – ANPP_{ck} times C;content_{ck}$$
    (3)
    $$MNRC = ANPP_{max } times N;content_{max } – ANPP_{ck} times N;content_{ck}$$
    (4)
    $$Maximum;N{text{-}}induced;C;gain = ANPP_{max } times C;content_{max } – ANPP_{ck} times C;content_{ck}$$
    (5)
    where ANPPmax, N contentmax and C contentmax refer to the value of ANPP, N content and C content at Ncr, respectively.Data synthesisTo evaluate N limitation and saturation on the TP more accurately, we searched papers from the Web of Science (https://www.webofscience.com) and the China National Knowledge Infrastructure (https://www.cnki.net). The keywords used by article searching were: (a) N addition, N deposition or N fertilization, (b) grassland, steppe or meadow. Article selection was based on the following conditions. First, the experimental site must be conducted in a grassland ecosystem. Second, the experiment had at least three N addition levels and a control group. Third, if the experiment lasted for many years, we analyzed data with multi-year average. Based on the above, we collected 89 independent experimental cases. Among these, 27 cases were located on the TP alpine grasslands, 25 in the Inner Mongolia (IM) grasslands and 37 in other terrestrial grasslands (detailed information sees Appendix S2: Table S1).We extracted ANPP data and N addition rate of these cases and estimated Ncr and ANPPmax (Appendix S2: Fig. S2). We then calculated NRF, N retention and C gain of each group of data for further analysis (Appendix S2: Table S2). Most of the 89 cases did not have data on N and C content. To facilitate the calculation, we summarized N and C content from 40 articles in the neighboring areas of the cases and divided the N and C content into seven intervals according to the N addition rate (Appendix S2: Table S3 and Fig. S3). The unit of N addition rate was unified to “g N m−2 yr−1”. All the original data were obtained directly from texts and tables of published papers. If the data were displayed only in graphs, Getdata 2.20 was used to digitize the numerical data. For the estimation of N retention and C gain of the TP at current N deposition rates and future at Ncr, we fitted the exponential relationship to the data from 27 cases on the TP, and then substituted N rates into the fitted equations (Eq. 6):$$y = a times left[ {1 – exp left( { – bx} right)} right].$$
    (6)
    We also included MAT, MAP, soil C:N ratio, fencing management (fencing or grazing) and grassland type (meadow, steppe and desert steppe) of the experiment sites for exploring the drivers affecting N limitation (Appendix S2: Table S1). When climatic data were missing from the article, MAT and MAP were obtained from the WorldClim (http://www.worldclim.org).Species were usually divided into four functional groups (grasses, sedges, legumes and forbs) to study the response of species composition to N addition in previous study41. We synthesized 13 TP experimental cases (including our field experiment) from the data synthesis and each case included at least three functional groups (detailed references see Appendix S2).Statistical analysisThere were 42 species in our field experiment. We divided them by family into eleven groups: Asteraceae (forbs), Poaceae (grasses), Leguminosae (legumes), Rosaceae (forbs), Boraginaceae (forbs), Caryophyllaceae (forbs), Cyperaceae (sedges), Labiatae (forbs), Primulaceae (forbs), Scrophulariaceae (forbs) and Others. Due to species in the group of Others contributed only 1.22% of AGB, we analyzed AGB and foliar stoichiometry among other ten families (Appendix S1: Table S1). In Namco steppe, forbs, grasses, sedges and legumes accounted for 78.0%, 7.4%, 8.2% and 5.2% of the AGB respectively (Appendix S1: Table S1 and Fig. S2). Such a large number of forbs suggested that our experiment was conducted on a severely degraded grassland.For our field data, two-way ANOVAs were used to analyze the effects of year, N fertilization rate and their interactions on species AGB. One-way ANOVAs were used to test the response of ANPP, BNPP, root:shoot ratio, species foliar C content, N content and C:N ratio to N addition rate. Duncan’s new multiple range test was used to compare the fertilization influences at each rate in these ANOVAs. Prior to the above ANOVAs, we performed homogeneity of variance test and transformed the data logarithmically when necessary. Simple regression was used to estimate the relevance among ANPP, NRF, N retention capacity and C gain with N addition rates.Structural equation modeling (SEM) was used to explore complex relationships among multiple variables. To quantify the contribution of drivers such as climate and soil to Ncr, ANPP, NRF and MNRC, we constructed SEM based on existing ecological knowledge and the possible relationships between variables. We considered environmental factors (MAT, MAP and soil C:N) and ANPPck as explanatory variables, and Ncr, NRF and MNRC as response variables. We included the ANPPck in the SEM rather than the ANPPmax because we wonder whether there was a relationship between ANPP in the absence of exogenous N input and the ecosystem N retention in the presence of N saturation. This has important implications for assessing N input. Before constructing the SEM, we excluded collinearity between the factors. In addition, Student’s t-test and one-way ANOVAs were performed to explain the effect of fencing management and grassland type on above response variables, respectively. The SEM was constructed using the R package “piecewiseSEM”42. Fisher’s C was used to assess the goodness-of-model fit, and AIC was for model comparison.Given the influence of extreme values in the data synthesis, we calculated the geometric mean of Ncr, NRF, N retention and N-induced C gain. All statistical analyses were performed with SPSS 26.0 and RStudio (Version 1.2.1335) based on R version 3.6.2 (R Core Team, 2019). More