Suppose we are given a dataset but not the capability of performing
some AB testing. We do some regression using X as predictor and Y as
response and get a model. Can we actually say something about the
causal relationship between X and Y?
No you can't, even when all variables are observed, see here for instance. If you are given only distributional information about the data (ie., you know the joint distribution of the observed variables), but no information about how the data was generated (a causal model), causal inference is impossible. In short, you need causal assumptions to get causal conclusions. You can get started on learning causal inference with the references here.
It is easy to understand why that is the case by constructing an example where different causal models entail the same observed joint probability distribution. Consider that you have observed the joint probability distribution $P(x, y)$ of two random variables. Here, imagine you have no sampling uncertainty---so you have perfect knowledge of $P(x, y)$, which entails perfect knowledge of the regression function and so on. To simplify things, consider that, in your data, $P(x,y)$ was found to be jointly normal with mean $0$, variance 1 and covariance $\sigma_{xy}$ (this is without loss of generality, you can always standardize the data). What can you say about the causal effect of $x$ on $y$ or vice-versa?
With only this information, nothing. The reason here is that there are several causal models that would create the same observed distribution, yet have different interventional (and counterfactual) distributions. Here I will show three of such models. Notice all of them gives you the same observed $\sigma_{xy}$, but their causal conclusions are different: in the first model $X$ causes $Y$, in the second model $Y$ causes $X$, and, in the third model, neither causes each other --- $X$ and $Y$ are both common causes of the unobserved variable $Z$.
Model 1
$$
X = u_{x}\\
Y= \sigma_{yx}x + u_{y}
$$
Where $U_{x} \sim \mathcal{N}(0, 1)$ and $U_{y} = \mathcal{N}(0, 1 - \sigma_{xy}^2)$.
Model 2
$$
Y = u_{y}\\
X = \sigma_{yx}y + u_{x}
$$
Where $U_{x} \sim \mathcal{N}(0, 1 - \sigma_{xy}^2)$ and $U_{y} = \mathcal{N}(0, 1)$.
Model 3
$$
Z = U_{z}\\
X = \alpha Z + U_{x}\\
Y = \beta Z + U_{y}
$$
Where $\alpha\beta= \sigma_{xy}$, $U_{z} = \mathcal{N}(0, 1)$, $U_{x} = \mathcal{N(0, 1- \alpha^2)}$ and $U_{y} = \mathcal{N(0, 1- \beta^2)}$.