Sampling locations and other information for suspended particulate matter (SPM) and surface sediments are shown in Figure 1 and Table 1. SPM samples (n = 18) and surface sediments (n = 8) were collected along a salinity gradient from the lower Pearl River and its estuary to the coastal South China Sea in the summer of 2011. SPM samples were collected from the surface (stations R1–R6) and the bottom (stations R1 and R2) water in the lower Pearl River, from three water depths (surface, middle, and bottom) and during three tidal periods (high tide, slack tide, and low tide) at station M located in the PR estuary, and at four water depths (surface, subsurface, middle, and bottom) at station S in the seawater of the coastal SCS (Figure 1). The depth of the sampling layers in the water column is given in Table 1. About 4–103 liters of water were filtered onto combusted (450°C, overnight) glass-fiber filters (Whatman GF/F, 0.7 μm, 142 mm diameter) using an in situ submersible pump. The pH, temperature, salinity, and depth were determined in situ by a Horiba instrument (W-20XD, Kyoto, Japan; Table 1). Surface sediments (top ca. 10 cm) were collected at all stations (Figure 1; Table 1) using a grab sampler. All samples were frozen immediately in liquid nitrogen and kept at −80°C in the laboratory before analysis.

The SPM samples (n = 18) and surface sediments (n = 8) were freeze-dried and extracted using a modified Bligh and Dyer method (Bligh and Dyer, 1959); the separation of core lipids and intact polar lipids followed the procedure described in Weijers et al. (2011b). Briefly, the total lipid extract (TLE) was obtained by ultrasonic extraction (10 min each, 6 times) of SPM (1 filter) or sediments (5 g) with a single-phase solvent mixture of methanol, dichloromethane (DCM), and phosphate buffer (2:1:0.8, v/v/v; pH 7.4). The TLE was separated over an activated silica gel column eluted with n-hexane/ethyl acetate (1:1, v/v) and methanol for CL and IPL, respectively. For GDGT quantification, a known amount of an internal C46 GDGT standard was added into the CL fraction or IPL fraction (Huguet et al., 2006). The CL fraction was directly measured. The IPL fraction was determined by measuring the CL after cleavage of the polar head groups via acid or base hydrolysis (Pitcher et al., 2009; Weijers et al., 2011b). Briefly, 1/3 IPL fraction (non-hydrolyzed IPL fraction) was directly condensed; another 1/3 IPL fraction was hydrolyzed (2 h) in 1.5 N HCl in methanol, which was called the acid-hydrolyzed IPL fraction (total IPL). DCM and MilliQ water were added, and the DCM fraction was collected (repeated 4 times). The DCM fraction was rinsed (6 times) with MilliQ water in order to remove acid and dried under N₂ gas. The last 1/3 IPL fraction was subjected to base hydrolysis (2 h) in a 1N KOH in methanol/H₂O mixture (95:5, v/v), which was called the base-hydrolyzed IPL fraction (phospho IPL). Together with the condensed CL fraction, the CL fraction, two fractions of IPL-derived core lipids and non-cleaved IPL fraction were dissolved in n-hexane/isopropanol (99:1, v/v), and filtered using PTFE filters (pore diameter of 0.45 μm). The acid-hydrolyzed IPL reflected total IPL, which, prior to hydrolysis, were attached to both phosphatidyl and glycosidic head groups; The base-hydrolyzed IPL represented IPLs with phosphatidyl head groups only (phospho IPL; Weijers et al., 2011b). Analysis of the non-hydrolyzed IPL fraction was performed to determine any carryover of CL into the IPL fraction.

GDGTs from all treatments were analyzed using high performance liquid chromatography-atmospheric pressure chemical ionization-tandem mass spectrometry (HPLC-APCI-MS/MS), which was performed as described by Zhang et al. (2012), using an Agilent 1200 LC equipped with an automatic injector and coupled to a QQQ 6460 MS; peaks were evaluated using Mass Hunter LC-MS manager software. Separation was achieved using a Prevail Cyano column (2.1 × 150 mm, 3 μm; Alltech Deerifled, IL, USA) with n-hexane (solvent A) and a mixture of n-hexane/isopropanol 90/10 (v/v; solvent B). The (M+H)⁺ ions of each core isoprenoid GDGT (m/z 1,302, 1,300, 1,298, 1,296, 1,294, 1,292) were monitored via selected ion monitoring (SIM) mode (Schouten et al., 2007).

Indices based on the fractional abundance of GDGTs were calculated as follows:

(Schouten et al., 2002)

(Zhang et al., 2016)

with the GDGT numbers corresponding to the GDGT structures in Figure S1. Note that RI₂ was originally developed for the current study and modified from RI₁ (Zhang et al., 2016), in which the fractional abundance of crenarchaeol was replaced by GDGT-4 in order to eliminate the influence of crenarchaeol on weighted average number of cyclopentane rings in GDGTs.

The SPM samples (n = 12) from station R1 (river water), station M (mixing water), and station S (seawater) were selected for the DNA analysis. The frozen filters were washed 3 times by phosphate buffered saline (pH 7.4). The supernatants were centrifuged under 11,000 g for 10 min. The DNA was extracted following the protocol of FastDNA SPIN Kit. The DNA samples were dissolved with a final dilution in 100-μL deionized water and preserved at −80°C until further processing. The DNA concentrations were quantified in duplicate with a Nano-Drop spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA). The quantitative PCR primers were Arch_334F (5′ACGGGGCGCAGCAGGCGCGA3′)/Arch_ 518R (5′ATTACCGCGGCTGCTGG3′) for total archaeal 16S gene quantification (Bano et al., 2004), GII-554-f (5′GTCGMTTTTATTGGGCCTAA3′), and Eury806-r (5′CACAGCGTTTACACCTAG3′) for MG-II Euryarchaeota 16S gene quantification (Massana et al., 1997; Teira et al., 2004). The qPCR analysis was performed at 95°C for 30 s and 40 cycles at 94°C for 30 s, 55°C for total archaea and 53°C for MG-II Euryarchaeota for 30 s and 68°C for 1 min. Triplicate measurements were run for each sample and standard.

PCR bands of 16S rRNA gene and MG-II gene were amplified from SPM samples in Station M. They were recovered by a Gel Extraction Kit (omega) and sequenced on the 3730-sequencing platform. The sequences were annotated as the corresponding target genes, which demonstrated the specificity of the chosen qPCR primers. A dilution series of purified DNA from those positive clones were used as standards. A melting curve analysis was performed to demonstrate that the fluorescent signal obtained in a given reaction was consistent with the expected profile for specific PCR products. The R² values of standard curves were >0.99. The efficiency of each qPCR was between 87 and 99%.

SPM samples collected in river water, mixing water and seawater were selected to conduct MiSeq pyrosequencing targeting the archaeal 16S rRNA gene. In contrast to the qPCR primers, these primers targeted longer sequences to increase the precision of phylogenetic analysis. The primers were Arch_344F (5′ACGGGGCGCAGCAGGCGCGA3′) and Arch_915R (5′GTGCTCCCCCGCCAATTCCT3′; Gantner et al., 2011). The MiSeq sequencing was conducted on the MiSeq platform (2 × 250 PE, Illumina) at the Shanghai Personalbio Biotechnology (Shanghai, China). Mothur (version 1.29.2; Schloss et al., 2009) was applied to filter the raw pyrosequencing data. The selected sequences were analyzed using the QIIME standard pipeline (Caporaso et al., 2010). Taxonomy was assigned according to the Ribosomal Database Project (RDP) classifier 2.2 (minimum confidence of 80%; Cole et al., 2009). The GenBank accession numbers are PRJNA38421 for those archaeal 16S rRNA genes.

The satellite-derived SWT was determined with a spatial resolution of 4 km from the NOAA advanced very-high-resolution radiometer (AVHRR; version 5.2; http://www.nodc.noaa.gov/SatelliteData/pathfinder4km/). The June mean SWT was obtained from the daily averaged values of 30 days in June 2011 (sampling month). The annual mean SWT and winter mean SWT represented 8-year mean values of annual mean temperature (2004–2011) and monthly mean temperature (December–February), respectively, as the surface sediment (top ca. 10 cm) collected in this study might represent a deposition of 6–10 years, based on an estimation from Strong et al. (2012).

The TEX₈₆-derived temperature was calculated based on the calibration of Kim et al. [SST = 68.4 × log (TEX₈₆) + 38.6; (Kim et al., 2010)]. The CL-TEX₈₆ temperatures derived from either the SPM or the surface sediments were close to the satellite-based annual mean SWT in the river water and mixing water, station R and M, respectively; whereas CL-TEX₈₆ temperatures for the seawater station S were lower than the winter mean SWT (Figure 2). The correspondence between the CL-TEX₈₆ temperatures derived from SPM and sediment samples (R² = 0.70, P < 0.01) along the salinity gradient indicated that the TEX₈₆ signal in the sediment predominantly reflected archaea from the water column, which is consistent with previous studies in the PR estuary (Wang et al., 2015) and other coastal settings (Herfort et al., 2006; Zell et al., 2014).

Total IPL-TEX₈₆ temperatures in the water column were consistently lower than both the June mean SWT and in situ measurements, although the sampling season was summer (Figure 2). In the mixing water and seawater, the phospho IPL-TEX₈₆ was lower than the total IPL-TEX₈₆ whereas it was very close to the CL-TEX₈₆ (Figure 2; Table 1). Since phosphatidyl head groups can be degraded faster than the glycosidic head groups, the phospho IPL is considered to be a better reflection of the living microorganisms (Harvey et al., 1986; Schouten et al., 2010), which may explain the deviation between these IPL pools. Relatively rapid conversion of phospho IPL to CL may result in more similar ring distributions and thus TEX₈₆ values between phosphor IPL and CL. Furthermore, in the same study area, the variability in TEX₈₆ was suggested to be due to changes in archaeal community composition in the water column, in which the unusually low TEX₈₆-derived temperature in the coastal SCS was speculated to be linked to MG-II Euryarchaeota (Wang et al., 2015).

Since TEX₈₆ can be influenced by factors other than temperature, the ring index was proposed to evaluate the accuracy of TEX₈₆ in marine sediments (Zhang et al., 2016). Here, CL-, phosphor IPL-, and total IPL-TEX₈₆ values were plotted against RI₁ (Equation 2; Figure 3) using data derived from the SPM and surface sediment samples collected during the current and previous studies (Wei et al., 2011; Ge et al., 2013; Zhang et al., 2013; Wang et al., 2015). Most SPM and surface sediment samples from the open South China Sea plotted within the RI₁-TEX₈₆-confined zone [RI₁ = 3.32 × (TEX₈₆)² − 0.77 × TEX₈₆ + 1.59, ±2σ ~ 0.3; (Zhang et al., 2016)], whereas the majority of samples from coastal SCS and the PR estuary fell above the calibration zone (Figure 3), exhibiting higher ring index values than those from the open SCS. This implies that the GDGT pool in the water column of the PR estuary and coastal SCS comprised relatively more cyclopentane rings than predicted from the measured SWT. If it is assumed that GDGTs produced only by Thaumarchaeota underpin the relationship between SWT and TEX₈₆ (e.g., Schouten et al., 2002), then it follows that either Thaumarchaeota in the PR estuary respond differently to temperature than marine strains, or there is another source of cycloalkyl-containing GDGTs in the estuary.

To further assess the distribution of RI values in the SPM and to explore other contributor(s) to the cyclopentyl GDGT pool in the study area, mean values of CL-, total IPL-, and phospho IPL-ring index were examined (Figure 4). Note that crenarchaeol, a biomarker for MG-I Thaumarchaeota, was excluded from the ring index calculation (RI₂, Equation 3) in order to limit its overwhelming influence on the index. The re-defined RI₂ equation is more sensitive to the variation of cyclopentane-containing GDGTs that might be contributed by other archaea. Compared with the river water and seawater, the highest RI₂ value for either CL (avg. 0.39 ± 0.08) or IPL (avg. 0.48 ± 0.07 for total IPL; avg. 0.47 ± 0.10 for phospho IPL) occurred at station M in the mixing water (Figure 4; Table 1), suggesting that the PR estuary (station M) appeared to be a hot spot of production of GDGTs with cyclopentane moieties. Further confirmation came from the comparison of %GDGT 1–4 in different water settings (Table S1), which showed that the sum of the fractional abundances of the GDGTs with 1–4 cyclopentane moieties in the mixing water was significantly higher than that in the seawater or river water. In the mixing water station, the mean values of total IPL-RI and phospho IPL-RI were not significantly different; both were higher than the CL-RI (Figure 4). In the seawater station, however, the total IPL-RI (0.34 ± 0.07) was more elevated than the phospho IPL-RI (0.22 ± 0.01) and the CL-RI (0.24 ± 0.04; Figure 4). This ring index distribution pattern at station S corresponded to the TEX₈₆-temperature distribution in the SPM and sediment samples (Figure 2), suggesting that cyclopentane-containing GDGTs altered the TEX₈₆ record in the water column, and the vertical transportation of GDGTs from the water column appeared to be a predominant source to the sediment.

Previous studies reported that a significant proportion of MG-II Euryarchaeota was diversely present in the estuarine and coastal regions, including the Pearl River estuary (Liu et al., 2014; Wang et al., 2015), the Yangtze River estuary (Liu M. et al., 2011), and the Jiulong River estuary (Hu et al., 2015). In this study, the MG-II Euryarchaeota 16S rRNA gene averaged 5.4 ± 5.9 × 10⁸ copies L⁻¹ (n = 6) at the mixing water station, which was two to three orders of magnitude higher than that in river water station (avg. 1.5 ± 2.1 × 10⁵ copies L⁻¹, n = 2) and seawater station (avg. 4.9 ± 9.4 × 10⁶ copies L⁻¹, n = 4; Figure 4). Considering the heterotrophic life style of MG-II, which have been demonstrated by the former genetic analysis (Iverson et al., 2012; Li et al., 2015; Martin-Cuadrado et al., 2015) and cultivation experiment (Orsi et al., 2015, 2016), the high abundances of MG-II in the mixing water station seem to be due to the high phototrophs that enhanced by nutrients input from upper river (Gan et al., 2014).

The ratio of MG-II Euryarchaeota 16S rRNA gene abundance to total archaeal 16S rRNA gene abundance ([MG-II 16S]/[Archaea 16S]) in mixing water (avg. 0.25 ± 0.08, n = 6) was significantly higher than in the seawater (avg. 0.10 ± 0.004, n = 4; P < 0.01), whereas it was negligible (avg. <0.0001) in the river water (Figure 4). This observation was also supported by pyrosequencing analysis, which exhibited a linear correlation with the qPCR-based ratio of [MG-II]/[Archaea 16S] (Figure S2; Table S2). These results further confirmed that the PR estuary (mixing zone, salinity avg. 16.6) provided a habitat to sustain a natural enrichment of planktonic MG-II Euryarchaeota.

The presence of (more labile) phospho IPL-GDGTs implied in situ production of isoprenoidal GDGTs in the water column along the entire salinity gradient from the Pearl River to coastal SCS. The elevated phospho IPL-RI in the mixing water additionally suggests that higher relative proportions of GDGTs with 1–4 cyclopentane moieties were produced in the PR estuary. A linear regression analysis confirmed the positive relationship between phospho IPL-derived RI and the [MG-II 16S]/[Archaea 16S] ratio in the water column along the salinity gradient (R² = 0.72, P < 0.01; Figure 5A). Therefore, it is reasonable to hypothesize that MG-II Euryarchaeota preferentially synthesized GDGTs with 1–4 cyclopentane moieties in this region. Moreover, production of GDGTs-1 to -4 by MG-II Euryarchaeota could represent the missing source needed to explain the elevated value of RI and deviation from the RI-TEX₈₆ relationship driven by the preservation of Thaumarchaeota GDGTs in global core-top sediments.

To further constrain the relationship between MG-II Euryarchaeota and cyclopentane-containing GDGTs, linear regression analysis was conducted between the fractional abundance of GDGTs and the [MG-II 16S]/[Archaea 16S] ratio in the SPM along the salinity gradient. Results exhibited a significantly positive linear correlation between the [MG-II 16S]/[Archaea 16S] ratio and the fractional abundance of phospho IPL-based GDGT with cyclopentane moieties (Figures 5B–F). In contrast, we observed no correlation between the [MG-II 16S]/[Archaea 16S] ratio and the fractional abundances of phospho IPL-based GDGT-0 or crenarchaeol (Figures 5G,H; Table 2). Similar trends of the linear correlations were also observed between the [MG-II 16S]/[Archaea 16S] ratio and the CL- and total IPL-GDGTs with 1–4 cyclopentane moieties, with a less significant correlation between the [MG-II 16S]/[Archaea 16S] ratio and the total IPL-based crenarchaeol regioisomer (Table 2). Although MG-II Euryarchaeota were suggested to be an alternative source of crenarchaeol in the ocean (Lincoln et al., 2014a), our study showed an absence of a significant correlation between the distribution of MG-II Euryarchaeota and crenarchaeol (Figure 5H). However, members of MG-II Euryarchaeota living in the North Pacific Subtropical Gyre (i.e., those targeted by Lincoln et al., 2014a) may be different from those living in the coastal zone of the PR estuary. This is consistent with the phylogenetic distribution of MG-II reported by Wang et al. (2015), which showed diverse groups of MG-II living in this region.

By comparison of the R² values and slopes of the regression equations (Figures 5B–F), GDGT-1 exhibited not only the strongest correlation with MG-II Euryarchaeota in the study area, but also largest relative enrichment (i.e., slope of 14.06). If MG-II Euryarchaeota preferentially synthesized GDGT-1, additional contribution of GDGT-1 to the water column of the PR estuary and the coastal SCS would be reflected by an increase in RI values and a substantial decrease in TEX₈₆ values. Offsets in TEX₈₆ values have been similarly proposed to result from the decreased ratio of GDGT-2 to GDGT-3 in the surface sediments of this area (Wang et al., 2015), whereas the increased ratio of GDGT-2 over GDGT-3 in the deep-water column seems to be responsible for a warm bias of TEX₈₆-derived temperature in other marine environments (Taylor et al., 2013; Hernandez-Sanchez et al., 2014). In this study, however, no correlation is exhibited between GDGT-2/3 ratio and [MG-II 16S]/[Archaea 16S] ratio (Figure 5I). A recent study by Kim et al. (2015) suggested that coincident increases in GDGT-2 and the crenarchaeol regioisomer and decreases in GDGT-1 and GDGT-3 shifted TEX₈₆-derived temperatures toward higher values in the deep-water surface sediments of the Mediterranean Sea. In combination with our data (Figures 4, 5), these observations are consistent with the interpretation that planktonic Euryarchaeota have the potential to bias TEX₈₆ by changing the distribution of TEX₈₆-related GDGTs (especially those with 1–3 cyclopentane rings) in different marine environments.

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. The reviewer FS declared a shared affiliation, with no collaboration, with several of the authors to the handling Editor.