# MARINE MICROBIOME AND BIOGEOCHEMICAL CYCLES IN MARINE PRODUCTIVE AREAS

EDITED BY : Alejandro A. Murillo, Veronica Molina, Chris Harrod and Julio Salcedo-Castro PUBLISHED IN : Frontiers in Marine Science and Frontiers in Microbiology

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-276-3 DOI 10.3389/978-2-88963-276-3

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# MARINE MICROBIOME AND BIOGEOCHEMICAL CYCLES IN MARINE PRODUCTIVE AREAS

Topic Editors: Alejandro A. Murillo, European Molecular Biology Laboratory Heidelberg, Germany Veronica Molina, Universidad de Playa Ancha, Chile Chris Harrod, University of Antofagasta, Chile Julio Salcedo-Castro, Universidad de Playa Ancha, Chile

Citation: Murillo, A. A., Molina, V., Harrod, C., Salcedo-Castro, J., eds. (2019). Marine Microbiome and Biogeochemical Cycles in Marine Productive Areas. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-276-3

# Table of Contents

*05 Editorial: Marine Microbiome and Biogeochemical Cycles in Marine Productive Areas*

Alejandro A. Murillo, Verónica Molina, Julio Salcedo-Castro and Chris Harrod

*08 Temporal Variability in Net Primary Production in an Upwelling Area off Central Chile (36°S)*

Giovanni Testa, Italo Masotti and Laura Farías

*25 First Record of Microbiomes of Sponges Collected From the Persian Gulf, Using Tag Pyrosequencing*

Akram Najafi, Maryam Moradinasab and Iraj Nabipour

Oliver Müller, Lena Seuthe, Gunnar Bratbak and Maria L. Paulsen

*69 Diverse Bacteria Utilize Alginate Within the Microbiome of the Giant Kelp*  Macrocystis pyrifera

Jordan D. Lin, Matthew A. Lemay and Laura W. Parfrey

Italo Masotti, Pilar Aparicio-Rizzo, Mariela A. Yevenes, René Garreaud, Lucy Belmar and Laura Farías

*123 Controls on the Fate of Dissolved Organic Carbon Under Contrasting Upwelling Conditions*

Mariana B. Bif, Dennis A. Hansell and Kimberly J. Popendorf

*134 Picoeukaryotic Diversity and Activity in the Northwestern Pacific Ocean Based on rDNA and rRNA High-Throughput Sequencing* Feipeng Wang, Yuyuan Xie, Wenxue Wu, Ping Sun, Lei Wang and Bangqin Huang

*147 Microbial Communities in Coastal Glaciers and Tidewater Tongues of Svalbard Archipelago, Norway*

Eva Garcia-Lopez, Irene Rodriguez-Lorente, Paula Alcazar and Cristina Cid

*160 Microbial Community Structure and Functionality in the Deep Sea Floor: Evaluating the Causes of Spatial Heterogeneity in a Submarine Canyon System (NW Mediterranean, Spain)* Sara Román, Rüdiger Ortiz-Álvarez, Chiara Romano, Emilio O. Casamayor

and Daniel Martin

# Editorial: Marine Microbiome and Biogeochemical Cycles in Marine Productive Areas

#### Alejandro A. Murillo1 \*, Verónica Molina2 , Julio Salcedo-Castro2 and Chris Harrod3,4

1 Structural and Computational Biology, European Molecular Biology Laboratory, Heidelberg, Germany, 2 Departamento de Biología, Universidad de Playa Ancha, Valparaíso, Chile, 3 Instituto de Ciencias Naturales Alexander von Humboldt, Universidad de Antofagasta, Antofagasta, Chile, 4 Núcleo Milenio INVASAL, Concepción, Chile

Keywords: coastal areas, upwelling systems, carbon cycle, nitrogen cycle, microbial oceanography, high nutrients, microbial diversity, coastal ecosystems

**Editorial on the Research Topic**

#### **Marine Microbiome and Biogeochemical Cycles in Marine Productive Areas**

The global ocean covers more than 70% of the Earth's surface and is remarkably heterogeneous. Marine productive areas, and coastal ecosystems comprise a minor fraction of the ocean in terms of surface area, yet have an enormous impact on global biogeochemical cycles carried out by microbial communities, which represent 90% of the ocean's biomass (Alexander et al., 2011). These marine productive areas also vary considerably in terms of primary and secondary production, and are classified as highly productive or class I (>300 gC m−2 yr−1 ), moderately productive or class II (150–300 gC m−2 yr−1 ), to low productive or class III (<150 gC m−2 yr−1 ). Work to date has largely focused on cycling of carbon and macronutrients such as nitrogen, phosphorus, and silicate: other important elements such as sulfur or trace elements have been less studied, reflecting associated technical and logistical issues. Increasingly, these marine areas, and the taxa that form their ecosystems, are subject to significant anthropogenic pressure, impacting marine life and recycling of energy and nutrients (Galton, 1884; Hasler, 1969; Jickells et al., 2017). A key example is that of cultural eutrophication, where agricultural run-off leads to nitrogen and phosphorus enrichment of coastal ecosystems, greatly increasing productivity resulting in algal blooms, deoxygenation of the water column and seabed, and increased greenhouse gas emissions (Bouwman et al., 2005), with direct local and global impacts on nitrogen and carbon cycles. However, the runoff of organic matter from the mainland to coastal ecosystems is just one of a series of pressing threats stressing microbial communities due to global change. Climate change has also resulted in changes in the cryosphere, as glaciers and permafrost melt, resulting in intensified marine stratification, while shifts of the redox-state in different biomes are rapidly reshaping microbial assemblages at an unprecedented rate (Altieri and Gedan, 2015; Breitburg et al., 2018; Cavicchioli et al., 2019; Hutchins et al., 2019).

Global change is, therefore, affecting key processes including primary productivity, CO2 and N2 fixation, organic matter respiration/remineralization, and the sinking and burial deposition of fixed CO2 (Hutchins et al., 2019). In addition to this, oceans are experiencing an acidification process, with a change of ∼0.1 pH units between the pre-industrial period and today, affecting carbonate/bicarbonate buffer chemistry. In turn, acidification has been reported to impact planktonic communities, principally through effects on calcifying taxa (Stillman and Paganini, 2015). There is also evidence for shifts in the production of key intermediary volatile products, some of which have marked greenhouse effects (e.g., N2O and CH4, reviewed by Breitburg et al., 2018), due to the increase in global temperature, ocean stratification and deoxygenation, driving

#### Edited and reviewed by:

François G. Schmitt, UMR8187 Laboratoire d'Océanologie et de Géosciences (LOG), France

#### \*Correspondence:

Alejandro A. Murillo<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Coastal Ocean Processes, a section of the journal Frontiers in Marine Science

Received: 16 July 2019 Accepted: 09 October 2019 Published: 22 October 2019

#### Citation:

Murillo AA, Molina V, Salcedo-Castro J and Harrod C (2019) Editorial: Marine Microbiome and Biogeochemical Cycles in Marine Productive Areas. Front. Mar. Sci. 6:657. doi: 10.3389/fmars.2019.00657

as much as ∼25–50% of nitrogen loss from the ocean to the atmosphere in the so-called oxygen minimum zones (Bertagnolli and Stewart, 2018) or anoxic marine zones (Ulloa et al., 2012), driven by microbial processes. Other products, that are typically toxic for the marine nekton, including reduced sulfur species such as H2S, have a negative impact for marine resources like fisheries and coastal aquaculture. While global change has accelerated, there has been a parallel increase in awareness of the complexity of marine ecosystems, and especially the fundamental role of microbes as drivers of ecosystem functioning (Cavicchioli et al., 2019). At the same time, calls have intensified for the need for representative ecosystem-based models that can guide effective management of marine ecosystems. An example of this is the call made by the United Nations, proclaiming the Decade of Oceans Science for Sustainable Development 2012–2030 (https:// en.unesco.org/ocean-decade).

This Research Topic is focused on studies characterizing the complexity of interactions seen in microbial communities, and their response to environmental variation across a comprehensive range of marine habitats. These extend from coral reefs stressors and microbialization (Meirelles et al.), through to coastal habitats including glaciers in a Norwegian archipelago (Garcia-Lopez et al.), to the deep-sea floor and the spatial heterogeneity of a submarine canyon system (Román et al.). Other authors examined the picoeukaryotic community dynamics in the Northwestern Pacific Ocean (Wang et al.), the influence of river discharges (Masotti et al.) and temporal variability of net primary production (Testa et al.) both in an upwelling system of the Eastern Subtropical South Pacific Ocean. The role of microbial communities in the carbon cycle is at the core of our Research Topic, including the fate of dissolved organic carbon in upwelling systems (Bif et al.), the utilization of onecarbon and methylated compounds in the coastal environment (Dinasquet et al.) and dissolved organic matter derived from permafrost thaw in Arctic fjords, influencing the microbial communities (Müller et al.). This diversity of environments and microbial processes provides a small glimpse of the multiple roles that microbial communities play in marine ecosystems. Other important aspects include the potential for

### REFERENCES

biotechnological applications regarding the discovery of new bio-active compounds in a mangrove soil community (Chen et al.), the first sponge microbiome from the Persian Gulf (Najafi et al.) and the microbiome of the ecological engineer giant kelp Macrocystis pyrifera, including the utilization of alginate (Lin et al.).

 This includes the identification of novel bioactive compounds, development of green-blue energy, effective and efficient sewage treatments, and the mitigation of CO2 production, as carbon burial, and other greenhouse gases. By increasing our understanding of the functioning of those natural processes that allow our planet to auto-regulate and respond to perturbations, we can apply such knowledge to surmount current and imminent environmental problems. Life is thought to have originated in the ocean and since then, to the actual state of the biosphere, marine habitats have been the support to an astonishing diversity of life forms and metabolic pathways. As such, increasing our understanding on how the largest biome on Earth functions is possibly the most important enterprise of the upcoming decades and can provide answers and solutions for many of the current challenges that humankind faces.

### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### ACKNOWLEDGMENTS

We thank all the contributors to this Research Topic.

Ulloa, O., Canfield, D. E., DeLong, E. F., Letelier, R. M., and Stewart, F. J. (2012). Microbial oceanography of anoxic oxygen minimum zones. Proc. Natl. Acad. Sci. U.S.A. 109, 15996–16003. doi: 10.1073/pnas.1205009109

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Murillo, Molina, Salcedo-Castro and Harrod. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Temporal Variability in Net Primary Production in an Upwelling Area off Central Chile (36◦S)

Giovanni Testa1,2,3, Italo Masotti 3,4 and Laura Farías 2,3 \*

1 Graduate Program in Oceanography, Department of Oceanography, University of Concepción, Concepción, Chile, 2 Department of Oceanography, Faculty of Natural and Oceanographic Sciences, University of Concepción, Concepción, Chile, 3 Center for Climate and Resilience Research, University of Chile, Santiago de Chile, Chile, 4 Facultad de Ciencias del Mar y de Recursos Naturales, Universidad de Valparaíso, Viña del Mar, Chile

The temporal variability of Net Primary Production (NPP) off central Chile (36◦S, 73◦W), an area subjected to seasonal coastal upwelling, was analyzed using monthly in situ 13C incubations within the photic zone, along with bio-oceanographic variables from a fixed time series station; and satellite NPP estimations (NPPE) from the Vertically Generalized Production Model between 2006 and 2015. NPP and NPPE rates varied from 0.03 to 18.29 and from 0.45 to 9.07 g C m−2 d −1 , respectively. Both rates were fairly well correlated with each other (r 2 = 0.61), but when these data were separated into two periods, higher r 2 value was found during winter (r 2 = 0.70) with respect to the rest of the year (r 2 = 0.24); the latter correlation was partially due to increased weekly NPPE variability during active and relaxed upwelling events. NPP rates along with other biophysical variables allowed for a division of the annual cycle into three distinct periods: September to January (high productivity, mean integrated NPP rates of 4.0 g C m−2 d −1 ), February to March (intermediate productivity, mean integrated NPP rates of 1.4 g C m−2 d −1 ), and May to August (basal level, mean integrated NPP rates of 0.5 g C m−2 d −1 ). NPP appeared to be partially controlled by nutrient inputs, either from upwelling (September-April) and river discharge (May-August), maintaining high NPP rates throughout the entire year, with an annual mean NPP rate of 1.1 kg C m−2 yr−1 . In this region, El Niño Southern Oscillation events did not appear to impact the NPP interannual variability.

#### Edited by:

Alejandro A. Murillo, EMBL Heidelberg, Germany

#### Reviewed by:

X. Antón Álvarez-Salgado, Consejo Superior de Investigaciones Científicas (CSIC), Spain Laura Lorenzoni, University of South Florida, United States

> \*Correspondence: Laura Farías<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Coastal Ocean Processes, a section of the journal Frontiers in Marine Science

Received: 07 February 2018 Accepted: 07 May 2018 Published: 29 May 2018

#### Citation:

Testa G, Masotti I and Farías L (2018) Temporal Variability in Net Primary Production in an Upwelling Area off Central Chile (36◦S). Front. Mar. Sci. 5:179. doi: 10.3389/fmars.2018.00179 Keywords: primary production rates, coastal upwelling, river discharge, El Niño Southern oscillation, central Chile

### INTRODUCTION

The Eastern Boundary Upwelling Systems (EBUS) are considered the most productive regions in the world's oceans and sustain the greatest portion of global fishery production (Pauly and Christensen, 1995; Cubillos et al., 1998b). This high productivity is mainly the result of the coastal upwelling process by which subsurface, nutrient-rich waters ascend to the surface in response to wind stress on surface water, fertilizing waters in the photic zone and promoting phytoplankton blooms. Lachkar and Gruber (2012) identified equatorward alongshore wind, eddy activity, the width of the continental shelf and mixed layer depth as key physical and environmental factors controlling Net Primary Production (NPP) in the EBUS. On the other hand, Messié and Chavez (2015) detected wind-driven macronutrient supply, iron provision, temperature, light and offshore export as main drivers of timespace NPP variability in the EBUS.

The coastal upwelling systems off Chile extend over a wide latitudinal range (∼19–40◦ S), along irregular coastline with diverse topographic features (including canyons and the continental shelf) and variable river runoff. NPP may therefore be temporally and spatially variable. Coastal upwelling off central Chile has been relatively well documented with regards to the hydrological properties and dynamics of upwelling waters (Montecino et al., 2006; Escribano and Schneider, 2007; Kämpf and Chapman, 2016). Additionally, biological properties (such as Chlorophyll-a distribution) and their roles in determining productivity have been investigated (Montero et al., 2007; Iriarte et al., 2012). However, a significant gap remains for understanding temporal NPP dynamics and rates, which in turn limits the knowledge of how NPP rates and their variability are maintained. Information on marine primary production has been centrally important in understanding the fluxes of energy and organic matter in the ocean and the major biogeochemical cycles (e.g., C, N, P).

NPP measurements off the Chilean coast have been performed using different techniques, including O2, 14C and 13C methods (Montecino et al., 1998; Daneri et al., 2000; Farías et al., 2009; Iriarte et al., 2012). Volumetric primary production rates have been derived from different methods, each one with intrinsic assumptions and biases, which lead to differences in these measurements (Cullen, 2001; Regaudie-de-Gioux et al., 2014). The oxygen evolution (light and dark incubation) method was one of the first techniques capable of determining gross primary production (GPP) and has been widely applied due to its simplicity. However, the 14C method (Steeman Nielsen, 1952) rapidly became the standard for the oceanographic community, and has been used to calibrate remote sensing algorithms (Behrenfeld and Falkowski, 1997a; Banse, 2002). This method was used to measure net primary production NPP, but has been gradually replaced by the use of stable isotope additions, including 13C (Slawyk et al., 1977) and 18O (Bender et al., 1987). The highest primary production estimates are derived using the 18O addition method, while 13C likely underestimates NPP rates, at least as compared to 14C (Regaudie-de-Gioux et al., 2014).

In addition, advances in optical satellite technology over the last three decades have facilitated semi real-time estimations of global NPP (total and fractioned) using multiple models (e.g., Behrenfeld and Falkowski, 1997b; Behrenfeld et al., 2005; Hirata et al., 2009). Despite rich biological productivity (2.57 million tons in 2012 according to FAO, 2014), both measured and estimated NPP rates along Chile's continental shelf are few and discontinuous (i.e., see Table 1 in Thiel et al., 2007) and therefore warrant an effort to combine these approaches.

Different temporal variability in coastal upwelling affects NPP along Chile's coast. As upwelling process occurs, subsurface water of equatorial origin (Equatorial SubSurface Water, ESSW) is upwelled to the surface, delivering waters with low temperature and dissolved O2 levels, but high salinity and nutrient content (Silva et al., 2009). Areas characterized by a semi-permanent upwelling (off Coquimbo, 30◦ S) show similar NPP rates throughout the year, while off central Chile (Concepción, 36◦ S) strong seasonal variability is observed (i.e., Table 1 in Thiel et al., 2007).

Off central Chile, upwelling events are mainly driven by South Pacific Subtropical Anticyclone seasonal fluctuations (Fuenzalida et al., 2008) that promote local southern winds concentrated in spring and summer (Sobarzo and Djurfeldt, 2004) and even synoptic variability with 3–12 day fluctuations leading active and relaxed upwelling events (Aguirre et al., 2014). In addition, upwelling events driven by the Biobío canyon (Sobarzo and Djurfeldt, 2004; Sobarzo et al., 2016) have been documented, as well as NPP enhancement due to the effect of geomorphological protection (Montecino et al., 2006; Henríquez et al., 2007).

The coastal area off Concepción is also influenced by river runoff, mainly from the Itata and Biobío rivers (Saldías et al., 2012; **Figure 1**). Seasonal freshwater discharge to the continental shelf (with a peak during the rainy June and July; Sobarzo et al., 2007) may affect the stability/stratification of the water column, as river runoff mainly affects salinity and coastal NPP with the input of substantial amounts of inorganic nutrients (Sánchez et al., 2008).

Furthermore, remote equatorial processes (i.e., El Niño Southern Oscillation, ENSO) affect the productivity and composition of the marine phytoplankton community in the northern Chile, resulting in a shift in thermocline depth and upwelling patterns (Ulloa et al., 2001; Escribano et al., 2004; Iriarte and González, 2004). The factors affecting NPP off central Chile and the effect of processes such as the ENSO cycle (El Niño (EN) and La Niña (LN) events) on the upwelling and the precipitation regime are poorly understood (Montecinos and Aceituno, 2003; Montecinos and Gomez, 2010), and information regarding biological productivity and phytoplankton phenology remains lacking (Corredor-Acosta et al., 2015).

This study investigates how upwelling events, river runoff, ENSO forcing, and other processes modulate NPP patterns in an upwelling region off Central Chile, using in situ NPP measurements and satellite estimated NPP rates. Bio-oceanographic data from the Centro de Investigación Oceanográfica en el Pacífico Sur-Oriental (COPAS) time series station, in situ NPP measurements from 2006 to 2013, and a decade's worth (2006–2015) of satellite NPP estimations (NPPE), were collected, calibrated and compared in order to infer the mechanisms that cause temporal variability in NPP off central Chile. In this way, this study represents one of the first temporal baselines registering mid-term NPP dynamics and will increase our understanding of how natural and anthropogenic changes in costal upwelling affect C, N and other biogeochemical cycles. Climate change is expected to contribute to a reduction in precipitation in central-south Chile with a concomitant decrease in freshwater discharge (Boiser et al., 2016), and to an increase of along coast winds, which will support upwelling off the coast of central Chile (Sydeman et al., 2014). These climate change scenarios and their impacts (some already observable) may affect NPP rates and timing. Consequently, we may observe alterations in the production of organic carbon and the trophic network, and therefore in heterotrophic organisms and fisheries (Neira and Arancibia, 2004), as well as changes to the magnitude of

greenhouse gas exchanges across the air-sea interface (Nevison et al., 2004; Cai et al., 2006).

### MATERIALS AND METHODS

### Study Zone and Sampling Strategy

The study area is located on the widest (approximately 60 km) continental shelf of Chile (excluding Patagonia) and framed by two submarine canyons (**Figure 1**), formed by the Itata and Biobío (Sobarzo and Djurfeldt, 2004). The area between 35 and 37◦ S is characterized by an upwelling front ∼25 km wide (Pizarro, 1999). The location of the upwelling front in relation to the coastline is dependent on the position of the major upwelling centers, Punta Nugurne (35◦ 57′ S) and Punta Lavapié (37◦ 08′ S; **Figure 1**), the width of the continental shelf, and the presence of eddies, filaments, meanders and coastal jet (Sobarzo and Djurfeldt, 2004; Letelier et al., 2009).

In situ observations from January 2006 to January 2016 were provided by the COPAS time series station (ST18), located at the center of the continental shelf, ∼18 nm from shore and a depth of 90 m (**Figure 1**). Monthly cruises (R/V Kay-Kay II) visited the site to collect continuous (temperature, salinity and Photosynthetically Active Radiation, PAR) profiles using a Conductivity Temperature Depth (CTD) probe (SeaBird 19 and SeaBird 25), along with discrete samples for O2, nutrients, Chlorophyll-a (Chl-a) and NPP rates at depths of 2, 5, 10, 15, 20, 30, 40, 50, 65, and 80 m using 10 L Niskin bottles. In addition, estimations were computed for:

### Upwelling Index

We computed the Upwelling Index (UI; m−3 s −1 ) as Ekman zonal transport per 1,000 m of coastline (Bakun, 1973), from 2006 to 2015, using Equation (1):

$$UI = \left(\frac{\tau\_{\mathcal{V}}}{\rho f}\right) \cdot 1000m \tag{1}$$

where τy (Pa; kg m−1 s −2 ) is the mean wind stress meridional component within the box 73–74◦W, 36–39◦ S, ρ represents the mean density of the water column at ST18 (1026.21 kg m−3 ), and f indicates the Coriolis parameter corresponding to the latitude of ST18 (8.67 ∗ 10−5 s −1 ). Daily τy data were obtained from the CoastWatch project (http://coastwatch.noaa.gov). In order to generate the 2006–2015 time series, Level 3 (L3) data were extracted from the Quik Scatterometer (QuikSCAT) mission (spatial resolution = 0.125◦ ) for 1 January 2006 to18 November 2009, and L3 Metop Advanced Scatterometer (ASCAT) data (spatial resolution = 0.25◦ ) was used for 19 November 2009 to 31 December 2015. We estimated a coefficient of determination (r 2 ) of 0.85, p < 0.05 (Fisher's p-value), observed during a period when two satellites overlapped (10 October 2009 until 18 November 2009; y-intercept: 0.0123 and slope: 0.8618).

Monthly standardized anomalies were calculated using Equation (2):

$$Anomaly = \frac{(\mathbf{x}N - \mathbf{X}N)}{\sigma N} \tag{2}$$

where xN is the discrete value at a certain time, XN is the mean climatological value and σN represents the climatological standard deviation.

### Rivers Runoff and Precipitation Rates

Data for Itata River daily runoff were collected at the Coelemu Station (72◦ 41′W, 36◦ 28′ S), and daily rainfall was recorded at a meteorological station at the Carriel Sur Airport (73◦ 4 ′W, 36◦ 46′ S, **Figure 1**). Both datasets were provided by the Chilean Dirección General de Aguas (DGA; http://www.dga.cl/).

### Oceanic Niño Index

ENSO was measured as the Oceanic Niño Index (ONI); data were obtained from the Climate Prediction Center (http://www.cpc. ncep.noaa.gov/). For cumulative analysis and considering that data did not cover the entire period of EN 2015-16, five EN and LN events were selected. We used Monte Carlo experiment (1,000 iterations) to assess statistical significance between cumulative NPPE during EN and LN events (Montecinos and Gomez, 2010).

### Net Primary Production Estimations

NPPE were obtained using the standard Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski, 1997b) and extracted using 8-day average estimations at a spatial resolution of ∼9 km (Ocean Productivity Page, http://www. science.oregonstate.edu/ocean.productivity/). Among many

TABLE 1 | Summary of error statistics, i.e., Absolute Average Error (AAE), Root Mean Square Error (RMSE), Bias Mean Normalized Bias (PBIAS), Mean Absolute Percentage Error (MPE) and r 2 coefficients for different satellite NPP areas.

NPP models, standard VGPM from SeaWiFS was considered as the model with higher correlation (r 2 = 0.66) with 610 in situ NPP measurements calculated for the California Current (Kahru et al., 2009).

In order to obtain representative satellite estimates, a NPPE time series was created with 8-day average estimations for a similar time period as used for in situ measurements. Next, basic statistics parameters were used to compare in situsatellite data for different pixel boxes (Swirgon and Stramska, 2015). Comparison revealed higher correlation between NPP and NPPE, expressed as r 2 , and lower error, within the boundaries of 72◦ 45′ -74◦W, 36–37◦ S (**Table 1**), a zone with a total area of 13∗ 103 km2 and approximately 55% located over the continental shelf. Analysis of variance (ANOVA) and linear regressions were performed using model II reduced major axis regression (Legendre and Legendre, 1998). Mean values for each pixel box were calculated, except in the case of the pixel box nearest to ST18, which used a single estimation. If a mean estimation described at least 33% of the pixels, it was considered valid.

### NPPE Weekly Variability

The magnitude of week-to-week (WW) NPPE variability was calculated using Equation (3):

$$WW = \boldsymbol{w}\_t - \boldsymbol{w}\_{(t+1)}\tag{3}$$

where wt represents the 8-day-average NPPE during time t and w(t+1) is the successive NPPE.

Wavelet analysis was used to establish dominant frequency bands and their recurrence in the time series, decomposing the signal and estimating its spectral characteristics as a function of time (Torrence and Compo, 1998). Continuous Morlet wavelet transform was applied using the MATLAB functions provided by Torrence and Compo (1998; http://paos.colorado.edu/research/ wavelets).

Agglomerative hierarchical cluster analysis (single linkage technique) was used to evaluate the degree of dissimilarity in monthly bio-oceanographic features.

### Photosynthetically Active Radiation

L3 satellite PAR data were obtained (using 8-day composites with a spatial resolution of ∼4 km) from the Aqua Moderate Resolution Imagining Spectroradiometer (MODIS) mission, within the limits 72–76◦W and 35–39◦ S. Photic Layer Depth (PLD) was estimated using 62 in situ instant PAR profiles collected using a Biospherical QSP 2150 sensor deployed together with CTD profiler. PLD was identified as the depth where downwelling photosynthetic photon flux was equal to 1% of the surface value. Satellite-determined PLD was calculated as PLD = ln(0.01)∗K −1 490 (Behrenfeld et al., 2005), using the 8-day composite Diffuse Attenuation Coefficient K490 (with a spatial resolution of 4 km) from the Aqua MODIS mission. The vertical attenuation coefficient for downward irradiance in the upper 20 m was calculated as described by Kirk (1983).

### Phytoplankton Sizes Classes

Calculations from Hirata et al. (2011) were used to estimate the percentage contribution of micro-, nano- and picophytoplankton to total satellite Chl-a estimations. L3 MODIS Chl-a estimations (8-day average at a spatial resolution of 4 km) were analyzed for the area 72◦ 45′ -74◦W, 36–37◦ S by setting the proportions of microphytoplankton >100% at 100%, and negative proportion of picophytoplankton at 0%.

### Net Primary Production Measurements

A total of 54 in situ NPP monthly measurements were obtained from July 2006 to December 2013 at ST18 (**Figure 1**). At least four measurements for the photic layer are available for every month, with 5–9 monthly measurements available from 2008 to 2013. Light Carbon Assimilations Rates (LCAR) were calculated using the 13C stable isotopic technique (Slawyk et al., 1977) on a mooring line. Samples were captured in polycarbonate 0.6 L bottles, starting at the surface and to a maximum depth of 50 m. The 13C tracer concentration was 3.6456 mg 13C mL−1 , equivalent to 0.5 mmol mL−1 , and 0.5 mL of tracer was added to each sample, to achieve approximately 10% enrichment. Immediately following tracer addition and prior to sunrise, bottles were deployed to sampling depths for in situ incubation. Incubation periods of 10–12 h (counted from the time of sunrise) were used to average phytoplankton physiological activities along the day. Incubation was terminated at dusk using gentle vacuum filtration (<100 mm Hg) through precombusted grade GF/F glass filters. Filters were dried at 60◦C for 24 h and then stored at constant temperature until laboratory analysis using continuous-flow isotope ratio mass spectrometry (Finnigan Delta Plus IRMS).

NPP rates were determined as LCAR (mg C m−3 t −1 ) (Raimbault et al., 1999; Fernández et al., 2005) according to Equation (4):

$$
\rho D I^{13} C = \left[ \frac{(\% R\_{POC} - R\_n) \* \left(\frac{POC}{12 \* V\_f}\right)}{\% R\_{DIC}} \right] \* 12 \tag{4}
$$

$$\%R\_{DIC} = \frac{\left(\frac{V^{13}C \ast^{13}DIC}{V\_b}\right) + DIC\_i \ast 0.01112}{DIC\_i - \frac{V^{13}C \ast^{13}DIC}{V\_b}} \ast 100\tag{5}$$

where Vf represents filtered volume and POC is the amount of particulate organic carbon recovered in the filter after incubation and measured by mass spectrometry (µg). The excess enrichment of the tracer after inoculation (To) is indicated by %RDIC, calculated using Equation (5). %RPOC indicates 13C enrichment in the filter after incubation, measured by the tracer mass. In Equation (5), V 13C indicates the volume of 13C added to the sample during the inoculation, while 13DIC accounts for the tracer concentration added to the sample (3.6456 mg 13C mL−1 ). DICi represents the initial amount of dissolved inorganic carbon in the sample before inoculation. For this study, a constant value of 26 mg C L−1 was selected and applied based on previous measurements for the area (Paulmier, unpublished data). Vb is the volume in the incubation flask (0.6 L). Volumetric assimilation rates (LCAR) were integrated for the photic layer throughout the photoperiod in order to estimate integrated NPP rates (g C m−2 d −1 ). In this way, our measurements do not include dark respiration, which may cause overestimation of NPP rates (Regaudie-de-Gioux et al., 2014).

### Physico-chemical Variables

Brunt-Väisälä frequency (BVF) based on temperature and salinity data were obtained using the Ocean Data Viewer (version 4.7.4). Dissolved oxygen was analyzed in triplicate for a 125 mL sample by Winkler titration using a Dosimat 665 with automatic photometric end-point detector (detection limit ∼2 µmol L−1 ). Fixed samples were analyzed between 12 and 24 h after collection and the percent saturation of dissolved oxygen was calculated based on Boyer et al. (1999).

Triplicate samples for nutrients were filtered onboard (using 0.7µm, grade GF-F glass filters for NO− 3 and PO3− 4 and 0.45µm membrane of cellulose acetate for Si(OH)4, stored in 15 mL screw-capped centrifuge tubes, and frozen (at −20◦C) until analysis. Until 2011, nutrient concentrations were determined using standard manual colorimetric techniques (Grasshoff et al., 1983) and a posteriori using the Seal AutoAnalyzer 3 (AA3). We only considered Si(OH)4 data from 2010 to 2016 because of possible contamination due to sampling strategy during the period 2006–2009. Chl-a was measured for triplicate 250 mL water samples, filtered using 0.7µm, grade GF/F glass fiber filters. Samples were frozen at −20◦C and later analyzed using fluorometry (Turner Designs 10AUTM) according to standard procedures (Parsons et al., 1984).

### Upwelling-Driven Nutrient Flux

Upwelling-driven nutrient flux (kmol d−1 ; with 1 kmol = 103 mol) from 2003 to 2016 was calculated according to Equation (6):

$$F = \text{UI} \cdot \text{C} \tag{6}$$

where UI is the mean upwelling index (m3 s −1 ) during the 3 days prior to field survey and C represents the mean concentration of subsurface nutrients (40–80 m; mmol m−3 ) at ST18. UI was calculated as volumetric Ekman zonal transport per 1,000 m of coastline; explained by the width of the Itata River mouth (∼1,000 m) and our goal of comparing the two fluxes over an area of similar order of magnitude. This simple model assumes: (1) a constant wind field with no changes in magnitude and direction (i.e., without considering wind curl and Ekman pumping); (2) an ocean with no horizontal or vertical boundaries (i.e., ignoring upwelling events induced by topography and trapped waves); (3) the replacement of the zonal transport by vertical advection (i.e., without considering horizontal advection); (4) an absence of diapycnal mixing; and 5. downwelling-favorable winds produce nutrient loss from the photic layer.

### Nutrient Export by River Discharge

Nutrient flux (kmol d−1 ) from Itata River discharge between 2003 and 2016 was calculated using Equation (7):

$$F = Q \cdot C \tag{7}$$

where Q represents mean Itata runoff (m3 s −1 ) during the 3 days prior to field survey and C is the nutrient concentration (mol m−3 ) at the Coelemu station (**Figure 1**). Approximately trimonthly NO− 3 and PO3− 4 concentrations were obtained from the DGA; NO− 3 was estimated using Kjeldahl methods; PO3− 4 was determined colorimetrically by spectrophotometry. In the absence of Si(OH)4 measurements in the Coelemu station dataset, a well-established stoichiometric ratio, Si(OH)4:NO− 3 = 7.745 (µM/µM; Yévenes, personal communication), was used to estimate Si(OH)4 concentrations from existing NO− 3 data in the Itata River. Our simple approach attempted to compare nutrient input from wind-driven coastal upwelling (NO− 3 , PO3− 4 and Si(OH)4) and from river discharge (freshwater source) in the study zone.

### RESULTS

### Annual Cycles and Their Interannual Variability

### Meteorological and Hydrological Forcing Agents

From 2006 to 2015, the upwelling index (UI) oscillated between −3,324 and 4,369 m3 s −1 1,000 m−1 coastline (**Table 2**; Supplementary Figure 1).The annual cycle based on monthly TABLE 2 | Summary of statistics values, i.e., number of samples (n), Minimum (Min), Maximum (Max), Mean and Standard Deviation (SD) for variables analyzed in annual cycles section. In situ mean values were calculated for the upper 10 m of the water column.

8indicates satellite estimation.

averages of UI time series revealed the predominance of upwelling favorable winds from September to April (66% of the year, with a maximum monthly mean of 1,271 m3 s −1 1,000 m−1 in January; see **Figure 2A**). By contrast, downwelling-favorable winds dominated between June and August and produced a negative UI value with a monthly mean of −87 m3 s −1 1,000 m−1 during August.

High Itata River runoff (>400 m3 s −1 ) was observed from June to September, with a maximum monthly mean of 683 m3 s −1 during August (**Figure 2B**). River discharges of less than 60 m3 s −1 were observed between January and April. The majority of Itata River runoff was explained by precipitation (r 2 = 0.33, p < 0.05, n = 418; **Figure 2C**).

Satellite surface PAR estimations (monthly means) were highest from November to February (>50 E m−2 d −1 ; **Figure 2D**), and lowest (between 14.89 and 22.05 E m−2 d −1 ) from May to August.

#### Oceanographic Variables

Temperature presented a clear seasonal variation (**Figure 3A**). From December to March, when a warming of the thin

surface layer (with values higher than 13◦C), caused by high levels of short wave radiation contrasted with the upwelling of cool subsurface water from the ESSW, created a marked thermal stratification. The intrusion of cold ESSW advected by coastal upwelling and demarked by an 11.2◦C isotherm was clearly observed from September to March, due to a favorable upwelling process (**Figure 3A**). Between May and August the water column presented fairly homogeneous conditions, with higher temperatures with respect to upwelling periods mainly associated with Subantarctic Water (SAAW; Silva et al., 2009).

Seasonal variation in salinity was also observed (**Figure 3B**). Low salinity was detected in a thin surface layer each winter and spring and coincided with high precipitation and river runoff. From June to October, a low salinity area associated with the Itata River plume often extended as far as the zone around ST18, creating notable haline stratification (**Figure 3B**). During 2014 and 2015, the lowest salinities were also observed in the spring (i.e., 20.34 in November 2014 and 25.70 in October

2015; Supplementary Figure 2). Salinity ≥34.0 was observed for the upper 20 m of the water column from June to September, with intermediate surface riverine influence between October and December, and sporadic events from January to April. Higher BVF values were observed between depths of 0 and 20 m, particularly from August to December (**Figure 3C**), and likely due to river discharge; subsurface BVF values appeared to be mainly associated with thermal stratification throughout the year, with the exception of late autumn (April and May). Upwelling of the nearly-homogenous ESSW between October and May caused high density homogeneity and lower BVF values between 30 and 80 m. During austral winter, however, elevated deep stability was observed.

Oxygen distribution in the water column (**Figure 3D**) presented higher O2 concentrations (>200µM) in the upper 15–20 m depth throughout the entire year, but subsurface water exhibited seasonal variation due to the intrusion of poorlyoxygenated ESSW from coastal upwelling. Subsurface hypoxic waters (with O2 levels lower than 50µM) were observable during the upwelling season from September to April; while between May and August this layer was deeper and more oxygenated (**Figure 3D**).

The surface layer showed reduced but not depleted nutrient concentrations, whereas the subsurface layer presented consistently higher nutrient concentrations (**Figures 3E–G**). Within the upper 10 m, NO− 3 and PO3− 4 showed minimum mean concentrations from September to January (with minimums of 6.66 and 0.75µM in October, respectively), intermediate values in wintertime (at 13.74 and 1.41µM in August, respectively) and higher values between February and May (with maximums of 17.46 and 1.96µM, respectively, observed in April). These trends were consistent with the plankton temporal dynamics (see below). In contrast, Si(OH)4 exhibited higher concentrations during late autumn and winter, and lower values (with a minimum mean value of 2.13µM in November) during spring and summer (**Figure 3G**). This pattern reflects the high levels of Si(OH)4 associated with river discharge and the development of diatom blooms in spring and summer.

Mean upwelling-driven NO− 3 and PO3− 4 fluxes were an order of magnitude higher than continental discharges from Itata River (**Table 2**). Conversely, mean oceanic Si(OH)4 flux (1,352 kmol d−1 ) was comparable to river input (853 kmol d−1 ). High standard deviation was partially explained by seasonal variability (**Figure 4**). Indeed, wind-driver nutrient flux was higher from September to April (>1,086, 141, and 1,041 kmol d−1 for NO− 3 , PO3− 4 and Si(OH)4, respectively), as compared with the rest of the year, whereas continental inputs reached maximum values between June and August (>226, 15, and 1,755 kmol d−1 for NO− 3 , PO3− 4 and Si(OH)4, respectively).

### Photic Layer, Chlorophyll-a and Primary Production Rates

Light penetration varied temporally, with a higher vertical attenuation coefficient (0.22 m−1 ) from October to January, as compared with February to April (0.18 m−1 ) and June to August (0.12 m−1 ). In this way, PLD showed a distinctive seasonal pattern, with deepening in winter (to approximately 35 m) as compared to ∼15 m during spring-summer, and intermediate values during late summer-autumn (∼25 m, **Figure 5A**). Climatology for satellite-derived PLD revealed a similar temporal pattern (deeper values from April to August) as compared to in situ measurements (**Figure 5A**), although satellite estimations appeared to overestimate PLD, especially from September to February and in May.

In situ Chl-a concentrations higher than 1 mg m−3 (with a marked subsurface peak through variable depths) were observed from September to May within the upper 20 m of the water column, while lower and more homogeneous values (ranging from 0.5 to 1 mg m−3 ) were detected during winter (**Figure 5B**). Maximum monthly means of surface Chl-a were observed during October (12.0 mg m−3 ), with minimum values (of 0.85 mg m−3 ) recorded in July. The highest monthly mean of satellite Chl-a (7.7 mg m−3 ) was registered during December (**Figure 5C**), while the lowest monthly mean (1.4 mg m−3 ) was observed in July.

Microphytoplankton contributed between 18.3 and 96.5% of Chl-a, with the greatest overall contribution to total Chl-a observed during spring-summer (October–January), with monthly mean contributions from 65.5 to 69.2% (**Figure 5D**). The largest phytoplankton size class demonstrated elevated values from February to April (59.1–66.9% Chl-a), decreasing to below 50% in winter. Throughout the annual cycle, the nanophytoplankton community made greater contributions to Chl-a than picophytoplankton, with highest monthly means recorded in July (31.3 and 22.3% Chl-a, respectively).

Throughout the year, LCAR was greater than 20 mg C m−3 d −1 in the surface layer; LCAR greater than 100 mg C m−3 d −1 was observed during the entire upwelling season (**Figure 5E**). LCAR greater than 400 mg C m−3 d −1 was observed in surface waters (5–10 m depth) in spring and early summer (October 2006, December 2007, January, September and November 2009, and January 2011 and 2013; Supplementary Figure 3), and exceptionally high values of 1,721 and 1,536 mg C m−3 d −1 were recorded during November 2006 and 2013 at a 5 m depth, indicating NPP hotspots. Excluding the exceptionally high LCAR event in November 2006, LCAR:Chla mean values of 19.5 and 31.8 mg C mg Chl-a−1 d −1 were obtained for non-favorable and favorable upwelling periods, respectively, indicating that the wintertime phytoplankton community presents comparable photosynthetic efficiency with groups characterizing the upwelling season.

Daily integrated LCAR for the photic layer varied between 0.03 and 18.23 g C m−2 d −1 . NPP measurements showed a clear seasonal cycle, with lower monthly mean from February to August (from 0.22 to 1.81 g C m−2 d −1 , **Figure 5F**) as compared to September-January (1.81 to 6.43 g C m−2 d −1 ). Monthly winter averages between 0.24 (July) and 0.84 g C m−2 d −1 (August) were observed.

Satellite NPPE fluctuated from 0.45 to 9.07 g C m−2 d −1 with the same seasonal variability as observed for Chl-a (**Figure 5F**). Between October and January, NPPE were higher than 3.9 g C m−2 d −1 (with a maximum monthly mean of 5.1 g C m−2 d −1 in December), while from February to June NPPE decreased to 1.1 g C m−2 d −1 . Basal productivity between 1.1 and 1.5 g C m−2 d −1 was observed using satellite estimations during the winter season. In September, mean NPPE rose to 2.9 g C m−2 d −1 with the activation of upwelling-favorable wind.

A linear regression calculated for in situ NPP and integrated Chl-a for the PLD and between satellite NPPE and Chl-a indicated significant determination coefficients of r 2 = 0.50 and r 2 = 0.86, respectively (**Table 3**). The regression for in situ NPP and satellite NPPE data within the limits of 72◦ 45′ -74◦W and 36–37◦ S revealed a r 2 of 0.61 (**Table 3**), but when this data were separated into two periods, higher r 2 values were found during winter (r 2 = 0.70) with respect to the rest of the year (r 2 = 0.24).

The annual NPP cycle off central Chile has been traditionally divided in two oceanographic periods: a favorable and nonfavorable upwelling season based on physical (Montecino et al., 2006; Sobarzo et al., 2007) and biological (Morales et al., 2007) criteria. However, when bio-oceanographic processes and variables, such as NPP, PLD and nutrients, are considered, a hierarchical cluster tree with a cophenetic correlation coefficient of 0.83 separates the annual NPP cycle into three

different periods (see **Figure 6**): a rising stage (high NPP), a falling stage (intermediate NPP) and a basal stage (low or basal NPP). The bio-oceanographic features observed for TABLE 3 | Summary of error statistics, i.e., number of observations (n), r 2 coefficient, and ANOVA F and p-values for biological variables (Chl-a, Chlorophyll-a; NPP, in situ Net Primary Production; NPPE, satellite-derived Net Primary Production Estimation).

\*indicates Chlorophyll-a integrated within the photic layer; 8indicates satellite estimation. Upwelling period: September to April. Winter: May to August.

the study area during these periods are summarized in **Table 4**.

### Interannual Variability of NPP and its Relationship With ENSO

We observed three EN (2006–07, 2009–10, 2015–16) and three LN (2007–08, 2010–11, 2011–12) events of varying magnitude and duration between 2006 and 2015 (**Figure 7**). Minimum ONI measurements (−1.4◦C) were found during January 2008, October 2010 and December 2010, while maximum ONI indexes (+1.3◦C) were observed during December 2009 and January 2010.

Interannual variability in NPPE was also observed, with higher cumulative NPPE during neutral year 2014–15, LN 2010–11, and LN 2011–12 (1.28, 1.21, and 1.15 kg C m−2 yr−1 , respectively; see **Figure 7C**). In contrast, the lowest values were observed during the neutral years 2012–13 and 2013– 14 (0.98 and 1.01 kg C m−2 yr−1 , respectively). Cumulative NPPE during EN and LN years yielded mean values of 1.07 and 1.14 kg C m−2 yr−1 , respectively. Mean NPPE was 1.1 kg C m−2 yr−1 . No statistical significance (Monte Carlo simulations, 1,000 iterations) was determined for NPPE between EN and LN conditions. A NPP wavelet power spectrum and normalized global wavelet spectrum revealed low variance and no statistical significance in interannual variability (**Figure 8**), and identified annual variation as the dominant frequency band.

### DISCUSSION

The continental shelf off Concepción is an important upwelling zone and has been described as one of the most productive areas within the Humboldt Current system (Montecino et al., 1998; Montero et al., 2007; Thiel et al., 2007; Daneri et al., 2012; Iriarte et al., 2012). However, no systematic measurements of NPP have been performed, and NPP rates showed a great variability (e.g., 0.1 and 25.8 g C m−2 d −1 ). Understanding the wide variation in NPP rates, NPP dynamics and factors controlling NPP is vital to improve our knowledge of C, N, and P cycles in central Chile.

At 36◦ S, upwelling occurs over one of the widest continental shelves along the southern part (18–41◦ S) of Humboldt EBUS. This shelf accelerates benthic pelagic coupling and thus the

FIGURE 6 | Dendrogram presenting hierarchical clustering of the degree of similarity between monthly bio-oceanographic characteristics. We used monthly mean surface (0–10 m) values for temperature, salinity, nutrients (nitrate, phosphate and silicate) along with monthly means for satellite and in situ net primary production, satellite and in situ photic layer depth, satellite Chlorophyll-a and upwelling index.

TABLE 4 | Summary of bio-oceanographic conditions during three different periods of the year off central Chile (mean ± standard deviation).

Based on cluster analysis, we considered April as transitional month with high variability. In situ mean values were calculated for the upper 10 m of the water column. 8indicates satellite estimations.

supply of regenerated N, P and even Fe from sediment to bottom water (Farías et al., 1996; Graco et al., 2006; Fennel, 2010), which in turn controls NPP (Ryther and Dunstan, 1971). Thus, the Atlantic EBUS are also characterized by a wider continental shelf and higher NPP as compared to the Pacific, where the Humboldt EBUS presents higher productivity than the California

satellite-derived NPP values during El Niño, La Niña, and Neutral conditions (E). Green points in NPP time series panel indicate in situ measurements.

EBUS (Lachkar and Gruber, 2012); this is partially explained by the presence of continental shelves in northern Peru and central Chile (Thiel et al., 2007). Nevertheless, central Chile is characterized by its intricate geometry, with gulfs, bays, islands and peninsulas generating high time-space NPP variability. A NPP spatial pattern was observed by Testa (2017) and found by Morales et al. (2013) to coincide with Chl-a distribution.

Each upwelling system is therefore considered unique in terms of its productivity. Off central Chile, the magnitude of fertilization in the photic layer appears to be significant and quasi-permanent throughout the year, not only due to the coastal upwelling present during 66% of the year, but also as a result of nutrient influx from rivers (see below). These nutrients can sustain a large population of photosynthetic organisms and in turn present a marked impact on subsequent trophic levels (Quiñones et al., 2010). We estimated mean annual NPP at 1.1 kg C m−2 yr−1 (**Figure 7D**), similar to the value determined by Daneri et al. (2000), who estimated an annual mean of ca. 1 kg C m−2 yr−1 , but lower than estimates for Benguela and Canary EBUS (Lachkar and Gruber, 2012). Cubillos et al. (1998b) estimated the primary production required to sustain central-southern Chile coastal pelagic fisheries at 0.65 kg C m−2 yr−1 , suggesting that no productive capacity limits the regional fisheries within this zone.

### Seasonal Variability in NPP

The biological annual cycle has been statistically divided into three periods (**Figure 6**; **Table 4**):

### Rising Period

From September to January, optimum conditions for the phytoplankton production are observed due to the

synchronization of upwelling activation, optimal surface PAR values, thermal stratification and high subsurface nutrient input. During this period, NPP rates as high as 18.23 g C m−2 d −1 have been observed, with a mean integrated NPP of 4.0 g C m−2 d −1 . The strong light attenuation resulting in shallow PLD is likely to represent the main abiotic factor limiting NPP during the rising period. Nutrients, although occasionally reaching very low concentrations, were never depleted. Satellite estimations reveal a major proportion of Chl-a (and NPP) attributable to microphytoplankton (diatoms) during this period. In situ observations (Montero et al., 2007; Morales et al., 2007; Thiel et al., 2007; Iriarte et al., 2012; Anabalón et al., 2016) have revealed that chainforming diatoms (Skeletonema costatum, Thalassiosira spp. and Chaetoceros spp.) dominate the phytoplankton community during periods of high NPP. A majority of NPP generated by diatoms is consumed by microheterotrophs (Thiel et al., 2007; Vargas et al., 2007), with only a small fraction (0– 10%) of the C fixed by photosynthesis directly assimilated by copepods via the traditional herbivore food chain. Vargas et al. (2006) observed a negative relationship between egg hatching success, naupliar survival and diatom ingestion. Low direct herbivore zooplankton pressure on diatom communities may permit a proliferation of phytoplankton within this system. NPP rates of up to 18.3 (this study), 19.9 (Daneri et al., 2000) and 25.8 g C m−2 d −1 (Montero et al., 2007) have been observed, with the latter two determining the release of O2 during photosynthesis based on changes in bulk O2 concentration.

### Falling Period

During February and March, PAR values remained adequate for supporting the phytoplankton community and the input of nutrients by coastal upwelling remained significant. Nevertheless, in situ NPP strongly decreased, with monthly mean values between 1.0 and 1.8 g C m−2 d −1 . Lower NPP rates and Chl-a values may result from a change in the physiological state of phytoplankton induced by biological (i.e., viral lysis; Brussaard, 2004) or physical components, such as decreases in ocean temperature and surface PAR and elevated turbulence (Garrison and Tang, 2014). The influence of top-down regulation by small ictiofauna species, particularly Engraulis ringens, should be considered. Cubillos et al. (1998a) observed that anchovy reached a maximum relative abundance between February and May, and stomach content analysis demonstrated that this specie is an omnivore and displays a preference for phytoplankton (Arrizaga et al., 1993).

### Basal Period

In situ NPP revealed a mean basal productivity from May to August of 12% of spring-summer productivity, and a NPPE of approximately 24% of spring-summer values for the same period. Basal productivity (averaging 1.01 and 0.46 g C m−2 d −1 for NPPE and in situ measurements, respectively) may be supported by nutrient influx associated with continental runoff and river discharge. Maximum river discharge is correlated with precipitation in winter (Sobarzo et al., 2007). River discharge to the coastal zone produced a plume extending to the fixed station under study (Saldías et al., 2012). Northern wind patterns observed during the study period are known to drive the downwelling of subantartic waters SAAW (Silva et al., 2009) and consequential mixing with fresh water inputs. In this way, SAAW water appears to maintain relatively lower, but continued productivity. Sobarzo et al. (2007) show that wind was too weak to disrupt the buoyancy stratification around the river mouth at any point during the year. Indeed, continental freshwater discharge increased water column stability and inhibited winddriven vertical mixing. As a result, greater residence times in the water column were observed during winter (Montecino et al., 2006). More intense stratification, together with nutrients input from river discharge, may create favorable conditions for the growth of nano- and picophytoplankton communities (Iriarte et al., 2012).

During the wet season, higher stratification and a decrease in phytoplankton biomass (due to a reduced cell size) may create advantageous conditions for light penetration within the water column. In spite of lower surface PAR values, the 1% limit occurs at a greater depth. As a result, more homogenous Chl-a concentrations were observed within the upper 30 m during winter, with NPP resulting from the integration of a deeper layer. Comparable photosynthetic efficiency was found as compared to phytoplankton groups characterizing the upwelling season.In situ measurements of phytoplankton size classes (i.e., Collado-Fabbri et al., 2011; Iriarte et al., 2012; Velasco et al., in preparation) have shown that the empirical equations provided by Hirata et al. (2011) underestimate the nanophytoplankton community and overestimate microphytoplankton communities during the winter season, which can be attributed to elevated satellite Chl-a values within the area. Seasonal cycles of NO− 3 , PO3− 4 and Si(OH)4 appear to be partially controlled by coastal upwelling, continental discharge, alterations in water masses (ESSW and SAAW) and nutrient uptake by phytoplankton groups.

Poor light attenuation near ST18 during winter suggests that particulate matter transported by the river plume may be predominantly deposited within the river terrace; therefore they do not affect the waters of the continental shelf (Saldías et al., 2012). However, inorganic nutrients from riverine sources seem to reach a greater area, as shown by high near-surface Si(OH)4 concentrations observed at ST18 (located ∼27 km from the Itata River mouth). The study of seasonal continental and oceanic nutrient fluxes in this zone (**Figure 4**) reveals that during the wet season the Itata River is a primary source of Si(OH)4 (11.4 times higher than the contribution of wind-driven costal upwelling). Hickey and Banas (2008) suggested that watershedderived nutrients may help sustain productivity of the coastal ecosystem around the Columbia River during periods of weak upwelling-favorable winds and during periods of downwelling. These results indicate a smaller scale but similar dynamic within our study zone, and support an understanding of the zone off central Chile as one of the most productive in the world, even under non-upwelling conditions (Montecino et al., 2004). Hutchins et al. (2002) and Messié and Chavez (2015) identified an iron limitation on NPP in the northern zone of the Humboldt EBUS, particularly during austral winter. Itata River discharge may represent an important source of dissolved iron to our study zone (DGA, 2004; Salamanca and Pantoja, 2009), while helping to sustain high NPP near the river mouth during the non-upwelling season.

Surface salinity climatology (**Figure 3B**) differed from the results presented by Sobarzo et al. (2007), especially between September and December. We attributed this difference to surface salinity dynamics during 2014 and 2015 (e.g., surface salinity fluctuated between 20.34 and 32.96 from August to November 2014; Supplementary Figure 2). We recognized this period as anomalous, since 2014 and 2015 deviate from the observed trend of water column cooling and sea surface salinity increase at ST18, especially during the austral winter (Schneider et al., 2017). Cumulative Itata River runoff data from August to November 2014 reached decadal highs (Supplementary Figure 4). As a result, the cumulative NPPE value for November 2014 was one of the lowest of the decade (**Figure 7C**). Changes to river runoff timing may modify the division of the annual biological cycle as observed by cluster analysis. A potential variation in basal productivity magnitude and timing in river-influenced zones over the continental shelf plays a fundamental ecological role, as a greater abundance of E. ringens and Strangomera bentincki in early life stages (eggs and larvae) was observed during this season and at salinities between 32 and 34 (Soto-Mendoza et al., 2010). Additionally, increased stratification of the water column may favor a greater downward carbon flux (605 mg C m−2 d −1 ) as compared to other periods of the year (between 152 and 268 mg C m−2 d −1 ; Vargas et al., 2007).

### NPP Hot Moments and Weekly Variability

This is the first study of NPP variability over an extended time period in this area, combining measurements and estimates to complement established understanding of coastal upwelling (Montero et al., 2007; Sobarzo et al., 2007). The continental shelf is also subjected to other kinds of variability, and one cause of misalignments between measurements and estimates of NPP are the high frequency variations in the relaxed and active phases of upwelling events, which cause variations in the microbiological community and production rates. In situ NPP measurements for ST18 reveal the occurrence of events where LCAR were highly increased (or "hot moments" as described by McClain et al., 2003). Peaks as high as 18.2 g C m−2 d −1 (October 2006) do not appear to be reflected by NPPE. This may be explained by the fact that NPP peaks lasted less than 8 days and/or covered an area of less than 81 km2 , and NPPE for the closest pixel to ST18 is unable to capture this level of variability; in general the hot moment's phenomenon is poorly documented. Farías et al. (2015) observed that nitrous oxide hot moments in central Chile are positively correlated with phytoplankton biomass, although these events do not specifically occur with high wind stress or with colder water. Effective Chl-a accumulation requires an optimal window of 3–7 days of reduced winds (Farías et al., 2015). It may be necessary to carry out daily, weekly and intra-seasonal investigations to identify NPP hot moments and avoid underestimations of NPP for this region during the upwelling season. Synoptic processes including active and relaxed coastal upwelling events and coastal waves may result in poor correlation between in situ and satellite-derived NPP during the upwelling season (**Table 3**). Monthly measurements are unable to capture variability, and therefore generate a standard deviation that explains estimated and measured NPP rates. Week-to-week (WW) temporal variability presented in **Figure 9** varies between −5.33 and 3.15 g C m−2 d −1 , with a mean value ± SD of 0.72 ± 0.74 g C m−2 d −1 (n = 454). Absolute magnitude for maximum (>1.1 g C m−2 d −1 ) WW monthly mean variability was observed from November to January, while lower (<0.32 g C m−2 d −1 ) estimations were calculated between May and July (**Figure 9B**). Daneri et al. (2012) calculated up to a six-fold difference in WW gross primary production at a near-coast station within our study zone during the austral summer. We consider that a higher standard deviation in NPP and subsurface nutrients during spring-summer may also be caused by intense mesoscale activity (Morales et al., 2007), since a higher number of oceanic eddies were observed in the Peru-Chile current system for this period (Letelier, 2010). On the other hand, lower WW variability in NPP during the austral winter (**Figure 9**) may explain the higher coefficient of determination between in situ and satellite NPP (**Table 3**) for this period.

### ENSO Effect on NPP Interannual Variability

Modest interannual variability was detected over the study period, with the highest NPPE variability observed during neutral years between 2006 and 2015 (**Figure 7**). The maximum cumulative PAR observed in 2014–15 (Supplementary Figure 4) may explain the highest observed NPPE values. On the other hand, minimum cumulative NPPE observed during 2012–2013 may be explained by a deficit in subsurface nutrient input (Supplementary Figure 2).

No statistical significance was found for NPPE between EN and LN conditions at 36◦ S. Iriarte and González (2004) suggested that the ENSO cycle affected NPP and water column properties in northern Chile (23◦ S) during the EN 1997–98, as compared with the non-EN period, driving a shift toward a major proportion of pico- and nanophytoplankton and decreasing NPP rates. Contrarily, Ulloa et al. (2001) and Pizarro et al. (2002) found no significant alteration in phytoplankton biomass in coastal waters close to this latitude during EN 1998. A similar ambiguity was observed at 36◦ S during the 2006–2015 ENSO cycles, which were less strong than previous events (i.e., EN 1997–98). In situ monthly surveys struggle to effectively capture the ENSO effect on bio-physicochemical properties in central Chile, while the use of remote sensing data may be a useful tool to increase the robustness of these observations. Corredor-Acosta et al. (2015) revealed that the effect of ENSO-related climate variability on Chl-a was relatively weak in the coastal zone, and this study observed no statistical significance for NPPE between EN and LN conditions, perhaps due to the distance from the equator. Nevertheless, a study over a longer time series and covering more intense EN and LN events is recommended in order to improve understanding of the ENSO effect on NPP off the coast of central Chile.

### CONCLUSION

The analysis of a decade of bio-oceanographic information allowed for the identification of the main seasonal patterns modulating the annual biological cycle, which can be divided into three stages (i.e., rising, falling and basal).

The magnitude of fertilization in the photic layer appears to be significant and quasi-permanent throughout the year, not only due to the coastal upwelling present during 66% (from September to April) of the year, but also as a result of nutrient influx from rivers (from May to August). Specifically, the Itata River discharge during austral winter may partially support basal productivity around the river mouth, as a result of nutrient input and the influence of the fresh water plume on water column stratification.

Weekly and synoptic variability intensely influence the NPP annual cycle, especially during the upwelling season, which may explain the poor correlation between in situ and satellite-derived NPP for this season.

The ENSO cycle does not appear to impact annual NPP estimates for the study area, which was estimated at 1.1 kg C m−2 yr−1 .

This research presents relevant insights regarding NPP temporal variability in an upwelling area. The NPP baseline provided by this study may be useful in gaining improved understanding of how NPP responds to natural and humaninduced changes in wind patterns and rivers runoff regimes. Changes in the NPP annual cycle, magnitude and timing may alter the carbon cycle and heterotrophic population dynamics and therefore present important implications for biogeochemical and ecological equilibriums.

## AUTHOR CONTRIBUTIONS

LF provided in situ biogeochemical data from the fixed time series station; IM provided the VGPM estimations; GT supplied satellite data and the three authors wrote and analyzed data for the manuscript.

### FUNDING

This research was funded by the Center of Excellence FONDAP-CONYCIT N◦ 15110009, FONDECYT N◦ 1161138 (LF) and GT was supported by CONICYT-PFCHA/Doctorado Nacional/2017-21170561.

### REFERENCES

### ACKNOWLEDGMENTS

The authors thank the biogeochemistry group of the Universidad de Concepción for their help and support. We also thank F. Tapia, M. Sobarzo, C. Morales, R. Escribano, W. Schneider, and reviewers for their valuable comments and suggestions.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars. 2018.00179/full#supplementary-material

wind intensification in coastal upwelling ecosystems. Science 345, 77–80. doi: 10.1126/science.1251635

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Testa, Masotti and Farías. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# First Record of Microbiomes of Sponges Collected From the Persian Gulf, Using Tag Pyrosequencing

Akram Najafi1 , Maryam Moradinasab2 and Iraj Nabipour 1 \*

1 The Persian Gulf Marine Biotechnology Research Center, The Persian Gulf Biomedical Sciences Research Institute, Bushehr University of Medical Sciences, Bushehr, Iran, 2 The Persian Gulf Tropical Medicine Research Center, The Persian Gulf Biomedical Sciences Research Institute, Bushehr University of Medical Sciences, Bushehr, Iran

The Persian Gulf is a special habitat of marine sponges whose bacterial communities are under-investigated. Recently, next-generation sequencing technology has comprehensively improved the knowledge of marine sponge-associated bacteria. For the first time, this study aimed to evaluate the diversity of the Persian Gulf sponge-associated bacteria using tag pyrosequencing in Iran. In this study, 10 sponge samples from 6 different taxonomic orders were collected from the Persian Gulf using SCUBA diving. The diversity of the bacteria associated with the marine sponges was investigated using the 16S rRNA gene PCR-tagged pyrosequencing method. A total of 68,628 high-quality sequences were obtained and clustered at a 97% similarity into 724 unique operational taxonomic units (OTUs), representing 17 bacterial phyla. Cyanobacteria was the most abundant phylum in the sponges, followed by Proteobacteria, Chloroflexi, Acidobacteria, and Actinobacteria. Other phyla were detected as minor groups of bacteria. Bacterial community richness, Shannon, and Simpson indices revealed the highest diversity in sponge S11 (Dictyoceratida sp.) compared to other sponges. This study showed a diverse structure of bacterial communities associated with the Persian Gulf sponges. The dominance of Cyanobacteria may suggest an ecological importance of this phylum in the Persian Gulf sponges.

Keywords: sponge, symbionts, bacterial diversity, 454 pyrosequencing, the Persian Gulf

## INTRODUCTION

Marine sponges (phylum Porifera) are known as the oldest multicellular animals (metazoans) (found more than 600 million years ago) (Lee et al., 2011; Verhoeven et al., 2017) and represent ecologically important reef builders in benthic communities' worldwide (Bayer et al., 2014; Graça et al., 2015). In the last decade, sponges have attracted research interests because of their symbiotic relationships with a wide range of microbial communities including bacteria (Lee et al., 2011; Giles et al., 2013; Gao et al., 2014a), archaea (Zhang et al., 2014), Cyanobacteria (Gao et al., 2017; Regueiras et al., 2017), and fungi (Maldonado et al., 2005). The sponge-associated microorganisms can constitute up to 40–60% of the total sponge biomass (Gao et al., 2014a; Graça et al., 2015; Gaikwad et al., 2016). They may play crucial roles in sponge survival in the marine ecosystem including recycling of nutrients such as nitrogen and sulfur (Montalvo et al., 2014; Zhang et al., 2014), removing metabolic waste (Jackson et al., 2012), and producing bioactive secondary metabolites (Graça et al., 2015).

#### Edited by:

Alejandro A. Murillo, EMBL Heidelberg, Germany

#### Reviewed by:

David William Waite, University of Auckland, New Zealand Fabiano Thompson, Instituto Alberto Luiz Coimbra de Pós-Graduação e Pesquisa de Engenharia (COPPE), Brazil

> \*Correspondence: Iraj Nabipour<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 05 February 2018 Accepted: 18 June 2018 Published: 06 July 2018

#### Citation:

Najafi A, Moradinasab M and Nabipour I (2018) First Record of Microbiomes of Sponges Collected From the Persian Gulf, Using Tag Pyrosequencing. Front. Microbiol. 9:1500. doi: 10.3389/fmicb.2018.01500

Extensive studies have been conducted to investigate the bacterial communities associated with different sponge species using both culture-dependent and culture-independent techniques (Giles et al., 2013; Jeong et al., 2014). As is common with most environments, <1% of bacteria present in sponge tissues can be successfully cultivated (Nam et al., 2011; Jackson et al., 2012). In the past two decades, culture-independent methods (mainly based on 16S rRNA gene) have led to a deeper understanding of the microbial diversity in sponges (Alex and Antunes, 2015). Numerous sponge-associated bacteria have been identified using culture-independent techniques such as denaturing gradient gel electrophoresis (DGGE) (Li et al., 2006), fluorescent in situ hybridization (FISH) (Bayer et al., 2014), terminal restriction fragment length polymorphism (TRFLP) analyses (Zhang et al., 2006; Lee et al., 2011), PCR cloning and sequencing (Kennedy et al., 2008). The 454 tag pyrosequencing is a next-generation sequencing (NGS) technology that provides a faster and simpler way to analyze the microbial communities associated with marine sponges (Webster et al., 2010; Lee et al., 2011; Schmitt et al., 2012; Gao et al., 2014a; Gaikwad et al., 2016). This new method enables hundreds of thousands of nucleotide sequences from multiple samples to be examined in a single 10 h reaction (Lee et al., 2011; Nam et al., 2011; Jeong et al., 2013, 2014).

The next-generation sequencing techniques have revealed the presence of more than 25 different bacterial phyla and 2 archaeal lineages in marine sponges around the world (Moitinho-Silva et al., 2014; Rodríguez-Marconi et al., 2015; Verhoeven et al., 2017). Members of the phyla Actinobacteria, Acidobacteria, Cyanobacteria, Chloroflexi, Proteobacteria, Bacteroidetes, and Firmicutes have been described in association with different marine sponges (Jackson et al., 2012; Giles et al., 2013; Bayer et al., 2014). However, the marine bacterial communities can vary in different sponges with respect to both microbial richness and diversity.

The Persian Gulf, a small, shallow, semi-enclosed body of water bordered by the Arabian Peninsula and Iran, is a unique and greatly underexplored marine ecosystem. There are about 55 sponge genera recorded in the Persian Gulf (Najafi, 2012). However, next-generation sequencing technology has not been applied to identify the sponge-associated bacterial communities in Iran and the Middle East.

The present study was aimed at characterizing the bacterial community associated with the marine sponge species collected from the Persian Gulf, Iran using 454 pyrosequencing.

## MATERIALS AND METHODS

### Sponge Sampling

Sponge sampling was performed in May to September 2016 at a depth of 2–3 m offshore Bushehr, Persian Gulf, Iran by SCUBA diving. In this study, the sponges were living in an area where it was exposed to light. Sponge samples were placed in sterile plastic Ziploc bags containing seawater and immediately transported to the Persian Gulf of Marine Biotechnology Research Center. Sponge tissues were rinsed with 0.22-µm-membrane-filtered seawater (FSW) to remove exogenous materials and loosely attached microbes (Jackson et al., 2012). The samples were stored at −80◦C until further processing.

### Ethics Statement

In this study, the sponges collected did not involve endangered or protected sponge species. No specific scientific research permission was required to collect the sponges from the Persian Gulf.

### Sponge Identification

Sponge taxonomic identifications were confirmed with a combination of multilocus DNA markers. The cytochrome oxidase subunit 1 (COI) and partial 28S rDNA fragments (ITS) were amplified using specific primers as previously reported (Becking et al., 2013). PCR amplifications were carried out on a thermal cycler PeQlab, peqSTAR 96X Universal Gradient, Germany under the following conditions: 94◦C for 30 s; followed by 35 cycles of 94◦C for 5 s; 50◦C for 5 s; 72◦C for 12 s; followed by 72◦C for 1 min (Becking et al., 2013). PCR products were purified and sequenced by Macrogen Inc. (Seoul, Korea).

Also, morphological and spicule examination was carried out by Dr. Yusheng M. Huang (National Penghu University of Science and Technology, Taiwan). List of sponge species collected from different locations of the Persian Gulf is shown in **Table 1**.

### Metagenomic DNA Extraction From Sponges

Frozen sponge tissues were defrosted and washed with sterilized and filtered seawater. Then, they were cut into small pieces (about 1 cm3 ) and ground to fine powder under liquid N2 using a sterile pestle and mortar (Jackson et al., 2012; Bayer et al., 2014). DNA was extracted using a hexadecyltrimethylammonium bromide (CTAB) method. Briefly, a subsample of approximately 100 mg of each sponge tissue was suspended in lysis buffer [100 mM Tris, 100 mM EDTA, 1.5 M NaCl (w/v), 1% CTAB (w/v), 2% SDS (w/v)]. Then they were disrupted in the presence of factors such as proteinase K and extraction buffer containing chloroform: isoamyl alcohol (24:1). Samples were washed with a solution of phenol-chloroform in a few steps for DNA purification. Finally, DNA was precipitated with sodium acetate (3 M, pH 5.2) and isopropanol, then washed in 70% ethanol, dried and redissolved in TE buffer (Jackson et al., 2012; Schmitt et al., 2012). Metagenomic DNA was qualified by agarose gel electrophoresis [1% w/v agarose in Tris-acetate-EDTA (TAE) buffer]. The quantitative assessment of the isolated DNA was carried out using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The high-quality DNA was stored at −20◦C until use (Schmitt et al., 2012; Jasmin et al., 2015).

### Pyrosequencing of Barcoded 16S rRNA Gene Amplicons

In this study, the universal primers 27F and 518R were used to amplify a ∼ 400 bp fragment of the bacterial 16S rRNA gene targeting the V1 to V3 hyper-variable regions. These regions were amplified using primer sets (V1-27F: 5′ - X-MID-GAGTTTGATCMTGGCTCAG-3′ and V3-518R: 5′ -X-MID-WTTACCGCGGCTGCTGG-3′ ), in which X indicates the

TABLE 1 | List of sponge species collected from different locations of the Persian Gulf.

adapter sequences and MID (multiplex identifier) shows the different oligomers comprised of 10 nucleotides to tag different samples for barcoded pyrosequencing (**Table 2**; Jeong et al., 2013, 2014). This approach allowed for the mixing of multiple samples in parallel and re-sorting the sequences into order (Nam et al., 2011; Jackson et al., 2012). PCR amplification was performed in a volume of 50 µL containing 3 mM of MgCl2, 0.2 mM of dNTPs, 2.5 U of Pfu Turbo DNA polymerase (Stratagene, La Jolla, CA, USA), 1X Pfu reaction buffer, 0.1µM of each pair of barcoded primers, and 20 ng of metagenomic DNA (Lee et al., 2011; Gao et al., 2014a). PCR was conducted in a Thermal Cycler (Applied Biosystems ABI Perkin Elmer 9600 GeneAmp) using the following conditions: an initial denaturation at 95◦C for 5 min, followed by 35 cycles of denaturation at 95◦C for 60 s, annealing at 55◦C for 60 s, extension at 72◦C for 60 s, and a final extension at 72◦C for 5 min (White et al., 2012). The PCR amplicon libraries were purified using the NucleoSpin R Gel and PCR Clean-up (Macherey-Nagel, Germany) and quantified using NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Pyrosequencing was performed through a GS FLX Titanium system (454 Life Sciences) according to the manufacturer's instructions (Roche, Germany) by a commercial sequencing provider (Macrogen, Seoul, Korea).

### Processing of 454 Tag Sequences Data

In the present study, the sequences generated from pyrosequencing were analyzed as previously reported (Jeong et al., 2014). The low-quality sequences were filtered from the raw reads using Trimmomatic v0.30 (Gaikwad et al., 2016). Briefly, sequences with a read length of less than 172 bp or with mismatches on primer or barcode, a quality score of less than 25 with ambiguous bases N and homopolymers longer than 6 nucleotides were removed from further analysis (Gao et al., 2014a; Moitinho-Silva et al., 2014; Rodríguez-Marconi et al., 2015; Gaikwad et al., 2016). Final sequences that passed the quality checks were then assigned to respective samples based on the barcodes and subjected to Denoiser to increase the accuracy of the sequence processing (Gao et al., 2014a). The sequences were merged into one file and clustered into operational taxonomic units (OTUs) using the Quantitative Insights Into Microbial Ecology Toolkit (QIIME), version 1.8.0. (Caporaso et al., 2010). Chloroplast and mitochondria sequences and chimeric reads were excluded from downstream analyses. In this study, the chimeric reads were removed through the CD-HIT-OTU program (http://cd-hit.org). Taxonomic assignment of representative OTU sequences was performed using the UCLUST (version 1.2.22) taxonomy assigner method (Edgar, 2010) against the SILVA release 119 database as a reference. To have more information about the strains, each OTU was compared to the most closely related 16S rRNA gene sequences from the NCBI nucleotide databases using BLAST search. Phylogenetic analysis was inferred by using the Maximum Likelihood method based on the Tamura-Nei model (Tamura and Nei, 1993). Also, the evolutionary analysis was conducted using MEGA version 7.0 (Kumar et al., 2016).

### Analysis of Microbial Richness and Diversity

Taxonomic abundance was presented in the phylum, class, order, family, and genus. Alpha diversity metrics were computed

TABLE 2 | The MID barcodes for the amplification of 16S rRNA genes.

(observed species, Chao1, Shannon and inverse Simpson) among sponge samples using the QIIME package with a step size of 100 and 100 repetitions per step. These indices were presented to evaluate the richness and evenness of the associated bacteria within each sponge sample (Caporaso et al., 2010). To show whether the number of reads used in the analysis was sufficient in identifying species/OTU, rarefaction curves were calculated using the QIIME script alpha\_rarefaction.py. (Caporaso et al., 2010). Good's coverage index was calculated as C = 1-(s/n), where "s" is the number of unique OTUs and "n" is the number of individuals in the sample (Naim et al., 2014). The beta diversity among the microbial communities in different sponges was evaluated using UniFrac analysis and the QIIME package (Caporaso et al., 2010). The phylogenetic tree was constructed with the FastTree program to cluster the samples by an unweighted-pair group method with arithmetic mean (UPGMA) using average linkages (Nam et al., 2011; Giles et al., 2013). Also, principal coordinate analysis (PCoA) plots were provided using the QIIME to visualize the effect of the microbial community on structuring the diversity in different sponges (Caporaso et al., 2010).

### RESULTS

### Sponge Taxonomic Identification

In this study, 10 sponge samples were collected from three different sites in the Persian Gulf to evaluate their bacterial diversity (**Table 1**). Taxonomic identification of the sponges using a combination of multilocus DNA markers (28S rDNA and COI mtDNA) showed that the sponges belong to families Suberitidae, Pseudoceratinidae, Chondrillidae, Chalinidae, Halichondriidae, unclassified Dictyoceratida, and Irciniidae (**Table 1**).

### Bacterial Richness and Diversity Analyses

A combined total of 134,495 raw pyrosequencing reads of the bacterial 16S rRNA gene fragment comprising 52,835,851 bases were obtained from the sponge samples. Trimming and quality filtering of the raw reads derived 68,628 high-quality sequences. These denoised sequences were clustered at a 97% similarity into 724 unique OTUs. The highest number of OTUs was obtained from the sponge Dictyoceratida sp. (S11), representing 165 OTUs and minimum of 29 OTUs from the sponge Suberites diversicolor (S02).

In this study, 7.43% of 16S rRNA gene fragments were unassigned OTUs at the phylum level. To have more detailed confirmation of the OTU representative sequences corresponding to unassigned, we extracted the sequences from each OTU and performed a BLAST search to check the sequences corresponding to the conserved region of the target region (16S rDNA). But the closest known sequences had less than 89% similarity rate (Supplementary Table 1) and the phylogenetic tree was not informative. Therefore, these OTUs were not included in subsequent analyses. The total number of reads retrieved and OTUs from each sponge are shown in **Table 3**.

Rarefaction plots were constructed based on OTUs at a 97% sequence similarity cut-off value by the QIIME package. Alpha rarefaction curve showed that a reasonable number of reads have been used in analysis and identifying species/OTU. The sequencing depth of the sponge samples indicated the microbial communities were very well sampled (**Figure 1**). However, additional reads may be required to discover more OTUs for the samples such as S02 to show their total bacterial diversity. OTU based alpha diversity measures and Chao1 estimation of species richness revealed the highest richness of bacterial species in the sponge S11 and lowest in the sponge S02 (**Table 3**). Microbial community diversity, Shannon and Simpson indices displayed the highest community diversity and obviously distinguished the sponge S11 when compared to other sponge species collected from the same location (**Table 3**). In this study, mean and s.d. Expected richness (Chao1) was 80.242 ± 37.647, non-parametric Shannon (H′ ) was 3.28 ± 1.070 and Simpson (D) was 0.774 ± 0.141. Also, the obtained average of the Good's coverage index was 99.8 % ± 0.001 for all the sponge species.

### Taxonomic Composition of Bacterial Pyrosequencing Reads

Taxonomic assignment of the sequences of each OTU (68,628) was classified in the domains Bacteria (92.57% of the total dataset). Altogether, 17 bacterial phyla were recovered from the sponge samples. In the present study, 72.96% of the bacterial reads from the sponge samples were affiliated with two dominant phyla, Cyanobacteria (44.22%) and Proteobacteria (α-, β-, and γ- classes) (28.74%) (**Figure 2**). While the other reads belonging to Chloroflexi (8.67%), Acidobacteria (7.13%), Actinobacteria (4.72%), Bacteroidetes (2.04%), Gemmatimonadetes (2.04%), the candidate phylum TM7 (0.81%), Planctomycetes (0.69%), Deferribacteres (0.47%), Nitrospirae (0.19%), Firmicutes (0.10%), BD1-5 (0.07%), Tenericutes (0.06%), Armatimonadetes (0.02%), TM6 (0.02%), and Chlorobi (0.01%) were detected to be the minor groups in the sponge communities. Sponge S11 was obviously distinguished from the other sponges, in terms of containing taxa from 15 different bacterial phyla and candidate phyla. Other sponge samples contained 5 to 10 bacterial and candidate phyla. Also, the candidate phyla BD1-5 and TM6, and some phyla such as Firmicutes, Tenericutes, Armatimonadetes, and Chlorobi were only identified in sponge S11.

Diversity sponge-associated bacteria at the lower taxonomic levels showed that 45 classes and 87 orders were recovered from all datasets. Cyanophyceae (in Cyanobacteria),

TABLE 3 | The observed number of operational taxonomic units (OTUs) and estimations of richness (Chao1) and diversity index (Shannon, Simpson) for 16S rRNA libraries of the Persian Gulf sponge samples.

Gammaproteobacteria, Alphaproteobacteria were the major classes in the sponges, making up 43.78, 15.82, and 7.79%, respectively. The bacterial communities in sponges at the order level were heavily loaded with Synechococcales (43.78%) in Cyanobacteria and Rhodobacterales (16.32%) and Rhizobiales (9.42%) in Proteobacteria. At all the three levels of taxonomic classification, the bacterial community in the sponge S11 was more diverse than those associated with other sponges. Interestingly, sponge S11 has contained more rare bacteria compared with other 9 sponges. In addition, the bacterial compositions in the same sponges (e.g., S05 and S06) from similar sites did vary substantially (**Figure 2**).

### BLASTN and Phylogenetic Relationships of Highly Abundant OTUs

To have more information about the sponge-associated bacteria strains, the representative sequences of each OTU were compared against the nucleotide database in GenBank. In this study, the first 27 OTUs with average proportions of more than 0.5% among sponge samples are shown in **Figure 3**. The most predominant OTU was OTU\_denovo0 in the phylum Cyanobacteria, which accounted for proportions of 16.60% in the sponge-associated bacteria and shared the highest similarity (99%) to Synechococcus rubescens (**Table 4**). The second most dominant OTU was OTU\_denovo2 in the phylum Proteobacteria, which accounted for proportions of 5.30% and shared 93% identity with the closest relative Proteobacteria bacterium Sinobacterium caligoides (**Table 4**). The phylum Proteobacteria was divided into more than 100 OTUs mainly belonging to Gammaproteobacteria, Alphaproteobacteria, and Betaproteobacteria. Following these symbionts in terms of relative abundance were OTU\_ denovo1 in the phylum Cyanobacteria, OTU\_denovo4 in Acidobacteria (read count: 1614), and OTU\_denovo9 in Chloroflexi (read count: 1472) (**Figure 3**). Another highly abundant OUT in the phylum Cyanobacteria was related to OTU\_denovo11 that affiliated with Prochlorococcus marinus with an identity of 97%. Details of all cyanobacterial and proteobacterial sequences included in the study are summarized in **Table 4**. Also, the evolutionary analysis of these two dominant phyla involved 127 nucleotide sequences and are shown in **Figure 4**.

### Clustering of Sponges According to Bacterial Diversity

Weighted UPGMA tree from the Unifrac analysis was constructed to show the relationships between different sponges according to their bacterial communities (**Figure 5**). In this study, sponge S11 was always distinct from the others in all trees. Also, sponges S4 and S9 exhibited the same topology and the closest distance to each other. The association among other sponges was not noticeable in the Unifrac UPGMA tree topologies. Because the bacterial profiles in each sponge were highly different from each other. Further, the difference in the bacterial community structure of the sponge samples was evaluated using Principal coordinate analysis (PCoA) plot based on the unweighted unifrac distance metrics (**Figure 6**). Sponge

FIGURE 2 | Taxonomic classification of bacterial sequences retrieved from different sponge samples collected from the Persian Gulf at the phylum level. Refer to Table 1 for sponge abbreviation.

S11 formed a distinct clade and showed notable differences from other sponges. Three pairs of sponges belonged to the same species (**Table 1**). But each pair was clustered in different clades except for the species Chondrila sp. (S04 and S09). Also, only the sponges S02, S10, and S13 contained Cyanobacteria as a major phylum.

### DISCUSSION

Recently, marine sponges have been a major target of different studies due to their abundant and diverse microbial communities, ecological roles, production of novel bioactive natural compounds, and biotechnological significance (Menezes et al., 2010; Alex and Antunes, 2015). However, our understanding of bacteria-sponge interactions, nature, and diversity of bacteria associated with marine sponges is still incomplete. There are significant gaps in research on the bacterial composition, function, and maintenance of the symbiotic relationships (Menezes et al., 2010; Verhoeven et al., 2017). Deep sequencing approaches such as 454 tag pyrosequencing can be used to explore the microbial diversity of sponges with high efficiency rather than was possible by clone library construction and Sanger sequencing methods (Gao et al., 2014a; Moitinho-Silva et al., 2014).

### Taxonomic Richness and Bacterial Community Diversity

In this study, the number of unique OTUs (724) is in accordance with a study on 12 different marine sponge species from the Atlantic coast, where only 686 OTUs at 97% sequence similarity reported (Alex and Antunes, 2015). Here, the low diversity of OTUs in bacterial communities might be related to the bias of selected primers (V1–V3), in comparison to V3–V4 and V6 regions in other studies (Lee et al., 2011; Gaikwad et al., 2016; Souza et al., 2017). It will be led to inefficient amplification of the bacterial 16S rRNA genes (Gao et al., 2017). On the other hand, differences in sampling depth and the water temperature may affect the number of OTUs in our sponge-associated bacterial communities than other studies (Jackson et al., 2012; Alex and Antunes, 2015; Souza et al., 2017). Another possibility may be because of the difference in using various bioinformatics patterns to analyze the sponge data. In CD-HIT-OTU program, when defining the OUT, there is a step to remove a cluster with fewer reads, such as a singleton or doubleton, as noise. This is not the same OTU picking method in MOTHUR or QIIME. It looks like the total number of OTU appears to be low when removing low size clusters. In the present study, the cutoff X value defined as the low size cluster was applied as 7, and the cluster consisting of less than 7 reads were removed without being picked by OTU. At this time, the total number of reads removed is 46,318.

In this study, alpha rarefaction curve showed that a reasonable number of reads have been used in analysis and identifying species/OTU. The bacterial richness estimation reported in the Persian Gulf sponges ranged from 29 to 165 OTUs. It is concordant with the range observed (29–370 OTUs) in a study targeting 12 marine sponge species sampled from the Atlantic coast (Alex and Antunes, 2015). In contrast, higher bacterial richness (570–3,013 OTUs) was found in the sponge-specific bacterial communities in Irish waters (Jackson et al., 2012) and also in the sponges from the Red Sea (251–444 OTUs) (Moitinho-Silva et al., 2014). Also, (Thomas et al., 2016) showed contributing marine sponges to the total bacterial diversity of the world's oceans, with a bacterial richness accounted for 50–3,820 OTUs in each sponge (Thomas et al., 2016). In the present study, diversity index and number of OTUs were higher in sponge S11 compared to other sponges. Sponge S11 was obviously distinguished from the other sponges, in terms of containing taxa from 15 different bacterial phyla and candidate phyla. At all the three levels of taxonomic classification, the bacterial community in the sponge S11 was more diverse than those associated with other studied sponges. Interestingly, sponge S11 has contained more rare bacteria compared with other sponges.

Taxonomic identification of sponge S11 using a combination of multilocus DNA markers showed that the sponge belongs to the genus of Dictyoceratida sp. Many studies have considered the order Dictyoceratida as the largest producer of new marine natural products, contributing more than 20% of all spongederived novel compounds (Mehbub et al., 2014, 2016). Different bacterial communities associated with the order Dictyoceratida are reported to produce a wide range of natural compounds with a variety of biological activities (Thakur et al., 2005; Mehbub et al., 2014, 2016). Also, Sponge S11 formed a distinct clade in the weighted UPGMA tree from Unifrac analysis and showed notable differences from other sponges.

Interestingly, in this study, the bacterial compositions in the same sponges (S05 and S06) from similar sites did vary substantially. There was a 100 m distance between the sampling sites of these two sponges. This finding confirms the previous studies, indicating that two cohabiting sponges may have different bacterial signatures (Jasmin et al., 2015; Jeong et al., 2015).

In the present study, a small number of OTUs (7.43%) were recorded as unassigned at the phylum level, after quality filtering and removal of chimera. When alignment is performed with UCLUST in the Reference DB, it means that there are no more than 90% references to match (OTU representative sequence). Unassigned OTUs may be a sequencing error element (Chimera, etc.) that could not be removed in the previous step (OTU picking step). Because there is no match result due to high cut off value, 90%. The frequency of unassigned OTUs was also maintained, after re-analysis of the OTU representative sequences corresponding to unassigned in a BLAST search for the conserved region of the target region (16S rDNA). These OTUs were not included in further analyses, but their presence is noteworthy. However, this number was much lower than those reported from Florida (White et al., 2012) and Indonesia (Cleary et al., 2013), where 36 and 34% of OTUs could not be assigned to any bacterial phylum, respectively.

### Community Composition of Sponge-Associated Bacteria

One of the highlights of our study was the high frequency of Cyanobacteria in the Persian Gulf sponges, contributing 44.22% of the total phylum-level diversity. Also, the presence of Cyanobacteria was confirmed in all the sponges studied with various abundance. Cyanobacteria were also the dominant phylum in the sponge samples reported from other tropical and subtropical regions (Alex et al., 2012; Gao et al., 2014a, 2017; Regueiras et al., 2017).

Within this phylum, our sponge samples contained a high proportion (97.37%) of free-living Synechococcus. BLAST search further revealed the dominant OTU in sponge samples and showed a high similarity with previously reported freshwaterspecific species "Synechococcus rubescens" isolated from the deep subalpine lakes (central Europe), Lake Biwa (Japan), Lake Balaton (Hungary), and the Baltic Sea (Ernst et al., 2003). The Cyanobacteria Synechococcus has considered as an autotrophic plankton community and a substantial fraction of marine primary production (Flombaum et al., 2013), as sponge feeding on Cyanobacteria has been extensively confirmed (Pile et al., 1996; Hadas et al., 2009). This genus of marine bacteria has widely distributed in many ocean regions, covering both polar and high-nutrient waters (Flombaum et al., 2013). Furthermore, the cyanobacterial Synechococcus lineage is believed to have several

TABLE

4

of

all

cyanobacterial

and

proteobacterial

sequences

included

in

the

study.

(Continued)

(Continued)

(Continued)

Frontiers in Microbiology | www.frontiersin.org 11

TABLE 4 |

 was determined based on a 16S rRNA gene NA,Notavailable.

TABLE 4 |

ecotypes that are adapted to different environmental conditions including light, temperature, nutrients, and chlorophyll a concentration (Flombaum et al., 2013). It is possible that the core OTUs illustrate bacterial sponge ecotypes that are matched to the niche sponges and are probably environmentally transmitted (Schmitt et al., 2012). As Cyanobacteria are the center of carbon fixation and provide necessary nutrients to photosynthetic sponge hosts (Taylor et al., 2007), the high abundance of Cyanobacteria further indicated the specific roles of photosynthetic bacteria and their profitability in sponge biology (Lemloh et al., 2009; Alex et al., 2012).

The Persian Gulf sponges were recurrently exposed to light in low depth (<3 m) in this study. Therefore, the prominence of photosymbionts was predictable in these sponges. In spite of the important role of Cyanobacteria, they are typically absent in Antarctic sponges (Rodríguez-Marconi et al., 2015). It seems that environmental factors such as temperature, salinity or nutrient levels might impact the composition of bacterial community structures in different sponges (Giles et al., 2013; Cuvelier et al., 2014).

The co-evolution and functional aspects of spongecyanobacteria associations have not been revealed in details. Sponges may acquire their cyanobacterial symbionts by vertical, horizontal or combined transmission routes (Thacker and Freeman, 2012). The genome-level research on cyanobacterial symbionts of sponges showed that they have general and specific adaptations to life within the sponge host in comparison with free-living cyanobacteria (Gao et al., 2014b; Burgsdorf et al., 2015). Sponge symbionts have adapted mechanisms to actively seek out by their sponge hosts (Webster and Thomas, 2016). For example, cyanobacterial symbionts contains large number eukaryotic-like domains, such as ARs. These domains may be involved in avoiding digestion by the sponge host (Gao et al., 2017).

Lifestyle evolutionary and functional studies on other functions enriched and depleted in cyanobacterial symbionts of sponges compared to members of the closely related freeliving strains reveled the precise and smart adaptation of cyanobacteria to live in full of challenge sponge's intercellular environment. Sponge amoebocytes may not actively distinguish between food bacteria and their cyanobacterial symbionts (Webster and Thomas, 2016). Hence, the depleted genes involved in biosynthesis of LPS O antigen in cyanobacterial symbionts of sponges produces a defense mechanism against sponge predation and phage attack (Burgsdorf et al., 2015; Webster and Thomas, 2016). In functional profile studies of bacterial symbionts with sponges, common functions in similar niches were found, indicating functionally convergence of symbionts in the divergent hosts (Fan et al., 2012; Rua et al., 2015). The biological and ecological roles of these functional equivalences may be of general importance for the adaptation of cyanobacterial symbiont to the sponge host environment and other symbiotic interactions.

In the present study, Proteobacteria was the second most abundant phylum in the microbiome of sponge-associated bacteria. Our study is in accordance with other studies showing

proteobacteria as one of the most diverse phyla of spongeisolated bacterial communities, irrespective of the habitat (Menezes et al., 2010; Schmitt et al., 2012; Jeong et al., 2014, 2015; Gao et al., 2017). In Proteobacteria phylum, OTUs were mainly affiliated with Gammaproteobacteria, Alphaproteobacteria, Betaproteobacteria classes, respectively. Multiple studies have highlighted the presence of Gammaproteobacteria in marine invertebrates, such as sponges (Menezes et al., 2010; Giles et al., 2013; Graça et al., 2015; Rodríguez-Marconi et al., 2015), corals (Sun et al., 2014), and oysters (Garnier et al., 2007). In this study, Gammaproteobacteria accounted for 15.82% of the total bacterial community, including mainly isolates belonging to genus Nitrosococcus. Different species of the genus Nitrosococcus are well known as ammonia-oxidizing bacteria (AOB), inducing the process of nitrification in different sponges and removing the ammonia excreted by the sponge host (Hentschel et al., 2006; Gao et al., 2017). The combined action between AOB and nitrite-oxidizing bacteria (NOB), such as members of the phylum Nitrospirae might then be responsible for the conversion of ammonia to nitrate in sponges (Bayer et al., 2007; Han et al., 2013). In the present study, nitrite-oxidizing phylum Nitrospirae (0.19%) was identified as one of the minor phyla in the Persian Gulf sponges. Our study is in accordance to other pyrosequencing studies in which Nitrospirae constituted a small portion of reads, accounted for 0.01–3% among several sponge species (Webster et al., 2010; Lee et al., 2011; Bayer et al., 2014; Gaikwad et al., 2016).

Chloroflexi was ranked as the third most abundant group (8.67%) using our sequencing approach. In this study, sponge 04 (Chondrilla sp.) harbored the largest proportion of Chloroflexi in comparison to other sponges. The presence of this phylum was previously reported in the Mediterranean sponge Chondrilla nucula using a clone library of 16S rRNA gene sequences (Thiel et al., 2007). The Chloroflexi is one of the most common and

diverse bacterial phyla associated with a wide range of sponges, with many sponge-specific lineages detected (Schmitt et al., 2011; Hardoim et al., 2012; Jeong et al., 2013; Gao et al., 2017). Different studies have revealed the important role of Chloroflexi in nutrition, defense (Hardoim et al., 2012) and carbon fixation through photosynthesis in marine sponges (Brück et al., 2010), indicating its autotrophic lifestyle. In the present study, the Chloroflexi OTUs were closer to the autotrophic lineages.

Actinobacteria was another phylum inhabiting the Persian Gulf sponges. This phylum has been widely reported in marine sponges (Schmitt et al., 2012; Giles et al., 2013; Bayer et al., 2014; Cuvelier et al., 2014; Naim et al., 2014). Different studies have been considered Actinobacteria as an important source of bioactive natural products (Izumikawa et al., 2010; Pimentel-Elardo et al., 2010; Abdelmohsen et al., 2012), protecting the sponge hosts against pathogens (O'Connor-Sánchez et al., 2014). There is a possibility that this phylum may provide new opportunities for novel marine drug discovery.

The candidate division TM7, also known as phylum candidatus Saccharibacteria, is a highly ubiquitous and uncultured phylum of bacteria, described through environmental 16S rRNA gene sequence and genome data only (Ferrari et al., 2014). The existence of this phylum has widely been reported from different sponges (Webster et al., 2010; Lee et al., 2011; Schmitt et al., 2012; Gao et al., 2014a; Montalvo et al., 2014; Alex and Antunes, 2015; Gaikwad et al., 2016). In the present study, TM7 was found at very low abundance in some sponge species. It highlights the importance of deep sequencing technology for detection of the sponge-associated uncultivated bacteria and rare microbial groups in sponges. Otherwise, they might have been not discovered by other routine molecular methods (Lee et al., 2011). However, because of the lack of cultivated representatives and minimal genomic sampling knowledge on the metabolism and biological activities of this enigmatic group has been remained unclear (Ferrari et al., 2014).

The candidate phylum Poribacteria is a sponge-specific phylum that has been widely detected and described in a variety of sponge species (Lafi et al., 2009; Cleary et al., 2013). This candidate phylum has typically been reported in high microbial abundance (HMA) sponge microbiomes (Hochmuth et al., 2010) and considered as "indicator species" for these group of sponges (Bayer et al., 2014). However, unexpectedly, no Poribacteria was found in sponges of the present study.

One reason for the lack of Poribacteria in the Persian Gulf sponges could be related to the use of different primers and bioinformatics pipelines (Souza et al., 2017). Also, some studies have shown that Poribacteria has not been associated with the orders such as Halichondrida, Dictyoceratida, and Haplosclerida (Lafi et al., 2009; Jeong et al., 2013). Noteworthy, these orders have been among the six orders found in this study. In addition, in contrast to other studies conducted on sponges (Kennedy et al., 2008; Montalvo et al., 2014; Alex and Antunes, 2015), the Verrucomicrobia and Spirochaetes phyla were not observed in the bacterial community of the Persian Gulf sponges.

Our study showed that more than 90% of the microbial groups observed in the Persian Gulf sponges were also represented in other studies conducted on the seawater samples (Lee et al., 2011; Schmitt et al., 2011; Gao et al., 2014a; Alex and Antunes, 2015). Nitrospirae, Tenericutes, Armatimonadetes, and the two candidate phyla BD1-5 and TM6 were the only bacterial communities reported exclusively in association with the marine sponges.

Our results showed that the representative OTUs sequences were related to sequences mainly from marine environments such as seawater and marine sediments. Interestingly, sponges had a very low proportion of the isolation source in the closest relative bacteria. Also, some sequences were related to other marine invertebrates such as abalone, octocoral, sea cucumber, oyster, ascidian, starfish, etc. This finding supports the hypothesis of possible environmental acquisition and/or horizontal transmission of bacteria (Moitinho-Silva et al., 2014; Alex and Antunes, 2015). It seems that innate immune system in sponges is responsible for differentiation between symbionts and food microbes (Müller and Müller, 2003). As a result, some overlap between the bacteria in the surrounding seawater and marine sponges would occur. A further deep-sequencing approach is needed to improve our knowledge about the nature of the bacterial

### REFERENCES

specificity among the marine sponges (Alex and Antunes, 2015).

### CONCLUSION

In this study, evaluation of 16S rRNA gene amplicon tag pyrosequencing data showed a complex structure of the previously uncharacterized bacterial communities associated with the Persian Gulf sponges (class Demospongiae). OTU-based description of the bacterial communities exhibited altogether 17 different bacterial phyla, containing 14 formally described phyla and three candidate phyla. The dominance of Cyanobacteria may suggest an ecological importance of this phylum in the Persian Gulf sponges. More specifically, most of the bacterial symbionts were previously described as significant participants in the carbon, nitrogen cycle, and chemical defense of the studied sponges. Sponge S11 was highly diverse in comparison to other studied sponges and contained more rare bacteria. More research is needed to fully understand the composition of the Persian Gulf sponge-associated prokaryotic communities, specifically the functionality of these specific microbiomes as an important part of the marine ecosystem.

### AUTHOR CONTRIBUTIONS

IN and AN designed the research. AN and MM performed the experiments. AN analyzed the data. AN and IN wrote the manuscript with contributions from all authors.

### ACKNOWLEDGMENTS

This work was partly supported by the Persian Gulf Marine Biotechnology Research Center, Bushehr University of Medical Sciences under grant number 1667 and the Iran National Science Foundation (Research Chair Award No. 95/INSF/44913). The authors are grateful to Dr. Yusheng Huang for morphological identification of the sponge samples and also, Dr. Stephen A. Jackson for advice and critical reading of the manuscript.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2018.01500/full#supplementary-material

(Rio de Janeiro: Série Livros, Museu Naciona; Universidade Federal do Rio de Janeiro), 165–171.

the influence of host phylogeny and environment. FEMS Microbiol. Ecol. 93: fiw204. doi: 10.1093/femsec/fiw204

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Najafi, Moradinasab and Nabipour. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Streptomyces nigra sp. nov. Is a Novel Actinobacterium Isolated From Mangrove Soil and Exerts a Potent Antitumor Activity in Vitro

Can Chen1 , Yanghui Ye1 , Ruijun Wang1 , Yinglao Zhang2 , Chen Wu3 , Sanjit C. Debnath1 , Zhongjun Ma1 , Jidong Wang4 \* and Min Wu1 \*

1 Laboratory of Marine Microbial Resources Utilization, Ocean College, Institute of Marine Biology, Zhejiang University, Hangzhou, China, 2 Biomedical Research Program, School of Life Sciences, Anhui Agricultural University, Hefei, China, 3 Institute of Hydraulic and Marine Engineering, School of Hydraulic and Environmental Engineering, Zhejiang University of Water Resources and Electric Power, Hangzhou, China, 4 Department of New Drug Screening, Zhejiang Hisun Pharmaceutical Co., Ltd., Taizhou, China

#### Edited by:

Alejandro A. Murillo, EMBL Heidelberg, Germany

### Reviewed by:

Joachim Wink, Helmholtz-Zentrum für Infektionsforschung, Germany Polpass Arul Jose, Hebrew University of Jerusalem, Israel

#### \*Correspondence:

Jidong Wang<EMAIL_ADDRESS>Min Wu<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 07 December 2017 Accepted: 26 June 2018 Published: 18 July 2018

#### Citation:

Chen C, Ye Y, Wang R, Zhang Y, Wu C, Debnath SC, Ma Z, Wang J and Wu M (2018) Streptomyces nigra sp. nov. Is a Novel Actinobacterium Isolated From Mangrove Soil and Exerts a Potent Antitumor Activity in Vitro. Front. Microbiol. 9:1587. doi: 10.3389/fmicb.2018.01587 A new bacterial strain, designated 452T , was isolated from the rhizosphere soil of the mangrove Avicennia marina in China. As determined, its cell wall peptidoglycan contained LL-diaminopimelic acid; MK-9(H8) and MK-9(H6) were the major isoprenoid quinones; and iso-C16:0 (31.3%), anteiso-C15:0 (16.9%), and iso-C15:0 (12.5%) were the major cellular fatty acids (>10.0%). Phylogenetic analysis based on the 16S rRNA gene sequence revealed that strain 452T formed a distinct lineage in the clade of the genus Streptomyces, and was closely related to S. coerulescens DSM 40146T (99.6% sequence identity), S. bellus DSM 40185T (99.5%), and S. coeruleorubidus DSM 41172T (99.3%). The DNA-DNA relatedness between strain 452T and these type strains ranged between 29.3 and 42.3%. Based on the phenotypic, chemotaxonomic, and phylogenetic features, the strain 452T is considered to represent a novel species of the genus Streptomyces, for which the name Streptomyces nigra sp. nov. is proposed. The type strain is 452T (=KCTC 39960T = MCCC 1K03346T ). Further, strain 452T extracts exhibited a pronounced antitumor activity against human cancer cell lines A549, HCT-116, and HepG2, but not against normal human colon cells CCD-18Co. Active substances in the fermentation broth of strain 452T were isolated by bioassay-guided analysis, and then purified using a macroporous resin, silica gel, sephadex LX-20 column, and semi-preparative high-performance liquid chromatography (HPLC). Eight proline-containing diketopiperazines, namely, cyclo(Pro-Ala), cyclo(Pro-Gly), cyclo(Pro-Phe), cyclo(Pro-Met), cyclo(Pro-Val), cyclo(Pro-Leu), cyclo(Pro-Tyr), and cyclo(L-Leu-trans-4-hydroxy-L-Pro), were identified by electrospray ionization mass spectrometry (MS) and nuclear magnetic resonance (NMR). The compounds displayed different levels of cytotoxicity. The highest cytotoxicity was exhibited by cyclo(Pro-Ala) and cyclo(Pro-Met) against A549 cells, and cyclo(Phe-Pro) and cyclo(Pro-Ala) against HCT-116 cells, with average IC50 values equal to 18.5, 27.3, 32.3, and 47.6 µg/mL, respectively. The diversity of diketopiperazines and other chemicals produced by 452T

**43**

was further investigated using gas chromatography (GC)-MS and liquid chromatography (LC)-MS. The analysis revealed 16 types of metabolites with antitumor activity and 16 other types of diketopiperazines. Hence, extracts of the newly identified strain may be used a starting material for the development of antitumor agents.

Keywords: Streptomyces nigra sp. nov., polyphasic taxonomy, antitumor, bioactive metabolites, diketopiperazine

### INTRODUCTION

fmicb-09-01587 July 17, 2018 Time: 18:6 # 2

The genus Streptomyces contains over 770 species with valid published names at the time of writing1 . Numerous species from this genus produce antibiotics and several other biologically important compounds, including herbicides, antiparasitic agents, immunosuppressants, and other compounds of industrial interest (Behal, 2000 ˇ ; Chater, 2006; Bull and Stach, 2007; Schinke et al., 2017; Ser et al., 2017). Streptomyces are widely distributed in terrestrial ecosystems, especially in the soil. Specific marine environments are pronouncedly different from the terrene environments, and stimulate the production of novel metabolites by marine microorganisms and the evolution of countless new microbes (Simister et al., 2016; Dudek et al., 2017; Kamjam et al., 2017; Petro et al., 2017; Picard, 2017). Marine environments are considered to be a potential source for the discovery of novel Streptomyces that produce bioactive compounds (Sunagawa et al., 2015; Christensen and Martin, 2017; Wei et al., 2018). In recent years, over 20 novel species of the genus Streptomyces have been identified in marine and mangrove environments2 . The mangrove occupies millions of hectares of coastal areas; it is one of the most dynamic environments in the world, and the habitat of various flora and fauna of terrestrial, freshwater, and marine species (Law et al., 2017). Recently, utilization of the mangrove microorganism resource has received a lot of attention, subsequently leading to the discovery of novel Streptomyces species (Hong et al., 2009; Lee et al., 2014).

The genus Streptomyces is an important biological resource for the exploration and discovery of antitumor substances.  These drugs are produced by S. peucetius and S. nogalater (Nobili et al., 2009; Zhang Z.G. et al., 2015; Lu et al., 2017). The polypeptide metabolite dactinomycin is an antibiotic with antibacterial and antifungal activity. Its antitumor activity was confirmed by Juliano and Stamp (1978).  Mithramycin is an aromatic polyketide that contains two deoxysugar chains (a disaccharide of two D-olivoses, and a trisaccharide consisting of D-olivose), a D-oliose, and a D-mycarose, which inhibits testicular cancer and is produced by S. argillaceus (Fernandez et al., 1998; Bilyk and Luzhetskyy, 2016; Novakova et al., 2018). Streptomyces also produce 2,5-diketopiperazines (also known as 2,5-dioxopiperazines or cyclic dipeptides), a class of organic molecules, in which the two nitrogen atoms of the 6-membered ring of piperazine form amide linkages (Huang et al., 2014). They are the minimum natural cyclopeptides and possess important biological activities, such as antitumor, antiviral, antifungal, antibacterial, and antihyperglycemic activities (Martins and Carvalho, 2007; Borthwick and Da Costa, 2017). Remarkably, these bacteria also produce mitomycin C and enediyne compounds (including calicheamicins, neocarzinostatin, and lidamycin), which have unique structural characteristics and display pronounced antitumor-activity (Hata et al., 1956; Mao et al., 1999; Liu et al., 2005; Lee, 2009).

In the current study, we isolated a novel member of the genus Streptomyces from the rhizosphere soil of the mangrove Avicennia marina in China, evaluated its antitumor activity in vitro, and characterized the active compounds in its extract.

### MATERIALS AND METHODS

### Sample Collection and Strain Isolation

Strain 452T was isolated from a sample of a rhizosphere soil collected during the spring of 2016 in the mangrove Avicennia marina forest of Zhangzhou (24◦ 20', 117◦ 45'), Fujian Province (China), and stored at 4◦C until use. The soil sample were suspended in sterile water and diluted in a tenfold series. The dilutions were then spread onto modified ZoBell 2216E agar plates and incubated at 28◦C for 5 days. After incubation, a brown colony that secreted dark brown pigment was selected and purified by repeated streaking onto marine agar 2216 (MA; Difco). The strain was named 452T . It was routinely cultured on MA and preserved at −80◦C in marine broth 2,216 (MB; Difco) supplemented with 20% (v/v) glycerol.

### Genomic and Phylogenetic Analyses

Genomic DNA was extracted for amplification in a polymerase chain reaction as described by Hong et al. (2009). The DNA G + C content was determined by reversed-phase highperformance liquid chromatography (HPLC) (Mesbah et al., 1989) using the genomic DNA of Escherichia coli K-12 and salmon sperm DNA (Sigma) as calibration standards. The 16S rRNA gene of strain 452T was amplified using the universal primers 27F (5'-AGAGTTTGATCMTGGCTCAG-3') and 1492R (5'-TACGGYTACCTTGTTACGACTT-3') (Weisburg et al., 1991). The amplification products were cloned into the pMD19-T vector (TaKaRa) and sequenced. The determined 16S rRNA gene sequence (1,490 bases) was analyzed in pairwise sequence