In astrophysics, controlled experiments are typically impossible, and it is then necessary to make the most of observational data.Other disciplines that are in a similar predicament --- from epidemiology to economics --- increasingly leverage causal inference methods.This is however not yet the case in astrophysics. In this contribution, we apply causal discovery for the first time to an important open problem in astrophysics, namely the possible coevolution of supermassive black holes (SMBHs) and their host galaxies.We make use of a comprehensive catalog of observed galaxy properties, on which we apply the Peter-Clark (PC) algorithm to obtain a single completed partially directed acyclic graph (CPDAG), representing a Markov equivalence class over directed acyclic graphs (DAGs). We test the robustness of our analysis by randomly subsampling our dataset and showing that we recover similar results.We suggest a physical explanation for the causal structure that we learned in terms of the hierarchical assembly pathway of SMBHs.