I think that the answer is yes, and it follows from the following claim~:

If $f:[0,\infty)\to \mathbf R$ is smooth function with first derivatives vanishing fast enough at infinity (say for simplicity that $f$ is in the Schwarz class), then the matrix $$\widetilde A_{k,j} = \varepsilon f(\varepsilon(k+j))$$ belongs to the trace class with norm independant from $\varepsilon$.

How is this related to your question? Well, take $f(x) = 4x e^{-x^2}$ and $\varepsilon = \sqrt t$. Then using that $e^{-t n^2} - e^{-t (n+2)^2} = f(\sqrt n) + 0((t+t^2n^2) e^{-tn^2})$, you see that your matrix $A$ and the matrix $\widetilde A$ have a difference uniformly bounded in the trace class (for example because the $\ell^1$-norm of the coefficients of the difference is bounded). And by the claim $\widetilde A$ belongs to the trace class uniformly in $\varepsilon$.

Why is the claim true? It follows from Peller's theorem quoted in the question, and from the elementary inequality $\| \varphi \|_{L^1(\mathbf T)} \leq c \sqrt{ \|\varphi\|_{L^2} \|(1-z)\varphi\|_{L^2}}$ for every function $\varphi$ on the unit circle $\mathbf T$ with Lebesgue measure.

Indeed, Peller's theorem characterizes, up to uniform constants, the trace class norm of a Hankel matrix $A_{k,j} =\widehat \varphi(k+j)$ by the norm of $\varphi = \sum_{k \geq 0} \widehat \varphi(k)z^k$ in the Besov space $B_{1,1}^1$. Namely if $w \in C_c(0,\infty)$ is a smooth function satisfying $\sum_{n \in \mathbf Z} w(2^n x)=1$ (say that $w$ is supported in $[1/2,2]$ for the sequel) and $W_n(z) = \sum_k w(k/2^n) z^k$ for $n>0$ and $W_0(z) =1+z$, $$ \|\widetilde A\|_{S^1} \simeq \sum_{n \geq 0} 2^n \| W_n \ast \varphi\|_{L^1(\mathbf T)}$$

Now if $\widehat \varphi(k) = \varepsilon f(\varepsilon k)$ then $\|W_n \ast \varphi\|_2$ is dominated, up to a constant and small error that I forget here, by $\sqrt \varepsilon (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2)^{1/2}$. Similarly $\|(1-z)W_n \ast \varphi\|_2$ is dominated by $2^{-n} \sqrt \varepsilon (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2)^{1/2} + \varepsilon^{3/2} (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f'|^2)^{1/2}$. (the first term comes from the derivative of $w(\cdot/2^n)$ and the second from the derivative of $f(\varepsilon \cdot)$). By the "elementary inequality", we get the bound $$\|W_n \ast \varphi\|_{L^1} \lesssim C\sqrt \varepsilon (\sqrt \varepsilon + \sqrt{2^{-n}}) \left(\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2+|f'|^2\right)^{1/2}+\textrm{ small error}.$$

Then at least if one forgets the error, one gets for $\delta>0$ $$ \sum_n 2^n \|W_n \ast \varphi\|_{L^1} \lesssim \left(\int (1+x^{2+\delta})(|f(x)|^2+|f'(x)|^2)dx\right)^{1/2}.$$

I think that the answer is yes.

Edit: as noticed in the comments, the answer would be obviously yes if the matrix $A$ was positive, since in this case its trace norm (denoted $\|A\|_1$) would be equal to its trace which is less than $1$. But $A$ is not positive: when $t$ goes to $0$, $A_{i,j}/4t$ goes to to $1+i+j$, which is not a positive matrix since every $2\times 2$ submatrix has determinant $(1+2i)(1+2j)-(1+i+j)^2<0$.

One therefore has to be more careful, and the fact that $\sup_t \|A\|_1 <\infty$ follows from the following claim: (end of Edit)

If $f:[0,\infty)\to \mathbf R$ is smooth function with first derivatives vanishing fast enough at infinity (say for simplicity that $f$ is in the Schwarz class), then the matrix $$\widetilde A_{k,j} = \varepsilon f(\varepsilon(k+j))$$ belongs to the trace class with norm independant from $\varepsilon$.

How is this related to your question? Well, take $f(x) = 4x e^{-x^2}$ and $\varepsilon = \sqrt t$. Then using that $e^{-t n^2} - e^{-t (n+2)^2} = f(\sqrt n) + 0((t+t^2n^2) e^{-tn^2})$, you see that your matrix $A$ and the matrix $\widetilde A$ have a difference uniformly bounded in the trace class (for example because the $\ell^1$-norm of the coefficients of the difference is bounded). And by the claim $\widetilde A$ belongs to the trace class uniformly in $\varepsilon$.

Why is the claim true? It follows from Peller's theorem quoted in the question, and from the elementary inequality $\| \varphi \|_{L^1(\mathbf T)} \leq c \sqrt{ \|\varphi\|_{L^2} \|(1-z)\varphi\|_{L^2}}$ for every function $\varphi$ on the unit circle $\mathbf T$ with Lebesgue measure.

Indeed, Peller's theorem characterizes, up to uniform constants, the trace class norm of a Hankel matrix $A_{k,j} =\widehat \varphi(k+j)$ by the norm of $\varphi = \sum_{k \geq 0} \widehat \varphi(k)z^k$ in the Besov space $B_{1,1}^1$. Namely if $w \in C_c(0,\infty)$ is a smooth function satisfying $\sum_{n \in \mathbf Z} w(2^n x)=1$ (say that $w$ is supported in $[1/2,2]$ for the sequel) and $W_n(z) = \sum_k w(k/2^n) z^k$ for $n>0$ and $W_0(z) =1+z$, $$ \|\widetilde A\|_{S^1} \simeq \sum_{n \geq 0} 2^n \| W_n \ast \varphi\|_{L^1(\mathbf T)}$$

Now if $\widehat \varphi(k) = \varepsilon f(\varepsilon k)$ then $\|W_n \ast \varphi\|_2$ is dominated, up to a constant and small error that I forget here, by $\sqrt \varepsilon (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2)^{1/2}$. Similarly $\|(1-z)W_n \ast \varphi\|_2$ is dominated by $2^{-n} \sqrt \varepsilon (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2)^{1/2} + \varepsilon^{3/2} (\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f'|^2)^{1/2}$. (the first term comes from the derivative of $w(\cdot/2^n)$ and the second from the derivative of $f(\varepsilon \cdot)$). By the "elementary inequality", we get the bound $$\|W_n \ast \varphi\|_{L^1} \lesssim \sqrt \varepsilon (\sqrt \varepsilon + \sqrt{2^{-n}}) \left(\int_{2^{n-1}\varepsilon}^{2^n \varepsilon} |f|^2+|f'|^2\right)^{1/2}+\textrm{ small error}.$$

Then at least if one forgets the error, one gets for $\delta>0$ $$ \sum_n 2^n \|W_n \ast \varphi\|_{L^1} \lesssim \left(\int (1+x^{2+\delta})(|f(x)|^2+|f'(x)|^2)dx\right)^{1/2}.$$