Large-scale, persistent patterns of variability, such as ENSO and the Southern Annular Mode (SAM), have significant impacts on regional weather and climate around the globe.  Characterising the behaviour of and interactions between individual teleconnections may identify regimes of enhanced predictability, while understanding how these relationships are likely to change in response to anthropogenic warming is an important input in adaptation scenarios. Data-driven methods provide a complementary approach to dynamical and model-based teleconnection studies. However, traditional approaches tend to be built on underlying assumptions of stationarity and do not adequately capture the multiscale nature nor the secular trends associated with these modes. Recently developed non-stationary clustering methods attempt to remedy this deficiency by approximating the full dynamics in terms of a set of locally stationary regimes, with transitions between regimes governed by a latent set of persistent regime affiliations and relevant exogeneous variables, hence allowing the impacts of external forcings on the quasi-stationary states to be modelled. Here, we apply a variant of this approach to infer the causal relations among the major Northern and Southern Hemisphere teleconnection patterns from reanalysis data, highlighting the corresponding physical mechanisms. We then compare the structure of the resulting causal networks to those obtained by applying the method to output from the CMIP5 historical simulations in order to assess the extent to which the observed dependencies are reproduced in model runs. Finally, we discuss the evolution of the inferred networks under the RCP4.5 and RCP8.5 scenarios out to 2100.