The dust used in this study originated from a source of Asian dust events crossing over the NWPO (Shao et al., 2011; Tan et al., 2013), Mu Us Desert (37.92^°N, 107.11^°E). The soil samples were collected and artificially modified following Zhang et al. (2018) to simulate the aging process of dust particles in the atmosphere. Briefly, soil particles smaller than 20 µm were selected and mixed with synthetic cloud water. The treated dust was weighted and stored dry at -20°C.

The 7-day microcosm experiment was performed from April 26 to May 3, 2015, during the spring cruise in NWPO onboard the R/V DongFangHong 2. In the Kuroshio-Oyashio transition region, surface seawater was collected at B1 station (37.93°N, 146.98°E) from a depth of 5 m, using the water sampler mounted on CTD (SeaBird SBE 911). For the ambient baseline condition (hereafter ‘day 0’), after prefiltered through 200-μm mesh to remove large zooplankton, seawater was immediately collected for size-fractionated Chl a, DNA and ancillary measurements as described below. For the microcosm experiment, prefiltered seawater was aliquoted into nine 20-L acid-washed polycarbonate carboys (round and clear plastic container for large volume culture, Nalgene). Three groups were prepared in triplicate: i) ‘control’, no dust addition; ii) ‘low’, 0.2 mg L^-1 of dust particles in incubated seawater; iii) ‘high’, 2 mg L^-1 of dust. We added such amount of dust, aiming to simulate the deposition flux to the surface of the NWPO in mild (2 g m^-2) and strong (20 g m^-2) dust events, respectively (Iwamoto et al., 2011; Shi et al., 2012). The carboys were placed in an on-deck incubator with circulating seawater to keep near-ambient surface temperature. Neutral density filters were used to cover the carboys, reducing light intensity to ~40% of the incident. Temperature and light intensity were monitored three times per day. Size-fractionated Chl a and nutrient concentrations were measured on daily basis, while seawater samples for morphological and metabarcoding analysis were collected on day 0, 2, 4 and 7.

For each carboy, 300 mL incubated seawater was sequentially filtered through 3 µm (Φ47 mm, Millipore) and 0.2 µm polycarbonate filters. Filters were wrapped in aluminum foil and stored at -20 °C. After 24-h extraction by 90% acetone (v/v) at -20°C in dark, the pigment content on each filter were quantified by measuring the absorbance of supernatant using a Trilogy fluorometer (Turner Designs). The total Chl a concentration was the sum of the two size fractions.

The leachable inorganic nutrients in the artificially modified dust were determined in the laboratory before the cruise. Briefly, the dust (10 mg L^-1) was added into prefiltered (0.45 μm) oligotrophic seawater, and placed in ultrasonic bath at 0 °C for 30 min. The leaching solution was then filtered through acid-washed cellulose acetate membranes and the filtrate was immediately stored in acid-washed high-density polyethylene bottles at -20°C for nutrient analysis. In the field, 200 mL ambient seawater and 200 mL incubated seawater from each carboy were collected following the same sampling protocol. NO₃ ^-, NO₂ ^-, PO₄ ^3- and SiO₃ ^2- concentrations were measured by Auto Analyzer3 continuous-flow analyser (SEAL Analytical).

Subsample (500 mL) from the ambient seawater (day 0) and incubated seawater (day 2, 4 and 7) were fixed with 1% neutral Utermöhl’s solutions and stored at 4 °C in dark. Before microscopic examinations, fixed samples were settled for 48 h in dark and then concentrated to 10~50 mL. Phytoplankton were identified under microscope (YS100, Nikon) and cell count was obtained using Sedgwick-Rafter chamber.

On day 0, day 2, day 4 and day 7, 1 L sample water were sequentially filtered through 3 µm and 0.2 µm polycarbonate membrane (Millipore, Φ47 mm). The filters were immediately stored in 1 mL 100 mM EDTA, and then frozen in liquid nitrogen until DNA extraction.

Membrane filters of 3 μm and 0.2 μm from the same incubated samples were pooled together for DNA extraction. DNA extraction followed the method described by Yuan et al. (2015). Briefly, membrane filters were cut with sterile scissors and homogenized using microbeads in sample homogenizer (Bioprep-24, Allsheng). Homogenate in 2 mL DNA extraction buffer (100 mM EDTA, 1% SDS, and 10 mg mL^-1 proteinase K) was incubated at 56°C for 24 h. DNA was extracted using the CTAB method, purified with the DNA Clean & Concentrator kit (Zymo Research), and eluted in 25 μL 10 mM Tris-HCl. DNA concentrations and purities were examined on a NanoDrop 2000 Spectrophotometer (Thermo Scientific). DNA integrity was determined using 1% agarose gel electrophoresis.

The V4 region of the eukaryotic 18S rRNA gene was amplified with the primer pair 18SNCF3 (5’- CGCGHAAATTRCCCAATCY -3’) and 18SV4R2 (5’- CTWACTTTCGTTCTTGATYAA -3’). Thermal cycling was as follows: 95°C for 5 min, followed by 5 cycles of 95°C for 10 s, 56°C for 30 s, and 72°C for 30 s, and then 25 cycles of 95°C for 10 s, 52°C for 30 s, and 72°C for 30 s with a final extension at 72°C for 10 min. The V4 region of prokaryotic 16S rRNA gene was amplified with primer pair 515F (5’-GTGYCAGCMGCCGCGGTAA-3’) and 806R (5’- CAAYAAATCYRAGAATTTCACCTCT-3’) (Caporaso et al., 2011), following the thermal cycling: 95°C for 5 min, followed by 25 cycles of 95°C for 30 s, 56°C for 30 s, and 72°C for 40 s, with a final extension at 72°C for 10 min. PCR reactions were performed in a total volume of 25 µL using Ex Taq (Takara), with 20~300 ng template gDNA. All amplicons were purified and paired-end sequenced using Illumina MiSeq PE300 platform (GENEWIZ).

Raw reads were subjected to quality filtration by removing reads of low quality, adaptor contamination, high N content and low complexity to obtain clean reads with Cutadapt (v1.9.1; Martin, 2011), Vsearch (v1.9.6; Rognes et al., 2016) and Qiime (v1.9.1; Bolyen et al., 2019). The remaining high-quality paired-end reads were assembled into metabarcodes, and non-singleton ones subsequently clustered into OTUs (Operational Taxonomy Units) with 97% similarity using USEARCH (v7.0.1090; Edgar, 2013). For taxonomic assignment, the representative OTU sequences of 18S rDNA were classified using sintax algorithm (Edgar, 2016) based on PR2 database (Guillou et al., 2012), with a confidence threshold set at 0.8. The representative OTU sequences of 16S rDNA were classified using the RDP classifier (Ribosomal Database Program, v2.2; Wang et al., 2007) based on the Greengene database (DeSantis et al., 2006). The ecological functions of bacterial communities were further annotated using FAPROTAX (Louca et al., 2016).

In addition to the total community, abundant, moderate and rare subcommunities were analyzed separately. OTUs with a relative abundance greater than 0.1% of the overall sequences in a single library (i.e. sample) were considered as ‘abundant’ OTUs, those with a relative abundance less than 0.01% were considered as ‘rare’ OTUs, and the rest were considered as ‘moderate’ OTUs (Campbell et al., 2011). Moreover, we further classified all OTUs into six categories (Xue et al., 2018): Define OTUs that were abundant OTU in all samples as ‘always abundant’ OTUs; OTUs that were rare OTU in all samples as ‘always rare’ OTUs; OTUs that were always moderate OTU in all samples as ‘always moderate’ OTUs; OTUs that were abundant or moderate OTU in all samples as ‘conditional abundant’ OTUs; OTUs that were rare or moderate OTU in all samples as ‘conditional rare’ OTUs; and the rest were considered as ‘conditional abundant and rare’ OTUs.

Net nitrogen conversion index (NCEI, μg μmol^-1) was calculated using the equation below (Zhang et al., 2018):

where t represented the incubation duration (7 days), Chl a_i represented the concentration of Chl a at the ith day (μg L^-1), Chl a₀ represented the initial concentration of Chl a on the 0^th day in the control (μg L^-1), ΔN represented consumption value of dissolved inorganic nitrogen (DIN) during incubation (μmol L^-1).

For alpha-diversity, the ‘vegan’ package in R (v4.2.2) was used to calculate diversity indices including Chao 1, Pielou’ evenness and Shannon-Wiener index (Oksanen et al., 2013). Beta-diversity was visualized in Non-Metric Multidimensional Scaling (NMDS) with Bray-Curtis dissimilarity distance calculated in PRIMER 6 (Clarke and Gorley, 2009), and the significance of differences among groups was tested using ANOSIM. Beta diversity was partitioned into turnover and nestedness components by applying the function beta.pair within the R package (‘betapart’) following Jiao et al. (2017).

To characterize the differential response of different taxa, the fold change of relative abundance between treatment and control were calculated. Heatmap with hierarchical clustering was generated by ‘heatmap’ and ‘ggplot’ package in R. Complete linkage and Euclidean distance were used in hierarchical clustering. To further grouping the time series pattern, soft clustering was performed in Mfuzz.

To ensure that the network analysis was supported by sufficient data, control and low dust treatment were merged into ‘Con-Low’ group to characterize planktonic microbiome that were not or weakly affected by dust additions. Low and high dust treatments were merged into the ‘Low-High’ group to characterize the community affected or strongly affected by dust additions. Spearman’s correlations between the relative abundance of OTUs were calculated using the ‘rcorr’ function in the ‘Hmisc’ package, and co-occurrence networks were constructed based on these correlations. OTUs detected in over 30% of all samples with robust correlation (|spearman rho| > 0.8 and p < 0.01) were selected for network analysis. Networks were visualized in Gephi (v0.9.4; Bastian et al., 2009) and topological parameters of networks were calculated in the ‘igraph’ R package (Csardi and Nepusz, 2006). Spearman’s correlation between topological parameters and environmental factors were calculated and visualized in circus plot using the ‘circlize’ package in R. Nodes featured with high degree (>120) and low betweenness centrality (<2000) were considered to be keystone species in the co-occurrence network.

The baseline conditions at B1 were featured as HNLC, with a N:P ratio (14:1) close to the Redfield ratio but low Chl a concentration 0.78 g L^-1 ( Supplementary Table 1 ). The >3 μm size fraction only accounted for 26.1% of the total Chl a.

Prior to the field experiments, we determined the leachable nutrient contents in the dust. The NO₃ ^-+NO₂ content reached 989.23 μmol g^-1, while the PO₄ ^3- (4.53 μmol g^-1) and SiO₃ ^2- (6.60 μmol g^-1) contents were relatively low ( Figure 1A ). By summing the baseline values and dust-derived ones, we calculated the amended nutrient concentrations in the dust treatments, yielding 3.9% and 39.0% increase of NO₃ ^-+NO₂ in the low and high dust treatments, respectively ( Figure 1B ). But the increases of PO₄ ^3- and SiO₃ ^2- concentrations were negligible, even in the high dust treatments, which were only 2.5% and 0.2%, respectively.

During the on-board incubation, nutrients and size-fractionated Chl a were monitored on a daily basis. Overall, nutrients showed a decline trend over the time in all groups ( Figures 2A–C ), with the greatest decrease found in the high dust treatment (13.6~43.0%). In spite of the decrease, the NO₃ ^-+NO₂ concentrations in the high dust treatment remained higher than those in control and low treatment, reaching a comparable concentration on day 7. The PO₄ ^3- and SiO₃ ^2- concentrations were roughly similar among groups in the first 4 days, but the differences were more pronounced towards the end of the incubation. N:P ratios in the control and low dust treatment varied in the range of 11-17, while those in the high treatment (16-18) maintained above Redfield ratio with a slight decrease detected on day 6 and 7. Although the N:Si ratios were elevated by dust-derived N, ratios remained below one regardless of treatments (Figure 2E).

The consumption of nutrients came with the increase of Chl a concentration and phytoplankton abundance ( Figures 2G–I ). The high dust addition promoted the total Chl a concentration throughout the incubation with the significant increase detected after day 4 (p<0.05), and the value peaked on day 7 being 293.6% of the control. But the enhancing effect of low dust addition on total Chl a was not significant until day 7 (p>0.05). The Chl a increase was mainly contributed by >3 µm size fraction, taking up 69.7% and 84.0% of the total Chl a in low and high dust treatment on day 7, respectively. The specific growth rate of larger phytoplankton (using either >3 µm size fraction Chl a or phytoplankton abundance) was positively correlated with available N in the incubation (R^{2 =} 0.2134, p<0.05, Supplementary Figure 3 ). However, such correlation was not observed in 0.2~3 µm size fraction. NCEI value further indicated the positive effect of the added N on phytoplankton growth, which was more pronounced in high dust treatment and >3 µm size fraction ( Supplementary Figure 4 ).

Microeukaryotic and prokaryotic plankton communities were examined parallelly, yielding 903,119 high-quality 18S rDNA metabarcodes and 2,510 OTUs, 464,927 high-quality 16S rDNA metabarcodes and 838 OTUs, respectively. The taxonomic composition was uneven across samples showing both temporal and intergroup variation ( Figure 3 ). The microeukaryotic communities were further analyzed in three categories, phytoplankton, heterotrophic protist and metazoan. the relative abundance of phytoplankton reads in dust treatments peaked on day 7, which was 2.05-2.16 folds of the control ( Supplementary Figure 5A ), in concert with the trend of Chl a concentration and cell density ( Supplementary Figure 7 ). Compared to the baseline condition, a shift from dinoflagellate dominance toward diatom dominance was observed in all groups, but dust addition particularly enhanced the relative abundance of diatom, which was increased by 120.5~137.7% relative to the control on day 7. Such increase was mainly contributed by genera Thalassiosira, Chaetoceros, Pseudo-nitzschia, and Navicula ( Supplementary Figure 8 ), accounting for 80.28~87.88% of the total diatom reads in dust treatments. The relative abundance of reads and cell density were positively correlated in diatom (R² = 0.6348, p<0.01) but not observed in dinoflagellate ( Supplementary Figure 7 ).

The relative abundance of heterotrophic protist in the dust treatments generally showed an increasing trend over time, while that in the control peaked on day 4 and fell back to a comparable level to the baseline ( Figure 3B ; Supplementary Figure 5B ). The relative abundance of Cercozoa in the low and high dust treatments was 2.34 and 1.71 folds of that in the control, respectively. This was mainly due to the increase of Cercozoa, in which filose taxa (subphylum Filosa) were predominant, taking up 29.2~39.7% of the total Cercozoa reads in dust treatments. The relative abundance of Cercozoa was positively correlated with diatom (R² = 0.7366, p<0.01, Supplementary Figure 9A ). Moreover, despite lower abundance, Labyrinthulomycetes was also stimulated by dust addition, being 148.14~277.78% of the control on day 7. Majority (>93.1%) of the metazoan reads were affiliated with copepod, followed by Urochodata (tunicate) and Hydrozoa. Slight increases were observed in all groups on day 2 and 4, but pronounced decrease occurred in dust treatments on day 7 with the relative abundance falling to 62.5~71.5% of that in the control. The variation in taxonomic composition of bacterioplankton was not as strong as that of microeukaryotes ( Figure 3D ). Gammaproteobacteria, Alphaproteobacteria and Flavobacteriia dominated across all samples (60.09~80.25%). The most remarkable change was observed in Saprospirae, of which the relative abundance on day 7 increased 5.2-8.3 folds relative to the baseline condition, and that in the low and high dust treatments was 156.74% and 144.01% of that in the control. The relative abundance of Saprospirae was positively correlated with that of diatoms (R^{2 =} 0.9066, p<0.001, Supplementary Figure 9B ). The ecological function cellulolysis detected in bacterioplankton was exclusively affiliated with Saprospirae ( Supplementary Figure 10 ).

To better characterize the differential response among taxa to dust addition and describe the temporal dynamics, we performed cluster analyses based on the fold change between the relative abundance of certain taxa in the dust treatments and that in the control. With the effect of dust, the variability of taxonomic profile in eukaryotes was significantly greater than that in prokaryotes, but functional profile was relatively stable in terms of prokaryotes, except for the function cellulolysis ( Figure 3E ; Supplementary Figure 11 ). Both lineage- and OTU-level clustering yielded 5 clusters ( Figures 3E, F ). Generally, cluster 1-3 pointed to the promotional effects of dust which was characterized by the higher relative abundance in the dust treatments on day 7, while cluster 4-5 showed the opposite trend. When considering temporal variations, more complex response was observed among lineages. Except for the metazoan, the majority of eukaryotic taxa fell into cluster 1-3. The dominant prokaryotic taxa, Gammaproteobacteria, Flavobacteriia and Saprospirae, were also found in cluster 1-3, while the rest of prokaryotes in cluster 4-5.

To investigate the dynamics of different subcommunities throughout the incubation, we categorized OTUs by their relative abundance in each sample. Both microeukaryotic and prokaryotic plankton communities were dominated by few abundant taxa, as 8.0~24.7% of OTUs per sample were classified as abundant taxa and accounted for 49.3~87.8% of reads ( Figure 4A ). Non-abundant taxa (i.e. moderate and rare taxa) contributed to the majority of community richness, representing 75.3~92.0% of OTUs per sample. Remarkably, the contribution of moderate taxa was significantly higher in prokaryotic community relative to that in eukaryotic one, taking up averagely 35.67% of the richness and 41.32% of the reads ( Supplementary Figure 12 ). Moreover, most of the control- or dust treatment-specific OTUs belonged to rare taxa ( Supplementary Figure 13 ).

To address the possible transitions among abundant, moderate and rare taxa, we further categorized OTUs by their relative abundance across all samples (see methods). A total of 617 18S rDNA OTUs (59.91%) representing 62.52% of reads, and 596 16S rDNA OTUs (70.87%) representing 70.85% of reads, were stable in terms of their relative abundance level throughout the incubation ( Supplementary Table 2 ). The rest OTUs were assigned as oscillating taxa (conditionally to be rare and/or abundant taxa) indicating 29.13~37.48% of the community had been undergone transition. Further quantification revealed that all three subcommunities were dynamic with the moderate taxa having the highest proportion of OTUs (55.71-62.26%) that transitioned to abundant or rare taxa ( Figure 4B ). Moreover, dust addition stimulated the transitions of rare taxa. Particularly, the proportion of eukaryotic rare taxa that transitioned to abundant taxa increased with dust addition. Similar trend was also observed in prokaryotic rare taxa that transitioned to moderate taxa.

Temporal and intergroup variation were also observed in alpha diversity but eukaryotic and prokaryotic communities showed different patterns. For eukaryotic communities, dust addition decreased alpha diversity from day 0 to day 4, followed by a pronounced increase of Pielou’s evenness and Shannon diversity towards day 7 ( Figures 5C–E ). For prokaryotic communities, dust addition promoted the alpha diversity (Chao 1, Pielou’s evenness and Shannon diversity) and functional richness from day 0 to day 4, followed by a decrease towards day 7 ( Figures 5F–H ). Although relative to the control, higher functional evenness and diversity with lower functional redundancy was observed in high dust treatment on day 4, the general trend over time was similar among control and dust treatments ( Supplementary Figure 12 ).

NMDS ordinations and ANOSIM analysis revealed the distinct separation of plankton communities by time ( Figures 5A, B ). On the same sampling day, the community similarity between dust treatments and control decreased ( Figures 5A, B ), indicating that dust also played a role in community divergence, which was greater in microeukaryotic plankton relative to that in the prokaryotes. This is also the case when analyzing the functional structure of prokaryotic community with much weaker community divergence (similarity >80%, Supplementary Figure 4A ). We further partitioned total beta diversity into turnover and nestedness, revealing that species replacement (turnover) played a major role in driving the community succession, and the percentage of turnover in eukaryotic communities was significantly higher than that in prokaryotic ones ( Supplementary Figure 15 ). All these results indicated that the sensitivity of microeukaryotic and prokaryotic plankton to dust addition was different.

Correlation analysis suggested that the relative abundance of some lineages was significantly correlated with inorganic nutrient conditions ( Supplementary Figure 16 ). Remarkably, the relative abundance of most eukaryotic lineages (18 out of 20) was significantly correlated with prokaryotic community, while 7 out of 20 prokaryotic lineages showed significantly correlation with eukaryotes. This pointed to that microeukaryotic plankton was more sensitive to cross-domain biotic factors, which was also revealed by Mantel test at community level ( Supplementary Table 3 ).

To further explore the potential biotic interaction in the planktonic microbiome affected by dust addition, we conducted a co-occurrence network analysis ( Figure 6 ). Compared to the Con-Low network, the Low-High network was more complex with 4% more nodes, 13% more positive edges and 30% more negative edges ( Supplementary Table 5 ). Although positive edges dominated both networks (>60%), the proportions of negative edges slightly increased (30%) due to dust addition. Among the increased edges, 8.4%, 75.2% and 16.4% were contributed by the increased eukaryote-prokaryote, eukaryote-eukaryote and prokaryote-prokaryote connections, respectively. Phytoplankton-bacteria edges dominated the eukaryote-prokaryote connections (>45.06%), and the number of connections between diatoms and their correlated taxa including Cercozoa, Labyrinthulomycetes, Copepoda and Saprospirae were higher in Low-High network ( Supplementary Table 6 ). We compared the unique node-level topological characteristics of the two networks and found that the degree and closeness centrality of the Low-High network were significantly higher than that of the Con-Low network (p<0.05), while betweenness centrality did not show significant difference ( Figures 6C–E ). Our results implied that microbial communities display more complex networks when exposed to increased dust addition treatments. A total of 44 OTUs belonging to 14 taxa in the Con-Low network and 87 OTUs belonging to 22 taxa in the Low-High network were identified as keystone taxa ( Figure 6F ). Only 5 keystone taxa were shared between two networks, belonging to Alphaproteobacteria (2 OTUs), Gammaproteobacteria (2 OTUs) and Syndiniales (1 OTU). Five OTUs belonging to diatoms were identified as keystone taxa.

The N:P ratio at B1 station was close to Redfield ratio, indicating that N was not the primary limiting factor ( Figure 2D ). Despite this, dust addition still stimulated the growth of phytoplankton by supplying a considerable amount of N and the promotion effect showed a dust loading-dependent pattern as evidenced by both microscopic analysis and Chl a concentration ( Figures 2F–I ). Our results are in concert with previous studies in Kuroshio-Oyashio transition region (Zhang et al., 2018, 2019a; b; Zhang et al., 2020). Dust-derived N has been recognized as a supplementation of the N demand for phytoplankton in this HNLC zone, when compared with the response of phytoplankton in LNLC and eutrophic zones (Zhang et al., 2019). We further used NCEI to quantify the N-utilization efficiency among groups ( Supplementary Figure 4 ). A higher NCEI value indicates that the positive effect of added N on phytoplankton growth is stronger. Although Chl a increased with dust loading, the NCEI value is supposed to be independent of N loading. However, we found that the NCEI value in high dust treatment was higher than that in low dust treatment, implying that other substances leached from dust may influence the N-utilization efficiency in phytoplankton. Phytoplankton at B1 is possibly limited by Fe and other trace metals, such as Cu and Zn (Yang et al., 2019), which could be alleviated by dust deposition and in turn facilitates the N uptake as shown in high dust treatment. Moreover, the Asian dust provides considerable amount of N but negligible P, making the N:P ratios in high dust treatment (16-18) maintained above Redfield ratio, which potentially caused depletion of dissolved inorganic phosphorus (DIP). In our results, the enhanced DIN uptake by phytoplankton suggested the utilization of a supplementary P source, such as dissolved organic phosphorus (DOP), which is more abundant relative to DIP in the oceans. A recent study by Jin et al. (2024) provides a possible explanation: atmospheric deposition can stimulate the biological utilization of DOP by depleting DIP due to increasing N. Under DIP deficiency, phytoplankton most commonly use alkaline phosphatase but also other hydrolases to convert DOP to DIP for cellular uptake (Lin et al., 2016) and the activity of alkaline phosphatase is induced by low DIP concentrations. Based on microcosm experiments and global ocean analysis, Jin et al. (2024) proposed a threshold ratio of Chl a to DIP [Log₁₀ (Chl a/DIP) > 1.20] to determine whether DOP is an important P source. The value of [Log₁₀ (Chl a/DIP)] at B1 is 2.31, pointing to that DOP utilization might play a role in enhancing DIN uptake by phytoplankton.

Although dust input induced significant phytoplankton biomass accumulation, the degree of impact was varied between size fractions. Size as one of the most important fitness-related traits, its distribution can be used as indicator of biological transitions, and even regime shift in ecosystem (Spanbauer et al., 2016). In response to dust input, phytoplankton size structure shifted towards larger cells over incubation time ( Figures 2G–I ). Correlation between size-fractioned Chl a and N supply indicates that dust addition favored the growth of >3 μm-sized phytoplankton, which can be partially explained by the higher nutrient uptake rate of larger phytoplankton under nutrient-replete conditions (Maranon, 2015). Although dust addition slightly promoted the growth of picoplankton, no time-dependent increase was observed ( Figure 2H ). Trophic interaction may provide a possible explanation: increased abundance of protist grazers and their efficient grazing on picophytoplankton under atmospheric deposition keep the picophytoplankton biomass under control, as found in the microcosm experiments carried out in South China Sea (Guo et al., 2012; 2014). Contrastingly, a previous study in the Southern Yellow Sea showed that the abundance of heterotrophic and mixotrophic protists were inhibited by Asian dust deposition (Chen et al., 2019). Thus, trophic structure should be considered when evaluating the effects of atmospheric deposition on marine planktonic ecosystems. In addition, usually mesozooplankton (>200 μm) is excluded in microcosm experiments, which may technically lower the grazing pressure for micro- and nanophytoplankton.

Besides size-dependent response of phytoplankton to dust deposition, we further profiled the differential response at a fine taxonomic resolution covering both microeukaryotic and prokaryotic plankton. Communities were distinctly separated by time, and dust exerted stronger effect on community divergence towards the end of incubation ( Figures 5A, B ), indicating that both magnitude and duration of the deposition event played a role in modulating planktonic microbiome. Numerous incubation experiments in the NWPO find that the Chl a generally started to increase with 0-5 days after dust deposition (e.g. Zhang et al., 2019a, 2020), and subarctic NWPO showed the longest response time (up to 40 days) to high dust deposition events, compared to other regions (Meng et al., 2022). We have detected the change of community composition and structure since day 2 and characterized the temporal dynamics by clustering analysis ( Figures 3E, F ). Time- and lineage-dependent responses were complex within 7-day incubation.

Mantel test revealed that biotic factor shaped both bacterial and microeukaryotic community more than abiotic factors, particularly the cross-domain interaction significantly affected microeukaryotic community ( Supplementary Table 3 ). Dust deposition, as a disturbance to the surface ocean, may first affect the most sensitive lineages and its effects could pass on other organisms through interspecies interaction, and other components of the ecosystems (Battisti et al., 2016). Here, as the first attempt, we used co-occurrence network to portrait the biotic interactions responding to dust deposition. Ecological complexities can be depicted through networks, where species are represented as nodes and their interactions as connecting edges. This representation is essential for understanding the dynamics of species interactions within an ecosystem (Berry and Widder, 2014). Dust addition altered the network topology by increasing the number of edges, closeness centrality and the number of keystone taxa, resulting in a more densely linked network ( Figure 6 ). Network complexity and connectivity are considered to be associated with the sensitivity of microbes to environmental disturbances, and strongly correlated with stability (Zhou et al., 2010). Some insights have been revealed into the impact of environmental stressors, such as climate warming and salinity stress, on the microbial ecological network (Yuan et al., 2021; Li et al., 2023). But the induced changes in complexity and stability of network are different, depending on the type of stressor and temporal scale. Remarkably, our co-occurrence network analysis revealed that the number of diatom-related edges were increased by dust addition, suggesting that diatoms were actively interconnected with their symbionts, competitors, predators and other associated microbes. This confirms that dust deposition does not only induce a succession of phytoplankton in NWPO, which has drawn most of the attention, but also a progression of heterotrophic protists, metazoan, bacteria and other players in the planktonic ecosystem.

Although positive association dominated both networks, more negative associations were induced by dust addition ( Figure 6 ; Supplementary Table 4 ). Positive and negative associations primarily represent cooperation and competition, respectively (Yuan et al., 2021). Dai et al. (2022) studied the effects of nutrient supply on network association in marine bacterial communities and found that positive associations are more dominant in nutrient-scarce pelagic regions, while the proportion of negative associations are higher in nutrient-rich coastal sediment. High nutrient concentrations favor the actively growing bacteria which may inhibit the growth of competitors by producing metabolites (Ratzke et al., 2020), or cause competition for other essential nutrients. Based on this, we speculate that the slight increase of negative associations in dust-affected network could be partly explained by the nutrient supply from dust, and diatom is one of the strong competitors. Previous findings on negative association induced by nutrient supply are solely based on soil and marine bacterial studies (Ratzke et al., 2020; Dai et al., 2022), but whether this can be applied to interactions among microeukaryotes and cross-domain interactions remain elusive.

Dust-affected network also had greater number and more diverse keystone taxa. Keystone taxa regulate the network stability, which could be significantly reduced when keystone taxa were removed (Liu et al., 2022). Additionally, only 5 keystone taxa were shared between two networks, and several OTUs belonging to diatoms and their associated microbes became keystone taxa in the dust-affected network. This implies that dust addition may alter the microbial network by displacing keystone nodes, as they are associated with topological feature in co-occurrence network (Berry and Widder, 2014).