Before I try to answer the broader question let me try to clear up a potential confusion. The truncations that appear at the level of modelled distributions don't really change the distribution they represent so long as we always truncate to a positive order. This is because if you have $0 < \gamma_1 < \gamma_2$, $f \in D^{\gamma_2}$ and you define $g \in D^{\gamma_1}$ by truncating $f$ at order $\gamma_1$, you will still have $\mathcal{R} f = \mathcal{R} g$. This follows from the uniqueness part of the reconstruction theorem since for $\lambda < 1$,
$$(\mathcal{R} f - \Pi_x g(x))(\phi_x^\lambda) = (\mathcal{R}f - \Pi_x f(x))(\phi_x^\lambda) - \sum_{|\tau| \in [\gamma_1, \gamma_2)} \Pi_x (f_\tau(x) \tau) (\phi_x^\lambda).$$
The first term on the right hand side is bounded by a multiple of $\lambda^{\gamma_2} \le \lambda^{\gamma_1}$ by the reconstruction theorem for $f$ and the second term is bounded by $\lambda^{\gamma_1}$ by the definition of the model. Hence the distribution $\mathcal{R}f$ satisfies the bound required of the reconstruction of $g$ and we can appeal to uniqueness to conclude.
Now for the broader question.
For a fixed model, the corresponding reconstruction operator is a continuous map from the space of modelled distributions into the space of distributions (or actually even into $C^\alpha$, where $\alpha$ is the lowest degree of a tree in your problem). As a result, since the sequence of modelled distributions $f_n$ produced in the Picard iteration converge to some $f$, the sequence $\mathcal{R}f_n$ of their reconstructions converges in $C^\alpha$ to $\mathcal{R}f$.
The remaining question is why should $\mathcal{R} f$ be related to a solution of your original PDE. Here one of course needs to know that the model you work with (and hence its reconstruction operator) are suitably related to the data of your PDE. In what follows, I assume I work with an admissible model.
Note that $f$ is the solution of an equation on the space of modelled distributions. For example, in the case of $\Phi_3^4$ this equation is of the form
$$f = (K_\gamma + R_\gamma \mathcal{R}) (- f^3 + \Xi)$$
where $\gamma > 0$ and the truncations you mention are part of the definition of the cube on the space of modelled distributions as well as the definitions of $R_\gamma$ and $K_\gamma$. The operators $K_\gamma$ and $R_\gamma$ are defined in such a way that $\mathcal{R} K_\gamma g = K \ast \mathcal{R} g$ and $\mathcal{R} R_\gamma \mathcal{R} g = R \ast \mathcal{R} g$ where $G = K + R$ is the decomposition of the heat kernel into the singular part $K$ and smooth remainder $R$.
As a result, $$\mathcal{R}f = \mathcal{R} (K_\gamma + R_\gamma \mathcal{R}) (- f^3 + \Xi) = G \ast \mathcal{R} (-f^3 + \Xi)$$
so that we will have related $f$ to our original PDE in mild form if we can compute $\mathcal{R}(-f^3 + \Xi)$.
Here we need even more information about the model since the definition of an admissible model doesn't constrain how you define what happens for products. For simplicity, let me first say that the driving noise $\xi$ is smooth and work with the model given by the canonical lift of $\xi$. For this model $\mathcal{R}$ is multiplicative so that
\begin{equation}\label{eq:1}
\mathcal{R}(- f^3 + \Xi) = - (\mathcal{R} f)^3 + \xi
\end{equation}
which tells us that for smooth driving noise all of this recovered a way of solving the original PDE.
Of course, you are actually interested in $\xi$ being space-time white noise which is not smooth. You then proceed by considering the PDE driven by $\xi_\varepsilon = \xi \ast \rho_\varepsilon$ where $\rho_\varepsilon$ is a mollifier at scale $\varepsilon$. It turns out everything in the abstract solution theory was continuous in the driving model so if the canonical lifts of $\xi_\varepsilon$ converged as $\varepsilon \to 0$ then we would be done. This however does not happen and as a result one has to suitably renormalise the canonical lifts of the $\xi_\varepsilon$ to recover convergence of the models as $\varepsilon \to 0$. This then means that $\mathcal{R}_\varepsilon$ (the reconstruction of our "renormalised" model associated to $\xi_\varepsilon$) isn't multiplicative, which has the effect of introducing counterterms in our formula for $\mathcal{R}_\varepsilon(-f^3 + \Xi)$. These counterterms can be computed for suitable choices of renormalisation (see Proposition 9.10 in "A Theory of Regularity Structures" for this computation in the case of $\Phi_3^4$). This answer is already long enough so I won't go more into those details here.
