Marine bacterial richness increases towards higher latitudes in the eastern Indian Ocean

We investigated the bacterial community structure in surface waters along a 2500 km transect in the eastern Indian Ocean. Using high throughput sequencing of the 16S rRNA gene, we measured a significant latitudinal increase in bacterial richness from 800 to 1400 operational taxonomic units (OTUs) (42% increase; r2 = 0.65; p < 0.001) from the tropical Timor Sea to the colder temperate waters. Total dissolved inorganic nitrogen, chlorophyll a, phytoplankton community structure, and primary productivity strongly correlated with bacterial richness (all p < 0.01). Our data suggest that primary productivity drives greater bacterial richness. Because, N2‐fixation accounts for up to 50% of new production in this region we tested whether higher N2‐fixation rates are linked to a greater nifH diversity. The nifH diversity was dominated by heterotrophic Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria. We did not find any mechanistic links between nifH amplicon data, bacterial richness, and primary productivity due to the overall low nifH evenness in this region.

Mapping the diversity patterns and latitudinal ranges of organisms dates back to the 19 th century (Humboldt 1816), and has been crucial for understanding evolutionary and ecological processes that shape global biodiversity (Rosenzweig 1995). Generally, the tropics are believed to harbor greater biodiversity than the temperate regions (Barton et al. 2010), yet recent research is questioning this paradigm as research has shown that such gradients are not ubiquitous in the ocean (Chaudhary et al. 2017). The biogeographic coverage of marine microbial data is rather sparse (Ladau et al. 2013), and the existing distribution maps have been extrapolated across large areas through statistical modeling (Fuhrman et al. 2008;Ladau et al. 2013). Mapping microbial biogeography against key abiotic and biotic parameters should in principle allow for a better understanding of how bacterial communities respond to changes in the environment. This will, in turn, shed light on how biogeochemical cycles and ecosystem services, of which bacteria are critical components, are linked.
Geographic gradients of marine diversity are currently explained primarily by two, not mutually exclusive, ecological mechanisms: (1) diversity increases with increasing productivity (more resources support higher numbers of species) and (2) diversity increases with increasing temperature (higher temperatures catalyze the kinetics of biochemical and metabolic pathways; see references within Fuhrman et al. 2008).
A meta-analysis from Hillebrand (2004) proposed that microscopic analysis of unicellular marine eukaryotes showed weak latitudinal diversity gradients, with no expected gradients for bacteria. The underlying mechanisms are that their large population sizes, their high transportability, and their rapid generation times increase the chances to establish populations in new habitats (Finlay 2002). Pommier et al. (2007) and Fuhrman et al. (2008), on the other hand, showed that bacterial richness correlated inversely with latitude. They found that the diversity of many marine taxa increased toward the tropics, with temperature being the main ecological driver for bacterial diversity. The model presented by Ladau et al. (2013), however, predicted from environmental variables that bacterial richness should peak during the boreal and austral winters at temperate latitudes, with most bacterial diversity hotspots located along 308 latitude N and S, extending until 158N and 158S. Their data contrast the temperature driven mechanism postulated by Fuhrman et al. (2008), yet clear mechanism for the seasonal diversity trend remained to be clarified. Gilbert et al. (2012) also suggested that in their data the strong seasonal patterns (with highest bacterial diversity recorded in winter at 50815 0 N) are more important than trophic interactions.
With an estimated abundance of 10 4 -10 6 cells mL 21 , a great genetic diversity, and key roles in biogeochemical processes, marine bacteria are the foundation of the marine food web (Arrigo 2005;Sunagawa et al. 2015). The pelagic N-cycle in particular can constrain marine productivity, including the associated fixation and sequestration of CO 2 , since most of the open ocean is limited by the bio-availability of N (Falkowski 1997). The marine N-cycle encompasses a suite of complex biogeochemical heterotrophic and autotrophic assimilatory and dissimilatory pathways that transform a diverse range of organic and inorganic N compounds (Zehr and Kudela 2011). Heterotrophic N 2 -fixation has been suggested to support up to 50% of the primary productivity in the eastern Indian Ocean (IO; McInnes et al. 2014;Raes et al. 2014). Yet the main N 2 -fixing bacteria who play a key role in delivering new nitrogen remain unknown across the Leeuwin Current (LC) region in the understudied (Hood et al. 2016), eastern IO.
In this manuscript, we test the hypothesis whether (1) bacterial richness along a latitudinal gradient in the eastern IO increases toward the warmer tropical waters and (2) if bacterial richness can be linked to primary productivity. Because N 2 -fixation takes up a prominent role in this ecosystem, we also (3) test whether the diversity of the functional gene for N 2 -fixation (nitrogenase; nifH) increases with increasing N 2 -fixation rates.

Materials and methods
Underway samples from a clean ship's intake at 6 m depth were collected along a 2500 km transect from tropical waters near the Timor Sea ( 128S) to temperate waters ( 348S) aboard the R.V. Southern Surveyor in the southeastern IO from 07 September 2012 to 16 September 2012 ( Fig. 1a; voyage number SS2012T06). Underway temperature and salinity measurements along with biogeochemical observations (including NO 2 3 , NO 2 2 , NH 1 4 , Si, and PO 23 4 concentrations) were obtained along the transect. Dissolved inorganic nutrients were analyzed on a Bran 1 Luebbe AA3 HR segmented flow analyzer by the Commonwealth Scientific and Industrial Research Organization (CSIRO) hydrochemistry group, following standard spectrophotometric methods (Hansen and Koroleff 2009 to 0.01 lmol L 21 , for Si to 0.01 lmol L 21 , and for NH 1 4 to 0.004 lmol L 21 . All the data used to analyze the data for this paper are available at the PANGEA repository at DOI:10.1594/PANGAEA.882511.

16S Polymerase chain reaction (PCR), sequencing and bioinformatics
In this study, we present the first high throughput tagsequencing-based analysis of bacterial communities along the West Australian boundary current, the LC. Two liters of seawater was filtered through 0.22 lm pore Sterivex TM GP filter (Millipore V R , Massachusetts, Cat. # SVGPL10RC), with a peristaltic pump for DNA analysis. DNA was extracted from Sterivex filters following a modified phenol:chloroform:isoamyl-based DNA extraction protocol alongside extraction columns from the PowerWater DNA isolation kit (Mo Bio Laboratories, U.S.A.) (Appleyard et al. 2013). We note that the total microbial abundance, measured using flow cytometry, remained around 1.6 3 10 6 cells mL 21 across the whole transect, without a significant latitudinal gradient, as reported in Raes et al. (2016).
Bacterial diversity was examined via tag sequencing targeting the V1-V3 region 16S rRNA genes with the bacterial forward 27F and reverse 519R primer sets using Illumina MiSeq (Goodfellow and Stackebrandt 1991;Turner et al. 1999). After merging of paired-end reads using FLASH (Magoč and Salzberg 2011) and read quality control, including removal of chimeras (see Bissett et al. 2016) sequences were clustered into OTUs at 97% sequence similarity using the open reference OTU picking pipeline USEARCH 64 bit v8.0.1517 (Edgar 2010) with the UPARSE algorithm (Edgar 2013). See Supporting Information for detailed bioinformatics.

nifH PCR, sequencing and bioinformatics
Amplicons targeting nifH were amplified and sequenced at the Ramaciotti Centre for Genomics (UNSW Sydney, Australia). The nifH2 and nifH1 primers (Zehr and McReynolds 1989) modified with Illumina overhang adapter sequences (see Supporting Information) were used to amplify 359 bp region of the nifH gene. PCR reactions are outlined in the Supporting Information. Raw reads were trimmed (the last 10 bp and 100 bp from R1 and R2, respectively) using the trimfq function of Seqtk (https://github.com/lh3/seqtk) and resulting trimmed read pairs merged with FLASh (Magoč and Salzberg 2011) with a minimum overlap of 50 bp. Merged reads with length < 300 bp or > 380 bp were removed. The open reference OTU picking pipeline USEARCH 64 bit v8.0.1517 (Edgar 2010) was used to detect chimeric sequences with the de novo UCHIME2 algorithm (Edgar 2016) and to generate de novo OTUs at 95% similarity (Penton et al. 2013;Messer et al. 2015) and representative sequences using the UCLUST algorithm (Edgar and Flyvbjerg 2015). The Ribosomal Database Project's (RDP) Fra-meBot ) from the Fungene pipeline (Fish et al. 2013) was used to detect frameshift errors using Hidden Markov Models and the RDP's reference nifH protein database. Representative nifH sequences were aligned (using the ARB phylogenetic software package (Ludwig et al. 2004) using the database from the Zehr research group (http:// pmc.ucsc.edu/~wwwzehr/research/database/). See Supporting Information for more detailed in regards the bioinformatics.

Pigment analysis and primary productivity
Chlorophyll a (Chl a) extractions, along the transect from 6 m depth, were carried out following acidification according to Parsons et al. (2013) on a Turner Trilogy fluorometer. For Chl a analysis, 0.525 L of sample water was gently filtered with a pressure drop < 10 kPa on 25 mm Whatman GF/F filters. Four liters of sample water was also filtered via gentle vacuum filtration on 25 mm Whatman GF/F filters for pigment analysis. Samples were snap frozen and stored in liquid nitrogen until they were analyzed onshore using highperformance liquid chromatography according to the CSIRO methods, see chapter 2 in Hooker et al. (2005). Processed pigment data were analyzed for phytoplankton community composition using the diagnostic pigments of the dominant phytoplankton taxa groups (see Vidussi et al. 2001;Hirata et al. 2008;Aiken et al. 2009). These pigment data have been used to describe shifts in phytoplankton community structures along the LC system in earlier work from Raes et al. (2014). Primary productivity rates were calculated from bulk N assimilation rates (including NO 2 3 and NH 1 4 assimilation rates and N 2 -fixation rates) described in Raes et al. (2014) and unpublished data using a regional C : N ratio of 6.3 (n 5 31). The primary productivity rates along our transect are comparable to the surface production rates reported by Hanson et al. (2005) in the LC. Both pigment and primary productivity data are used in this manuscript to bridge conceptual links between bacterial richness and the eukaryotic primary producers.

Statistics
Bacterial richness was calculated as the number of OTUs observed per sample using the Vegan package in R Studio (Oksanen et al. 2017). The Chao index was also calculated using the SpadeR package in R Studio (Chao et al. 2016) and is used to constrain for the number of rare OTUs at each sampling station. The Chao equation is based on the notion that if a sample contains many singleton or doubleton OTUs, it is likely that more undetected OTUs exist. The Chao equation will consider these potential undetected OTUs and will consequently estimate greater or "true" species richness. Multi-dimensional scaling was done using PRIMER software package (Clarke and Gorley 2015).

Results and discussion
Our continental scale data from the eastern IO showed an increase in bacterial richness with increasing latitude (r 2 5 0.65, p < 0.0001, n 5 14; Fig. 1b). The bacterial richness data along our transect confirms the predicted diversity hotspots around 308 latitude from Ladau et al. (2013) and the data from Chaudhary et al. (2017), where the authors showed a decline in species richness near the Equator. Bacterial richness in the eastern IO peaked at 1400 OTUs around 348S and declined steadily to 800 OTUs toward the tropical waters of the Timor Sea at a rate of 29 OTUs per degree of latitude. Rarefaction curves, showing richness as a function of sampling effort (Fig. 1c), indicated that deeper sequencing efforts would only increase the apparent statistical differences in the observed bacterial richness between tropical and temperate waters. Chao diversity indices supported our previous statements of an increased bacterial richness toward temperate waters. Chao indices showed an expected overall higher bacterial diversity for all stations along the transect, with a lower number of estimated OTUs ( 1100 OTUs) in the tropical waters compared to the temperature waters ( 1600 OTUs; Fig. 1d).
Temperature and salinity along the transect showed an inverse relationship (Fig. 2a). Dissolved inorganic nitrogen (DIN; NO 2 3 1 NO 2 2 1 NH 1 4 ) concentrations decreased significantly from 0.3 lmol L 21 to 0.1 lmol L 21 with latitude toward the colder temperate waters (r 2 5 0.57; p < 0.001, Fig.  2b). The contribution of NO 2 3 to the total DIN pool was about 15% greater in temperate waters, whereas that of NH 1 4 and NO 2 2 increased by up to 20% toward the tropics. This DIN ratio shift has been previously correlated with a change from a pico to a micro and nano-phytoplankton community, an increase in total Chl a concentration (from < 0.2 lg L 21 to 0.6 lg L 21 ), and an increase in primary productivity from 180 lg C m 23 h 21 to 350 lg C m 23 h 21 from tropical to temperate waters . Bacterial richness was highly correlated with colder temperatures (Fig. 1b); higher-salinities (r 2 5 0.49), Chl a values (r 2 5 0.68), microplankton (r 2 5 0.73), nanoflagellate (r 2 5 0.65), abundance, and primary productivity rates (r 2 5 0.40); and lower picoplankton abundance (r 2 5 0.67; distance-based linear model, all p values < 0.01; see Figs. 2c,d, 3a-e).
The global data set from Fuhrman et al. (2008) showed a significant and positive correlation between bacterial richness and decreasing latitude (r 2 0.39) and increasing temperatures (r 2 0.34). Primary productivity and Chl a, however, were not found to be significantly related to bacterial richness on a global scale. Ladau et al. (2013) did show a strong and significant correlation between bacterial richness and distance from thermocline, day length, and PO 2 4 (Spearman r 0.36, 20.64, 0.52, respectively). Their model outputs, however, did not identify significant roles for DIN, Chl a, or primary productivity.
The findings of Pommier et al. (2007), Fuhrman et al. (2008), and Ladau et al. (2013) arise from analysis of global datasets which link physical drivers (such as temperature, salinity, and depth from thermocline) to the main correlations for diversity. These global findings are a good starting point to test which ecological mechanisms force bacterial diversity on finer spatial scales (e.g., in this manuscript, within one regional water mass). In general, high latitude oceanic environments experience low temperatures and high nutrient concentrations. In lower latitudes, a permanent stratified thermocline induces oligotrophic conditions, resulting in lower primary productivity rates (Martiny et al. 2013). However, while these trends may hold for the open ocean, there are many important exceptions to these general patterns associated with productive shelf seas and boundary currents, illuminating more clearly the regional and local mechanisms driving microbial biodiversity. One example is the LC, which carries warm tropical, oligotrophic, low salinity water toward temperate latitudes in the eastern IO (Cresswell 2004). The source waters from the LC spring from the Timor Sea, mix with the Eastern Gyral Current flowing South, and are responsible for suppressing upwelling along the west coast of Australia (Domingues et al. 2007). Here, we show that the bacterial richness in the LC was strongly correlated with lower temperatures, lower DIN, high Chl a, and high primary productivity. The findings from our temporal snapshot from the eastern IO show that richness increased with decreasing temperatures. We suggest that in early spring in this system, the effects of temperature on metabolic kinetics do not control bacterial richness. We propose that, because of the year round presence of the LC in the eastern IO, greater bacterial richness is likely to be driven by increased primary productivity. In previous work from Raes et al. (2014), we postulated that the temperate waters in the eastern IO have a longer food web than tropical waters in the Timor Sea. The longer food web is associated with a larger variety of unicellular eukaryote-derived dissolved organic materials such as exopolysaccharides (Fenchel 2008), which could sustain a larger microbial diversity (Amin et al. 2012). Heterotrophic bacteria are scavengers that utilize organic carbon produced by autotrophs. North, near the equator, we measured a greater abundance of autotrophic picoplankton and larger N 2 -fixers (Trichodesmium; Raes et al. 2014). In the south of this transect, we observed larger auto and mixotrophic planktonic size classes (nano and microplankton) growing on nitrate, coexisting with a relatively uncoupled microbial loop supported independently by heterotrophic N 2 -fixing proteobacteria (Wannicke et al. 2010;Raes et al. 2014). It is an interesting case of niche partitioning between NO 2 3 driven eukaryotic productivity (with higher productivity rates) and N 2 -fixing organisms existing within the same water volume. As a new hypothesis for this region, we propose that it may be this partitioned food web, which supports a greater bacterial diversity.
Claire Horner-Devine et al. (2003) showed that productivity influences the composition and richness of bacterial communities in aquatic mesocosms. The authors noted that a number of bacterial taxonomic groups displayed different richness responses (U-and hump-shaped relationships) along a gradient of primary productivity. Their results and those of others (Gasol and Duarte 2000) suggest that top-down effects are in general less important than bottom-up effects in regards the overall abundance and richness of bacteria (Claire Horner-Devine et al. 2003). The correlation between bacterial diversity and higher primary productivity rates has also been found across environmental gradients such as polar fronts (Wilkins et al. 2013) and oligotrophic gyres in the North Pacific Ocean (Jones et al. 1996). Our data from Chl a concentrations (r 2 5 0.68; p < 0.001); (c) Positive quadratic fit between bacterial richness and microplankton ratio (r 2 5 0.73; p < 0.001); (d) Positive quadratic fit between bacterial richness and nanoflagellate ratio (r 2 5 0.65; p < 0.001); (e) Negative quadratic fit between bacterial richness and picoplankton ratio (r 2 5 0.67; p < 0.001).
the eastern IO support a causal link between productivity and diversity gradients, which can help build a mechanistic understanding of patterns of global marine microbial diversity (Fraser et al. 2015).
Marine primary productivity and the N-cycle are closely linked in the LC system (Waite et al. 2007), as has been shown elsewhere (Zehr and Kudela 2011). In previous work, we showed that N 2 -fixation rate measurements revealed a particularly important role supporting up to 50% of regional primary productivity Raes et al. 2015). These estimates are similar to those reported in the oligotrophic subtropical North Pacific by Karl et al. (1997), where N 2 -fixation could explain > 40% of the new production. Along the transect bulk N 2 -fixation rates declined significantly from the South ( 6 nmol L 21 h 21 ) to the North ( 1 nmol L 21 h 21 ; r 2 5 0.59, p < 0.001; Raes et al. 2014). Our second hypothesis was therefore that higher N 2 -fixation rates in the South are supported by a greater nifH diversity, which consequently would relate to an overall higher bacterial diversity. In our data set, we used amplicon sequencing of the nifH gene to detect 83 unique OTUs (95% similarity; corresponding to 80,586 reads-non-subsampled). Rarefaction curves of the subsampled nifH amplicon reads plateaued in all samples after 400 reads (Supporting Information Fig. 1). Overall, we found three major phylogenetic nifH groups from: (1) mostly conventional nifH genes (Cluster I), (2) some nifH genes from non-molybdenum and non-vanadium nitrogenases (Cluster II), and (3) a variety of sulfate reducing bacteria and nitrogen-fixing bacteria (Cluster III, see Table 1; clusters derived from Zehr et al. 1998;Zehr et al. 2003;Langlois et al. 2015). The nifH sequences from Cluster III (sulfate reducers and clostridia nitrogen-fixing bacteria) have been suggested to be associated with the exoskeleton and/or the gut tract of crustacean zooplankton (Zehr et al. 1998).
The majority of the nifH reads and OTUs in our samples belonged to Cluster I, with representatives of the Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, and unicellular cyanobacteria clades. Our amplicon data showed a low evenness as one OTU from the cluster I nifH group dominated the nifH community (up to 95% at some stations). This dominant nifH OTU belonged to the gamma 4 proteobacteria and was very closely related to Vibrio (to 99% similarity). The subgroup gamma 4 proteobacteria has been found to dominate in other oligotrophic oceanographic regions such as the South Pacific Gyre (Halm et al. 2012). The gamma 4 proteobacteria (HM201363.1) found in Halm et al. (2012) was nearly identical, with a single nucleotide difference over 348 nucleotides, to the sequence we found in the eastern IO. The ubiquitous, likely heterotrophic gamma 4 proteobacteria reads, which made up to 65-95% of the nitrogenase enzyme diversity along our transect, showed no significant relationships with any of the measured parameters (latitude, temperature, salinity, DIN, or primary productivity). A metric multi-dimensional scaling did not show a separation between the temperature and tropical waters as we saw for the bacterial diversity (Supporting Information Cluster I: Conventional nifH genes from proteobacteria (Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria) and cyanobacteria. Cluster II: Alternative nifH genes from non-molybdenum-and non-vanadium containing nitrogenases. Cluster III includes sequences from sulfate reducers including Deltaproteobacteria and clostridia (gram-positive bacteria). Fig. 2). Nor did the overall nifH diversity show any significant trends with latitude, or any other physical or biogeochemical parameters. The higher N 2 -fixation rates did not result in a higher nifH diversity. We could therefore conclude, and reject our 2 nd hypothesis, that nifH diversity did not contribute to the overall increases in bacterial diversity at higher latitudes. The finding and presence of numerous Gammaproteobacteria nifH genes is consistent with the recordings of Gammaproteobacteria 16S rRNA phylotypes in open oceanic waters (Schmidt et al. 1991;Langlois et al. 2015). Messer et al. (2016) noted a seasonal change in the relative abundance of heterotrophic Deltaproteobacteria and Gammaproteobacteria (an increase in winter), along with a decline in N 2 -fixing cyanobacteria such as UCYN-A sequences in the tropical Timor and Coral Sea. Our samples were collected in late winter/early spring and link in with the relative high abundance of heterotrophic Gammaproteobacteria found by Messer et al. (2016) in a similar oligotrophic region. Our findings of Gammaproteobacteria dominating the nifH diversity in the eastern IO compliment the results that Gammaproteobacteria are cosmopolitan diazotrophs found throughout oligotrophic waters (Halm et al. 2012;Langlois et al. 2015).

Conclusion
We propose that the bacterial richness increase in the eastern IO toward temperate and more productive waters is driven in parallel by productive eukaryotes (NO 2 3 based) and heterotrophic proteobacterial N 2 fixers. Although N 2 fixation is an important source of new N in the eastern IO, we did not find a clear mechanistic link between nifH amplicon data, overall bacterial richness or primary productivity. The functional diversity of N-cycling enzymes such as nifH needs to be better explored in this ecosystem since N 2 fixation rates play a pivotal role in supporting regional marine productivity from the tropics to the temperate waters in the eastern IO.
Our data highlight that we need more information on how local environmental forcings can modulate the correlation between bacterial richness and productivity on a regional scale. Better knowledge on the shape of richness vs. productivity relationships and the controls of functional enzymes (such as nifH) will allow us to answer questions such as, how (functional) diversity relate to ecosystem resilience, and how climate change will impact both, in the understudied IO and other marine ecosystems.