$\textit{Lemma :}$ We have
\begin{align*}%$
\frac{1}{d!} \sum_{\lambda\in \text{Hooks}(d)} (-1)^{\ell(\lambda)-1} \, \dim \lambda \prod_{\Box \in \lambda}(x + c(\Box) )(y + c(\Box) ) = [c_d]\prod_{k = 1}^d (x + J_k) (y + J_k)
\end{align*}
where
\begin{align*}%$
J_1 = 0 \, ; J_2 = (1, 2) \, ; \qquad J_k = (1, k) + \cdots + (k - 1, k)
\end{align*}
are the Jucys-Murphy elements of the symmetric group $ \mathfrak{S}_d $, i.e. particular elements of the group algebra $ \mathbb{C}\mathfrak{S}_d $ that commute with each other and that generate the centre $ Z(\mathbb{C}\mathfrak{S}_d) $ in a certain way (symmetric functions in the JM elements generate it). Here we have noted $ [\sigma]f $ the coefficient of $\sigma$ in $ f \in \mathbb{C}\mathfrak{S}_d $, namely, if $ f = \sum_{ \sigma \in \mathfrak{S}_d } f_\sigma \sigma $, then $ [\sigma]f = f_\sigma $. Last, $ c_d $ is any $d$-cycle, for instance $ (1, 2, \dots, d) $.
$ $
$\textit{Proof.}$ We prove this last formula using the dual orthogonality of the characters of $ \mathfrak{S}_d $ : denote by $ \chi^\lambda $ such a character. Then, for all $ \sigma, \tau \in \mathfrak{S}_d $
\begin{align*}%$
\frac{1}{d!} \sum_{ \lambda \vdash d } \chi^\lambda(\sigma) \chi^\lambda(\tau) = \delta_{\sigma, \tau}
\end{align*}
We have used the fact that $ \chi^\lambda(\sigma) \in \mathbb{Z} $ to avoid writing a conjugation on one of the characters.
The next step is the formula (that one can prove with the Murnaghan-Nakayama formula)
\begin{align*}%$
\chi^\lambda(c_d) =: \chi^\lambda_{(d)} = (-1)^{\ell(\lambda) + 1} \mathbb{1}_{ \{ \lambda \in \mathrm{Hooks}(d) \} }
\end{align*}
We thus have
\begin{align*}%$
\delta_{\sigma, c_d} & = \frac{1}{d!} \sum_{ \lambda \vdash d } \chi^\lambda(\sigma) \chi^\lambda(c_d) \\
& = \frac{1}{d!} \sum_{ \lambda \in \mathrm{Hooks}(d) } (-1)^{\ell(\lambda) + 1} \dim(\lambda) \frac{\chi^\lambda(\sigma) }{\dim(\lambda)}
\end{align*}
Defining the normalised character by
\begin{align*}%$
\widehat{\chi}^\lambda := \frac{\chi^\lambda }{\dim(\lambda)} = \frac{\chi^\lambda }{ \chi^\lambda(id_d) }
\end{align*}
we then use the following (fundamental) formula that links content alphabets with Jucys-Murphy elements : for all symmetric function $ F $ in $d$ variables, for all $ \lambda \vdash d $
\begin{align*}%$
\widehat{\chi}^\lambda( F(J_1, \dots, J_d) ) = F(\mathcal{A}^\lambda)
\end{align*}
where $ \mathcal{A}^\lambda $ is the $\textit{content alphabet}$ defined by
\begin{align*}%$
\mathcal{A}^\lambda = \{ c(\square), \square \in \lambda \}
\end{align*}
My personal reference for this machinery is some articles of Biane (a priori probability theory...), but the people on this forum certainly know some better (more recent) exposition.
Using this and the fact that we have extended characters by linearity in the group algebra, i.e. $ \chi^\lambda(f) = \sum_{\sigma \in \mathfrak{S}_d} f_\sigma \chi^\lambda(\sigma) $ for $ f $ previously defined, one gets
\begin{align*}%$
[c_d]f & = \sum_{\sigma \in \mathfrak{S}_d} f_\sigma \delta_{\sigma, c_d} \\
& = \sum_{\sigma \in \mathfrak{S}_d} f_\sigma \frac{1}{d!} \sum_{ \lambda \vdash d } \chi^\lambda(\sigma) \chi^\lambda(c_d) \\
& = \frac{1}{d!} \sum_{ \lambda \vdash d } \chi^\lambda(f) \chi^\lambda(c_d) \\
& = \frac{1}{d!} \sum_{ \lambda \in \mathrm{Hooks}(d) } (-1)^{\ell(\lambda) + 1} \dim(\lambda) \widehat{\chi}^\lambda(f)
\end{align*}
The lemma is thus proven for the following particular element of the (centre of the) symmetric group algebra
\begin{align*}%$
F_{x, y}(J_1, \dots, J_d) := \prod_{k = 1}^d (x + J_k)(y + J_k)
\end{align*}
$ \square $
$ $
Now, to come back to your initial problem, we use the following Jucys-Murphy identity (or chinese restaurant theorem for the Ewens measure in probability theory)
\begin{align*}%$
\prod_{k = 1}^d (x + J_k) = \sum_{\sigma \in \mathfrak{S}_d} x^{C(\sigma)}\sigma
\end{align*}
with $ C(\sigma) $ the total number of cycles of $ \sigma $.
Taking the product in the group algebra, namely, the convolution product, one gets
\begin{align*}%$
\prod_{k = 1}^d (x + J_k)(y + J_k) & = \sum_{\sigma \in \mathfrak{S}_d} x^{C(\sigma)}\sigma \sum_{\tau \in \mathfrak{S}_d} y^{C(\tau)}\tau \\
& = \sum_{\sigma \in \mathfrak{S}_d} \left( \sum_{ \tau \in \mathfrak{S}_d } x^{ C(\tau) } y^{C(\sigma\tau^{-1})} \right) \sigma
\end{align*}
Taking the coefficient $ [c_d] $, we thus have
\begin{align*}%$
[c_d]\prod_{k = 1}^d (x + J_k)(y + J_k) = \sum_{ \tau \in \mathfrak{S}_d } x^{ C(\tau) } y^{C(c_d\tau^{-1})}
\end{align*}
You are concerned with the case $ x = y $, which gives
\begin{align*}%$
\frac{1}{d!} \sum_{\lambda\in \text{Hooks}(d)} (-1)^{\ell(\lambda)-1} \, \dim \lambda \prod_{\Box \in \lambda}(x + c(\Box) )^2 = \sum_{ \sigma \in \mathfrak{S}_d } x^{ C(\sigma) + C(c_d^{-1}\sigma) }
\end{align*}
If you want the coefficients of your polynomial, you get
\begin{align*}%$
\sum_{ \sigma \in \mathfrak{S}_d } x^{ C(\sigma) + C(c_d^{-1}\sigma) } = \sum_{k = 2}^{2d} a_{k, d} x^k, \qquad a_{k, d} := \#\{ \sigma \in \mathfrak{S}_d \, / \, C(\sigma) + C(c_d^{-1}\sigma) = k \}
\end{align*}
I personally doubt that this is simple to express. Maybe with some Hurwitz numbers or some free probability (the set defining $ a_{d + 1, d} $ defines the non-crossing permutations). You can try your recurrence relation with this new expression.
