Climate change in the Indo-Paciﬁc basin from mid- to late Holocene

Changes in climate mean state profoundly impact climate variability. Here, we quantify slow changes in the mean climate induced by the variations in the Earth’s orbit from mid- to late Holocene, and their feedback on the main modes of climate variability. We focus on the Indo-Paciﬁc system and show that mid-Holocene conditions favored the dominance of an equatorial dipole mode in the Indian Ocean (IO), independent of the El Niño Southern Oscillation (ENSO) and diﬀerent from the IO Dipole (IOD) observed today. Mean state changes induced a gradual shift to an IO basin mode that along with the IOD modulates most of the IO variance at present. The climate modes evolution and their connectivity changes are investigated over 6,000 years using a complex network methodology and principal component analysis. To characterize the nature of this transition, we explore the trajectory of the Indo-Paciﬁc climate by accounting for its spatiotemporal and multivariable dependency. The full trajectory of the system is explored from a dynamical system perspective by con-structing a state space representation. The manifold embedded in the 10 4 -dimensional state space provides a compact representation of the system evolution and points to a gradual shift of the basin of attractions in the tropics. This approach, together with a mean state analysis, reveals that a strengthening of the Walker circulation set the stage for a shift in modes in both basins.

Abstract Changes in climate mean state profoundly impact climate variability. Here, we quantify slow changes in the mean climate induced by the variations in the Earth's orbit from mid-to late Holocene, and their feedback on the main modes of climate variability. We focus on the Indo-Pacific system and show that mid-Holocene conditions favored the dominance of an equatorial dipole mode in the Indian Ocean (IO), independent of the El Niño Southern Oscillation (ENSO) and different from the IO Dipole (IOD) observed today. Mean state changes induced a gradual shift to an IO basin mode that along with the IOD modulates most of the IO variance at present. The climate modes evolution and their connectivity changes are investigated over 6,000 years using a complex network methodology and principal component analysis. To characterize the nature of this transition, we explore the trajectory of the Indo-Pacific climate by accounting for its spatiotemporal and multivariable dependency. The full trajectory of the system is explored from a dynamical system perspective by constructing a state space representation. The manifold embedded in the 10 4 -dimensional state space provides a compact representation of the system evolution and points to a gradual shift of the basin of attractions in the tropics. This approach, together with a mean state analysis, reveals that a strengthening of the Walker circulation set the stage for a shift in modes in both basins.
Keywords Indo-Pacific system · mid-to late Holocene · climate regime shifts · climate mean state · climate variability · complex networks · dynamical systems

Introduction
Tropical climate variability comprises local physical processes responsible for nonlocal, large-scale responses in climate fields, from rainfall to surface temperatures, with ample social and economic repercussions. The representation of these processes and their multiscale interactions in state-of-the-art global climate models (GCMs) remains challenging. Changes in solar radiation, albedo, volcanic activity and greenhouse gas concentrations can alter both the climate mean state, its modes of variability and their interactions, but disentangling causes, effects and interrelationships is no easy task.
Present-day interannual variability in the tropical Pacific is dominated by the El Niño Southern Oscillation (ENSO) (Dĳkstra, 2005), with global repercussions (Brönnimann et al., 2004;Yu et al., 2012;Ludescher et al., 2013;He et al., 2020), while the tropical Indian Ocean (IO) exhibits two main modes, the IOB (Indian Ocean Basin mode; (Klein et al., 1999)), and the IOD (Indian Ocean Dipole; (Webster et al., 1999;Saji et al., 1999)). Between 1980 and 2019, the IOB explained about 30% of the IO variance, and the IOD was responsible for another 10%. The IOB develops during boreal winter, peaks in boreal spring, decays in summer (Schott et al., 2009) and is understood to be a response to ENSO teleconnections that promote basin-wide Sea Surface Temperature (SST) anomalies in the IO through an atmospheric bridge and oceanic adjustment (Lau and Nath, 2003). A positive IOD, on the other hand, is characterized by cool SST anomalies and a shallow thermocline in the eastern IO, offshore of Java and Sumatra, and warm SST anomalies and a deep thermocline in the western IO; the opposite holds for its negative phase. The IOD initiates in boreal spring and peaks in September-November (Schott et al., 2009) and its impacts extend to Australia, east Africa, Japan and Europe (Ratnam et al., 2020).
ENSO, the IOD and the Asian monsoon underwent profound changes during the Holocene (e.g., Cobb et al. (2003); Lu et al. (2018); Abram et al. (2020)) according to paleo records. Mid-Holocene ENSO had weaker variance (Moy et al., 2002;Emile-Geay et al., 2016;Thompson et al., 2017;Grothe et al., 2020), likely in response to a different insolation (Abram et al., 2007;An et al., 2018). At the same time, paleo corals off Sumatra recorded long and intense upwelling events, rather independent from ENSO (Abram et al., 2007(Abram et al., , 2020. Finally, in the last 6000 years the Indo-Asian monsoon system changed significantly, with Holocene transient simulations suggesting strong drying trend and increasing variability over time (Braconnot et al., 2019a;Crétat et al., 2020). Further back in time, paleo reconstructions and model simulations of the Last Glacial Maximum (LGM) support the existence of an Equatorial Dipole Mode (EDM) modulating the IO in past climates (Thirumalai et al., 2019;DiNezio et al., 2020). The EDM was characterized by SST anomalies of opposite sign at the east and west side of the basin, similarly to the IOD, but shifted northward, towards the equator, compared to the modern analogue.
Here, we explore the transient climate of the Indo-Pacific basin over 6,000 years from the mid-Holocene to 1950 with in a simulation performed with the IPSL Earth system model (Braconnot et al., 2019b). Using complex networks, principal component analysis, and tools borrowed from dynamical system theory, we characterize the evolution of the Indo-Pacific climate, and uncover the causal link between these changes and slow variations in the Earth's orbit. This investigation uncovers a slow, non-abrupt shift from an Indo-Pacific system modulated primarily and independently by the EDM in the IO and a weaker-than-today ENSO in the Pacific, to the current one in which the variability is dominated by the strongly connected IOB -ENSO modes.
This paper is organized as follows: section 2 presents the climate simulation analyzed. Section 3 briefly describes the complex network framework and introduces a new framework for dimensionality reduction in climate science by means of state space embedding. Results are presented in sections 4, 5, 6. Section 7 concludes the work.

Model Simulation
We analyze a 6000-years long transient simulation realized with the IPSL model (Braconnot et al., 2019a,b). The ocean component has 2 • horizontal resolution with spatial refinement at the equator and in the Arctic, and 31 vertical levels, while the atmospheric module has a resolution of 2 • in longitude and 1.12 • in latitude, with 39 levels. The configuration includes the 11-layer ORCHIDEE hydrology model and a dynamical vegetation module (Braconnot et al., 2019b). After a spin-up that produces initial conditions in equilibrium with the external forcing, the orbital parameters and the atmospheric chemical composition (Berger, 1978;Otto-Bliesner et al., 2017) are updated yearly.
In this paper, our focus is on the Indo-Pacific region in the latitudinal and longitudinal ranges [30S -30N] and [40E -70W] and we consider surface temperature (ST), velocity potential at 200hPa (VP ), zonal and meridional velocity (U and V) at 850 hPa and precipitation (Precip).

External forcing
The external forcing from 6000 BP to 1950 includes two contributions: changes in atmospheric trace gases and the evolution of the Earth's orbit. Changes in the Earth's orbit -that refer to changes in eccentricity, obliquity, and longitude of the perihelion -resulted in a shift of ∼90 degrees in the position of the aphelion and vernal equinox (Joussaume and Braconnot, 1997) over the last 6000 years.
This shift in the orbital configuration impacted the seasonality and spatial distribution of incoming solar radiation at the top of the atmosphere (TOA) that decreased by ∼ 5% (Zhao et al., 2005) from mid-to late Holocene during boreal summer in the Northern Hemisphere (NH), and reversed for austral summer in the Southern Hemisphere (SH). The end result was a weakening (strengthening) of the inter-hemispheric energy gradient during boreal (austral) summer from mid-to late Holocene.
In Figure 1a-b, we show the evolution of the annual mean solar insolation averaged in both hemispheres along with the contributions in boreal (June to September) and austral (December to March) summer. Major differences are in summer, with nearly opposite changes across hemispheres (see Figure 1b). Changes in radiative forcing driven by trace gases were small through most of the Holocene (Joos and Spahni, 2008;Braconnot et al., 2019b), with a large, exponential drift since the industrial revolution ( Figure 1c).

Methodology
The framework adopted in this study allows for a comprehensive analysis of climate fields by considering their spatiotemporal and multi-variable dependence. It builds on concepts of dynamical system theory and machine learning. Here, we introduce a complex network framework useful to track changes in modes of variability and their connections and then we present a dynamical system-based analysis relying on concepts of state space embedding (Gibson et al., 2008;Cvitanović et al., 2016) to visualize and study transient climate change at spatiotemporal scales and across multiple fields simultaneously.

Complex network framework
The complex network framework was recently proposed in Falasca et al. (2020) as a dimensionality reduction scheme that is useful for analyzing long paleoclimate integrations. Given a generic spatiotemporal field Y ∈ R T,m it reduces its dimensionality Given the horizontal velocity vector V = ( , ) at 200hPa, the velocity potential is such that V = −∇ . is then a scalar field describing the divergent irrotational part of V. from a large number of grid cells into spatially contiguous partitions (hereafter domains), with << . Each domain is a set of time series with high average correlation in their dynamical evolution, quantified by changes in information entropy (Falasca et al., 2020).
Once domains are found and given a temporal window of length < , the average domain signal can be easily computed together with a weighted and direct network among domains, quantifying their connectivity patterns. By shifting the temporal window we can then evaluate changes in domains signals and in their linkages.
Given a domain , its signal ( ) is defined as the weighted cumulative anomaly of all time-series within that domain: where ( ) is a time series of length associated to grid cell with latitude and | | is the number of grid cells in .
The network is inferred by computing for each pair of domains and , their pairwise lagged-correlations , ( ) in a lag range ∈ [− , ] and by assess-ing their statistical significance. We account for the existence of autocorrelations (Guez et al., 2014) by adopting the Bartlett's formula (Box et al., 2011), and we solve the multiple testing problem employing the False Discovery Rate (FDR) proposed by Benjamini and Hochberg (1995). Two domains and are connected if there exists at least one significant correlation between the two at any lag in the range ∈ [− , ], denoted as , ( ). In this work we define the link directionality by the lag * which maximizes their lagcorrelation , ( ). If * = 0 the average variability of the two domains is considered as synchronous.
Average networks. Representing the linkages between domains is not trivial as links evolve over time. Nonetheless, it is often the case that two domains and are connected with a predominant directionality for a given period (i.e., 90% of the times → in the period from year 1 to 2 ). We leverage this characteristic to define the average network. Given a range of years ∈ [ 1 , 2 ] for which we inferred more than one network, for each pair of domains and , we focus on the most frequent link directionality. To each link is then associated a percentage of appearance, an averaged lag * and an average maximum significant correlation , ( * ) . Falasca et al. (2020) adopted ∈ [−12, 12] and a False Discovery Rate = 0.04. Therefore, in each time window, the number of false positive links is bounded to be at most 4%.

Dynamical system analysis
Climate is a high-dimensional system depending on multiple variables at both spatial and temporal scale. Currently, there is no framework allowing to study the system by accounting for its high-dimensionality in a comprehensive way. Here, exploiting tools from dynamical systems theory, we propose to study climate directly in its state space (Cvitanović et al., 2016). Consider a spatiotemporal climate system composed by fields Y , = 1, 2, ..., , embedded in a computational grid with points and with temporal length . For a given field Y , we first weight each time series by the cosine of its latitude and we center it to zero mean, then we standardize each field to unit variance. At each time step , the system is uniquely described by a state space vector In this formulation, the evolution of the system's dynamics is represented as a single-point trajectory in state space. The climate system is a dissipative dynamical system (Ghil and Lucarini, 2020), and therefore its dynamics should be confined to a lower (than the full state space) dimensional manifold. Quantifying the geometry of the underlying manifold allows for a comprehensive characterization of the system's dynamics.
Below we present two manifold learning algorithms: the well-known Principal Component Analysis (Hotelling, 1933) and newer Isomap (Tenenbaum et al., 2000), respectively for linear and nonlinear manifold learning.

Principal Component Analysis
Principal Component Analysis (PCA) (Hotelling, 1933) (commonly known as Empirical Orthogonal Functions (EOF) in climate science (Storch and Zwiers, 1999)) is a linear methodology for modal decomposition of spatiotemporal datasets. Consider a spatiotemporal dataset X ∈ R T,N with time series each of length . PCA aims at discovering the underlying manifold by fitting hyperplanes in directions explaining most of the variance. To do so we first compute the Gram matrix as G = 1 T−1 XX T ∈ R T,T (Bueso et al., 2020). The eigenvectors U ∈ R T,T of the Gram matrix G are the Principal Components (PCs) of the dataset . To each U we can assign an explained variance based on the ratio of its correspondent eigenvalue and the sum of all eigenvalues, i.e. / . To identify the correspondent spatial patterns, known as EOFs, it is sufficient to linearly regress each (standardized) PC on the dataset as 1 T−1 U T X. The low dimensional projection found by PCA preserves the variances as measured in the high-dimensional input X (Tenenbaum et al., 2000); similarly, instead of the covariance matrix it is possible to eigendecompose a distance matrix (containing the distances between each point in X) using the multidimensional scaling (MDS) algorithm (Borg and Groenen, 1997). Crucially, if the distance matrix is computed using Euclidean distances, the MDS and PCA solutions are equivalent (Tenenbaum et al., 2000).

Isometric feature mapping (Isomap)
If the system's dynamics lives on a curved manifold, Euclidean distances between points cannot capture the intrinsic distances along the manifold. Isomap, first introduced in Climate Science by Hannachi and Turner (2013), addresses this issue by adding an additional step to the MDS algorithm (Borg and Groenen, 1997). Given a dataset X ∈ R T,N , with time series each of length , and each centered to zero mean, Isomap first identifies the -nearest neighbors of each point in the manifold. By assuming that the manifold is locally flat in a radius of points, it computes the geodesic distances , between each couple of points and . The geodesic distance matrix D is then used as input to the MDS algorithm (Borg and Groenen, 1997). Given D , the double centered distance matrix is computed as A = − 1 2 JD J, with J = I − 1 ee T ; where I is the identity matrix of order and e is a vector of length containing all ones. The dimensionality of the dataset is then reduced by finding the eigenvectors of the double centered matrix A (i.e., solving A = U U T ). Each component (i.e., axis) V is Alternatively, it is possible to eigen-decompose the covariance matrix C = 1 N−1 X T X ∈ R N,N . In this case the eigenvectors V ∈ R N,N of C are spatial patterns. The projection of V onto X describes their temporal variability. The decomposition of the Gram matrix returns the same solution of the covariance matrix up to eigenvectors (Bishop, 2006;Bueso et al., 2020). then obtained by weighting the -th eigenvector U by the square root of its correspondent eigenvalue Tenenbaum et al., 2000). Similarly to PCA we can retrieve the associated spatial patterns by linearly regressing the (standardized) Isomap components on the original dataset as 1 T−1 V T X.
Differently from PCA, nonlinear dimensionality reduction methods do not have a direct way to compute the explained variance of the dataset. A tentative method, proposed by Tenenbaum et al. (2000) is to compute the residual variance as 1 − ( , ) 2 , where is the algorithm's estimate of manifold distances (i.e, Euclidean distance matrix for PCA and geodesic matrix for Isomap), is the matrix of Euclidean distances in the low dimensional embedding computed by each algorithm and is the Pearson correlation coefficient among all entries of and .
4 Mid-to late Holocene changes in variability in the Indo-Pacific system

Complex Network Analysis
The evolution of major modes of interannual-to-decadal variability and their connectivity patterns in the Indo-Pacific system is quantified through a complex network framework (see section 3). In this framework, modes (hereafter domains) correspond to spatially contiguous sets of sea surface temperature (SST) grid cells with highly correlated temporal evolution. Given a time window of years, the variability of each domain is quantified by its signal computed as the average monthly anomaly inside it.
Here, we focus on three domains in the Indo-Pacific region identified in Falasca et al. (2020) and quantify the evolution of the standard deviation of each signal along with their underlying correlation network in time windows of length = 100 years every = 3 years.
The three domains, shown in Figure 2a, correspond to an ENSO-like domain and two domains in the IO, covering respectively the western (WIO) and eastern (EIO) tropical part of the basin . Figure 2b shows that the standard deviation of the ENSO and WIO signals increases over time, suggesting that the thermodynamic feedback observed today (Lau and Nath, 2003) strengthened throughout the Holocene. The SST variability in the EIO domain, on the other hand, controls the evolution of the equatorial dipole index (EDI), computed as the difference of WIO and EIO signals. The variability of EIO domain weakens steadily from the mid-Holocene to about 3 kyr BP, and oscillates thereafter around a nearly constant value. These slow shifts in the variability translate into slow changes in domains connectivity, quantified, by the maximum significant To compare the different domains signals, we compute each one of them as the average time series inside it. It should be noted that, while all points inside each domain share high correlation values in temporal variability, the variance can show distinct regional differences. This is the case of ENSO, that, in the IPSL model, shows homogeneous variability in the tropical Pacific and larger variance only in the east and central Pacific. Therefore Figure 2b should not be used to compare magnitudes of the standard deviations, but rather their temporal changes. Given two domains, we plot only the link with the largest frequency which is also indicated (e.g., 99% for E→ WIO in panel (d) means that 99% of the links between those two domains are directed from E to WIO). The link color quantifies the average correlation between two domains and * the average lag.
correlation * in the [−12, 12] months lag range (Figure 2c). The increasing variability of ENSO throughout the Holocene translates into a stronger connectivity with the WIO domain, as shown by the small positive trend in the correlation values. More importantly, significant negative correlations between the EIO and WIO domains characterize the 6000 -4000 BP period, when an equatorial dipole controlled the IO variability. Over the same period, correlations between the ENSO and EIO domains were small and sporadic, indicating that the variability of the equatorial dipole was largely decoupled from that of the tropical Pacific. This decoupling is supported by work in coral reconstructions, showing that in the mid-Holocene the eastern IO was characterized by stronger than today upwelling events (Abram et al., 2007), together with reduced ENSO variability (Moy et al., 2002). Over time, the negative correlation between the WIO and EIO domains decreased and then changed to positive. The dominant mode of IO SST variability shifted from an equatorial dipole to a basin mode, and in concert the EIO variability became correlated with that of ENSO through its link with the WIO domain (red and black curves in Figure 2c after 3 kyr BP).
Figures 2d-f summarize the mid-to late Holocene connectivity changes between the ENSO, WIO and EIO domains in an average network. For a given time window and for each pair of domains we consider only the link directionality appearing with the largest frequency. The reported frequency therefore quantifies the stability of each link. The average lag * indicates the link directionality, while the mean correlation * for each pair of domains quantifies the strength of the connections. The strongest and most stable connectivity is between ENSO and WIO, with a robust lead of about 4 months from the Pacific to the IO domain. In the first 2000 years, 82% of the links between EIO and WIO are anticorrelated with the EIO leading. In the 4000 -2000 BP window, 82% of the links between those two domains are synchronous and correlations shift from negative to positive. Finally, in the most recent two millennia, 100% of the correlations between EIO and WIO are positive and maximized at zero lag, marking the transition to the IOB as dominant mode of variability.
Such dramatic, non-abrupt changes in climate modes are an exception rather than the rule. For example, a preliminary analysis focused on ENSO and North Pacific domains revealed quasi-stationary linkages (not shown here) and, therefore more subtle changes. In addition, Falasca et al. (2020) quantified changes in the global SST climate network through the evolution of the network density metric and showed that, on average, abrupt shits dominate changes in connectivities.

Main modes of interannual-to-decadal variability
We further compute the Empirical Orthogonal Functions (EOFs) (see PCA in section 33.2 and Storch and Zwiers (1999)) of the detrended monthly SST anomalies in the IO and tropical Pacific. The first EOF mode for the tropical Pacific and the first two for the tropical IO are shown in Figure 3 and 4 for the first and last 100 years of the simulation (a comparison with observations and the seasonality of the IO modes can be found in Appendix A, Figures 9-11).   (Figure 2b), which was highlighted by Crétat et al. (2020) as important to explain increased Indian monsoon variability. This pattern further identifies the central Pacific as the epicenter of these variability changes.
In the IO basin changes are more dramatic (Figure 4). In the most recent period (0.1-0 kyr BP) the first and second modes are the IOB and IOD, respectively, in agreement with observations (see Appendix A, Figure 10). In contrast, the mid-Holocene was subjected to a strong equatorial dipole mode (EDM) (see Figure 4 at 6-5.9 kyr BP) mainly independent from ENSO (see Figure 2c) and with a spatial distribution that matches that found by Thirumalai et al. (2019) in future projections and by DiNezio et al. (2020) during the LGM. The second EOF in the mid-Holocene portrays basin-wide same-sign SST anomalies ( Figure 4) that peak in boreal summer (Appendix A, Figure 11). As time passes, most of the variance shifts to the IOB, with the EDM as second mode. In current times, the EDM is replaced by the IOD with strong variability south of Sumatra rather than at the Equator (see Figure 4).
The EOF analysis highlights that the tropical Pacific variability strengthened over time, and that in the mid-Holocene the IO variability was modulated primarily by the EDM rather than by the IOB we observe today.
In summary, the complex network and EOF analysis highlight a dramatic change in Pacific and IO SST modes of variability and their connectivity patterns. These changes are non-abrupt and gradual and linked to the slowly changing Indo-Pacific system's mean state. This causal link and its multi-variable dependence is unveiled in the next two sections.

The Indo-Pacific as a high-dimensional dynamical system
To characterize the transition in SST variability, and motivated by recent progress in the study of high-dimensional systems, such as turbulence (Gibson et al., 2008), neuroscience (Gallego et al., 2020) and animal behavior (Ahamed et al., 2020), we construct a state space representation of the Indo-Pacific (i.e., [30S -30N , 40E -70W]) dynamics. For the last 6000 years we consider the evolution of a minimal set of variables that are key to tropical dynamics: ST, VP at 200hPa, U and V at 850 hPa and Precip, at yearly temporal resolution (see section 2).
At each time step t, the system is uniquely described by a state space vector X = [ST( , ), VP( , ), U( , ), V( , ), Precip( , )] ( ). In this analysis we retain trends induced by the external forcing to quantify drifts in dynamics. Given the model resolution (see section 2) and the variables considered, the reconstructed state space is 23,730-dimensional, and we visualize its high-dimensional dynamics by low dimensional projections as described in section 33.2. Figures 5(a,b) show the low dimensional projections using PCA and Isomap , respectively. For Isomap we set = 50, therefore assuming that the manifold is flat in a local radius of 50 years. Each point in Figure 5 is color-coded with the corresponding year from 6 kyr BP (mid-Holocene) to 0 kyr BP (year 1950). The explained variance of each principal component (PC) is ∼ 18%, ∼ 4%, ∼ 3% respectively. The residual variances (see section 3.2) of the first 3 components are ∼ 15% and ∼ 24% for PCA and Isomap, respectively. A dynamical shift is apparent, with different basins of attraction at the start and end of the period. This shift is gradual and linear at first order and, therefore linked to slow changes in the Earth's orbital configuration. The linearity (at least in the main components) is quantified by the strong similarities between PCA and Isomap projections, with correlation coefficients higher than 0.9 for the first three components.
To further investigate the spatial characteristic of these changes we compute the spatial patterns by linearly regressing each (standardized) component on all fields (see section 3.2). This is shown in Figure 6 for PCA and in Appendix B for Isomap. Each standardized component is reported in Figure 6a. The first PC captures the intensification of an ENSO-like pattern, in agreement with Paleo-records (Moy et al., 2002;Emile-Geay et al., 2016;Thompson et al., 2017;Grothe et al., 2020). It manifests as an increase in surface temperature in the tropics (see left panel in Figure 6b), driving stronger convection in the central Pacific as quantified by negative anomalies of VP at 200hPa (see left panel in Figure 6c) together with more precipitation in the Pacific The dimensionality reduction steps are performed via the Python Scikit-learn open source package (Pedregosa et al., 2011).
This analysis considers both trends and variability. For this reason here we refer to this pattern as "ENSO-like" and not to pure ENSO variability. sector (see left panel in Figure 6d) and a wind response to the warmer equatorial Pacific (Figure 6(e,f) ).
The second PC shows a steep trend in correspondence of warming over continental India, in the equatorial eastern Indian Ocean (IO) and in the south and central Pacific. The VP, Precip, U and V fields all point to changes in the equatorial Indian Ocean as important contributors to this signal. The warming of the eastern equatorial IO is followed by an increase in convection (middle panel in Figure 6c) and more precipitation (middle panel in Figure 6d). At the same time we see a strengthening of the equatorial westerly (middle panel in Figure 6e).
The third PC reveals a non-trivial trend coming from an Indian Ocean Dipole (IOD)-like variability, also apparent in the surface temperature and precipitation projections (right panel in Figure 6(b,d)).
This analysis identifies a clear transition of the system's dynamics from mid-to late Holocene emerging from at least 3 different regions: the central Pacific Ocean, the equatorial Indian Ocean and an IOD-like variability. Additionally, the time series in Figure 6 qualitatively show that, while the long-term dynamics (i.e., on the order of thousands of years) are constrained by orbital changes, the system's trajectory over periods as long as 500-year is dominated by internal vari- ability (see also Appendix B, Figure 13). This last result supports recent work showing that the internal variability of the climate system from the mid-Holocene to present times can mask any long-term trend for periods up to about 500 years (Braconnot et al., 2019a).

Mid-to late Holocene changes in the Indo-Pacific mean state
The dynamical shift is further investigated by considering changes in the IO mean climate over windows of 500 years. We examine the departure of 500-yr annual mean snapshots from the annual mean conditions averaged over 6-5.5 kyr BP. The annual changes are the residual of seasonal changes, which are briefly discussed for JJAS (June to September) and DJFM (December to March) in Appendix C. At 6-5.5 kyr BP, low-level wind convergence, deep convection and 200-hPa divergence prevail over the western Pacific, while opposite conditions characterize the region from the eastern Pacific into the west IO (Figures 7, 8: top panel). Most importantly, the annual mean view of the Walker circulation hides north-south displacements along with the seasonal cycle driven by insolation. Over the western-central Pacific, the largest 200-hPa divergence is in the northern hemisphere in JJAS (Appendix C, Figure 14: VP, top panel) and in the southern hemisphere in DJFM (Appendix C, Figure 15: VP, top panel). Importantly, mid-Holocene conditions promote an equatorial west-to-east surface temperature gradient of roughly 1 • K in the equatorial IO (Figure 7, top panel). The strengthening of the Walker circulation from mid-to late Holocene is reflected in the 200-hPa divergence over the central Pacific and convergence over the WIO (Figure 7) and is linked to atmospheric circulation changes in JJAS (Appendix C, Figure 14). At least two factors contribute to these changes. Firstly, the warming  (Figure 7). This translates into enhanced upper-level convergence and subsidence in the WIO, strengthens low-level westerly winds over the IO (Figure 8) and, over time, reduces the surface temperature gradient in the equatorial IO. The monsoon, on the other hand, is sharply damped by the orbitally-driven weakening of the inter-hemispheric energy gradient during boreal summer from mid-to late Holocene (Braconnot et al., 2019a), consistent with drying and warming trends over India (Figure 7 and Appendix C, Figure 14). The monsoon weakening results in weaker southerly winds over the Arabian Sea and in a southward shift of the zonal mean circulation, with westerlies weakening over land masses and strengthening along the equatorial IO ( Figure 8 and Appendix C, Figure 14). This in turn promotes low-level divergence in the WIO and low-level convergence over the central Pacific, hence contributing to reinforce the Walker circulation.
While changes in the atmospheric circulation over the Pacific and the IO are large, the oceanic connectivity between the two basins (Bracco et al., 2005;Sprintall and Révelard, 2014) does not vary significantly. The simulated heat transport (HT) through the Indonesian Throughflow (ITF) has a negative trend, with more heat entering the IO from the Pacific, but of very small magnitude (a change of ∼0.02 PW or less than 2% over 6000 years) (see Appendix D, Figure 16).
All the above highlights the amplification of Pacific-IO connectivity through the Walker circulation from mid-to late Holocene and translates in a reduced surface temperature gradient in the IO basin along the equator. As a result, the climate evolved towards current IO basin conditions favoring an IOB-like mode of variability.

Conclusions
We investigated the transient evolution of the Indo-Pacific climate system and its modes of variability over the last 6000 years. Gradual changes in the external forcing and especially in the Earth's orbital configuration induced slow, non-uniform changes in the climate mean state and, in turn, in its interannual modes of variability from the mid-Holocene to 1950. Through complex network analysis we showed that the mean state slow changes set the stage for profound shifts in variability by strengthening ENSO and the IO basin mode. We found that in the mid-Holocene, the strong SST gradient in the equatorial IO favored the EDM, which is independent of ENSO and similar, in its dynamics, to ENSO and the Atlantic Niño (Thirumalai et al., 2019;DiNezio et al., 2020). The EDM was present during the LGM (DiNezio et al., 2020) and according to our analysis, supported by coral paleo-reconstructions, remained the dominant IO mode of variability well into the mid-Holocene. The slow changes in the mean state damped the east-west surface temperature gradient that fuels the EDM inducing an increase in strength and relevance of the IOB around 4000 years ago, and the slow morphing of the EDM into the IOD we observe today. The state space representation of the Indo-Pacific climate system indicates that its trajectory explored different basins of attraction from the mid-Holocene to present times. The drift from one basin to the other was linear at the first order, and gradual, with the internal climate variability masking the mean state drift for hundreds of years (see Figure 13) (Braconnot et al., 2019a). The linearity of the transition was quantified by comparing linear and nonlinear manifold learning algorithms. The spatial signatures of the mean-state changes are a persistent warming over the tropical Pacific and a weakening of the Indo-Asian monsoon that drove a strengthening of the Walker circulation. This resulted in a more active convective cell over the Pacific Ocean and an intensification of the surface westerly winds in the IO basin, which dampened the SST gradient in the equatorial IO. In conclusions, slow changes in orbital (external) climate forcing caused gradual shifts in both mean state and connectivity in the tropical Indo-Pacific from mid-to late-Holocene. While gradual, they accounted for large and dramatic changes in surface temperature, wind and precipitation patterns, among other fields, and influenced some of the most populated areas of the planet. A fast departure from the basin of attraction in the first half of the 20th century is already evident if only surface temperature (see Figure 13) is considered and more rapid shifts in variability may follow.
Finally, the methodologies adopted in this study, borrowed from computer science and dynamical systems theory, allow to better visualize and understand transient climate change. They are powerful tools to delve into the connections between mean state and modes of variability, and to visualize their evolution over time and across multiple fields simultaneously. Most importantly, they allow for identifying both linear and nonlinear responses, and for evaluating their representation in Earth System models.
Appendix A: Sea Surface Temperature Modes: model evaluation, mid-to late Holocene evolution and Indian Ocean seasonality

A.1: Model evaluation
We qualitatively compare the Sea Surface Temperature (SST) variability in the Pacific and Indian Ocean (IO) basins from the IPSL model (latest 100 year period: 1851-1950) and the 1/4 • ERA5 reanalysis (Hersbach et al., 2020) in the 1980-2019 period, for which we have satellite coverage. The comparison is done by comparing the first Empirical Orthogonal Function (EOF) in the Pacific Ocean and the first two in the IO computed using monthly SST anomalies after removing seasonality and linear trend.
Under current climate conditions, the ERA5 interannual-to-decadal SST variability is largely dominated by ENSO ( 63% of the explained variance) in the tropical Pacific ( Figure 9). The IPSL modelled ENSO extends too much into the western tropical Pacific with almost uniform values in the tropical Pacific. The shape of the modelled ENSO is therefore not perfect, but, similarly to ERA5, the first EOF explains 60% of tropical Pacific variability. In the tropical IO, the first two EOFs for ERA5 are, respectively, the IO Basin (IOB) mode (Klein et al., 1999) and the IO Dipole (IOD) mode (Webster et al., 1999;Saji et al., 1999) (Figure 10). Biases in the model include negative anomalies south of Sumatra in the first EOF and in the north-eastern IO in the second EOF. Nonetheless, the IPSL model accurately captures the observed first two modes in the IO. Here, we investigate changes in seasonality of these three modes from mid-to late Holocene. The seasonal analysis of SST modes in the tropical IO is shown in Figure D3 for four 100 years long periods. Three of them are in the first 1000 years of the simulation (i.e., 6-5.9, 5.5-5.4 and 5.0-4.9 kyr BP) when the EDM experiences rapid changes in its variance. The last period covers the 0.1-0 kyr BP period and is used to point to the model ability to capture the observed seasonality of the IOB and IOD under current climate conditions. Interestingly, the EDM persistently peaks in boreal fall from 6 to 5 kyr BP, despite its decreased variance. The seasonality of the EDM is close to that of the IOD under current climate conditions. This suggests that EDM and IOD are driven by similar ocean-atmosphere feedbacks (e.g., Bjerkness) but that the modulation in the mean state of the zonal SST gradient leads to different flavors of the zonal dipole.
The seasonality of the IOB drastically changes from mid-to late Holocene. At 6-5.9 kyr BP, the IOB proves to be weak, with no clear basin-wide structure regardless of the season. As ENSO strengthens, the IOB structure rapidly emerges in boreal spring to become prominent almost throughout the annual cycle by late Holocene.  In Figure 12 we present the spatial patterns obtained using Isomap for dimensionality reduction. For the Isomap projection we set = 50, therefore assuming that the manifold is flat in a local radius of 50 years. In Figure 13 we present the two Principal Components of only the Surface Temperature (ST) dataset ( Figure 13a) and the state space representation of our system (Figure 13b). In this analysis we focus on 500 years-long windows and on the last 50 years. We make two observations: (i) while the long-term dynamics is constrained by orbital changes (see section 5) the system's trajectory in periods as long as 500-year is dominated by internal variability and no clear drift is identified; (ii) the state space dimensionality reduction differs from the one done on solely ST. Correlations between the first two PCs in surface temperature and state space embedding are 0.91 and 0.84 respectively. Importantly, the first PC in ST shows an exponential drift in the last 50 years, in correspondence to the large CO 2 increase in the trace gas forcing field (see Figure 1). This is not as obvious in the state space representation and can be seen only in the second component (see Figure A2 and Figure 6a).

Appendix C: Mid-to late Holocene changes in the Indo-Pacific mean state: seasonal dependence
Seasonal changes in the Indo-Pacific as departures of 500-yr mean snapshots in boreal summer and winter from conditions averaged over 6-5.5 kyr BP. In what follows, averages from June to September are indicated as JJAS and from December to March as DJFM. Importantly, the main mid-to late Holocene changes in yearly circulation highlighted in this paper (section 6) are dominated by the signals in JJAS. Figure 16 shows changes in the climatology of heat transport (HT) through in the Indonesian Throughflow (ITF) over the whole water column. The analysis is performed by computing the running mean for time windows of 500 years every 1 year. We find a statistically significant strengthening of heat transported from the Pacific to the IO from mid-to late Holocene. The changes, however, are small in magnitude ( 0.02 PW or less than 2% over the entire 6000 years period) indicating a rather small oceanic contribution to changes in the Pacific-IO dynamics, at least in the model analyzed and at the annual timescale.