Let us identify an involution $\sigma$ in $\mathfrak{S}_n$ with a set partition $\Pi_\sigma$ of $[n] := \{1,2,...,n\}$ into nonempty blocks of size at most two in the obvious way: we have $\{a,b\} \in \Pi_\sigma$ if and only if $\sigma(a)=b$. This is obviously a bijection between involutions and such set partitions; and the fixed points of $\sigma$ correspond exactly to the singletons of $\Pi_\sigma$.
Now to get a set partition $\Pi'$ of $[n+2]$ with block sizes at most two from such a set partition $\Pi_\sigma$ of $[n]$ corresponding to the involution $\sigma \in \mathfrak{S}_n$, we just have to decide what blocks the elements $n+1$ and $n+2$ belong to:
- they could each be singletons (one way);
- they could be in a size-two block together (one way);
- $n+1$ could form a block of size two with something that was a singleton in $\Pi_\sigma$, while $n+2$ is a singleton ($\mathrm{tr}(\sigma)$ many ways);
- $n+2$ could form a block of size two with something that was a singleton in $\Pi_\sigma$, while $n+1$ is a singleton ($\mathrm{tr}(\sigma)$ many ways);
- $n+1$ could form a block of size two with something that was a singleton in $\Pi_\sigma$, and $n+2$ could form a block of size two with a different singleton in $\Pi_\sigma$ ($\mathrm{tr}(\sigma)(\mathrm{tr}(\sigma)-1)$ many ways).
Adding these all together, we get $1+1+\mathrm{tr}(\sigma)+\mathrm{tr}(\sigma)+\mathrm{tr}(\sigma)(\mathrm{tr}(\sigma)-1)=\mathrm{tr}(\sigma)^2+\mathrm{tr}(\sigma)+2$ such set partitions for each involution $\sigma \in \mathfrak{S}_n$, proving your claimed identity.
EDIT:
Here is the formula for the more general case of $I_{n+m}$:
$$ I_{n+m} = \sum_{\sigma \in \mathrm{Inv}(\mathfrak{S}_n)} \sum_{j=0}^{m} \binom{m}{j} \cdot I_{m-j} \cdot \frac{\mathrm{tr}(\sigma)!}{(\mathrm{tr}(\sigma)-j)!}.$$
It can be proved in the same manner: for each involution $\sigma \in \mathfrak{S}_n$, we choose $j$ of the elements of $\{n+1,n+2,\ldots,n+m\}$ to join with singletons in $\Pi_{\sigma}$, and hence choose an ordered subset of $j$ singletons from $\Pi_{\sigma}$ to pair with these elements, and then we put an involution on the remaining $m-j$ elements of $\{n+1,n+2,\ldots,n+m\}$ in $I_{m-j}$ ways (with the convention $I_0 =1$).
EDIT 2:
Maybe one final repackaging of this result is worth recording.
Set $I_n(t) := \sum_{\sigma \in \mathrm{Inv}(\mathfrak{S}_n)} t^{\mathrm{tr}(\sigma)}$ for all $n \geq 0$ (with $I_0(t)=1$). Then for all $m \geq 0$,
$$ I_{n+m}(t) = \left(\sum_{k=0}^{m}\binom{m}{k} \, I_{m-k}(t) \, \frac{d^k}{dt^k} \right) \, I_{n}(t).$$
In particular,
$$ I_{n+1}(t) = \left(t+\frac{d}{dt}\right) I_{n}(t);$$
hence we also have,
$$ I_{n+m}(t) = \left(t+\frac{d}{dt}\right)^m I_{n}(t).$$