A Bayesian network (BN) is a probabilistic graphical model that consists of a directed acyclic graph (DAG), where each node is a random variable and attached to each node is a conditional probability distribution (CPD, and is used for knowledge representation and modelling causal relationships between variables. A BN can be learned from data using the well-known score-and-search approach, and within this approach a key consideration is how to simultaneously learn the global structure in the form of the underlying DAG and the local structure in the CPDs. Several useful forms of local structure have been identified in the literature but thus far the score-and-search approach has only been extended to handle local structure in form of context-specific independence and linear models. In this work, we show how to extend the score-and-search approach to Multivariate Adaptive Regression Splines (MARS), which are regression models represented as piecewise spline functions. MARS models capture non-linear relationships without imposing a non-linear representation over the global structure. We provide an effective algorithm to score a candidate MARS using the widely used BIC score and we provide pruning rules that allow the search to successfully scale to large sized networks without overfitting. Our empirical results provide evidence for the success of our approach to learning Bayesian networks that incorporate MARS relations.