Compute the fourth moment of multivariate Gaussian in Linear Regression

Let $X \in \mathbb{R}^d \sim N(0, \Sigma)$ be a multivariate Gaussian random vector. For any $\beta\in \mathbb{R}^d$, consider the regression model $$Y = X^\top \beta + \epsilon,$$ where $\epsilon \sim N(0,\sigma^2)$ is the random noise. How to compute expectation $\mathbb{E} [ Y^2 XX^\top ]$?

The vector $(X,Y)$ is Gaussian, because any sum is normal: $$\sum a_iX_i + bY= \sum c_i X_i + b\epsilon\sim\mathcal N(\mu,\sigma^2)$$ for some mean and variance. (assuming $\epsilon$ is independent of all other variables.)
If $(X,Y)$ is a zero-mean multivariate Gaussian, then we can use Isserlis' theorem: For some $i,j$ let us calculate the scalar $$E[Y^2X_i X_j]=E[Y\cdot Y\cdot X_i \cdot X_j]$$ $$E[Y\cdot Y\cdot X_i \cdot X_j]=E[YY]E[X_iX_j]+E[YX_i]E[YX_j]+E[YX_i]E[YX_j]$$
Everything here should be known: $E[X_iX_j]=(\Sigma)_{i,j}$, $$E[YX_i]=E[X_iX^T\beta ] + E[\epsilon X_i]=\sum E[X_iX_k]\beta_k+0,$$ and $$E[Y^2]=E[(X^T\beta+\epsilon)^2]=E[\epsilon^2]+0+\beta^T E[X X^T]\beta.$$