The anwer is yes, even if $p: \mathbb{R}^n \to \mathbb{R}$ is only a bounded continuous function.
Preliminary observations:
Let us first consider the case $p=0$.
Since $q$ is globally Lipschitz on $\mathbb{R}^n$, say with Lipschitz constant $L$, all solutions of the ODE $\dot x = q(x)$ on $\mathbb{R}^n$ exist globally and, if $(\varphi_t)_{t \in \mathbb{R}}$ denotes the induced flow in $\mathbb{R}^n$, each function $\varphi_t: \mathbb{R}^n \to \mathbb{R}^n$ is Lipschitz continuous with constant $e^{|t|L}$ (this follows from Grönwall's lemma).
Since we assumed $p=0$, the solution of your PDE is given by $$ f(t,x) = f_0(\varphi_t(x)) \tag{*}\label{1} $$ for each $t$ and $x$.
Hence, if $f_0$ decays faster than every polynomial, so does $f(t, \cdot)$ for each fixed time $t$.
Let us now phrase this in a more functional analytic language: For each integer $k \ge 0$, consinder the weight function $w_k: \mathbb{R}^n \to [1,\infty)$ given by $w_k(x) =1 + \lvert x \rvert^k$, and consider the Banach space $E_k$ of all continuous functions $f: \mathbb{R}^n \to \mathbb{R}$ which satisfies $w_k f \in C_0(\mathbb{R}^n; \mathbb{R})$; we endow the space $E_k$ with the weighted supremum norm $\|\cdot\|_{w_k}$ given by $$ \|f\|_{w_k} = \|w_k f\|_\infty. $$
It follows from $\eqref{1}$ that the Koopman group $(T_k(t))_{t \in \mathbb{R}}$ given by $$ T_k(t)f = f \circ \varphi_t $$ for each $f \in E_k$ is a well-defined one-parameter group of bounded linear operators on $E_k$. Moreover, the orbits of the group are continuous with respect to the weak topology on $E_k$ (use that (i) $E_k$ is isomorphic to $C_0(\mathbb{R}^n; \mathbb{R})$ via multiplication with $w_k$, that (ii) for each $x \in \mathbb{R}^n$, the mapping $t \mapsto \varphi_t(x)$ is continuous, and (iii) the Riesz representation theorem for measures on $C_0(\mathbb{R}^n; \mathbb{R})$). Hence, by a standard result from $C_0$-semigroup theory, the the group $(T_k(t))_{t \in \mathbb{R}}$ is even continuous with respect to the strong operator topology, i.e., it is a $C_0$-group on $E_k$.
It is trivial, but important to observe, that the groups $(T_k(t))_{t \in \mathbb{R}}$ act consistently on the spaces $E_k$.
Using these preliminary observations, we can show:
Theorem. Let $p$ be continuous and boundeddecay faster than every polynomial.
(i) For each $f_0 \in C_0(\mathbb{R}^n; \mathbb{R})$ there exists a unique mild solution $f: \mathbb{R} \to C_0(\mathbb{R}^n; \mathbb{R})$ of the Cauchy problem $$ \dot f = \langle q, \nabla \rangle f(t) + p, \qquad f(0) = f_0. $$
(ii) If $f_0$ decays faster than every polynomial and $f: \mathbb{R} \to C_0(\mathbb{R}^n; \mathbb{R})$ denotes the mild solution from (i), then for each fixed time $t \in \mathbb{R}$, the function $f(t)$ also decays faster (in space) than every polynomial.
Proof. We note that multiplication withFix $p$$k \ge 0$. Since $(T_k(t))$ is a bounded linear operator$C_0$-group on $E_k$ for eachand since $k$. Hence$p \in E_k$, by standard perturbation theory forit follows that, if $C_0$-semigroups$f_0 \in E_k$, the differential operator on the right hand side of the PDE generateshas a unique mild solution $C_0$-semigroup on$f: [0,\infty) \to E_k$ given by $$ f(t) = T_k(t)f_0 + \int_0^t T_k(s) p \, ds, $$ where the integral is a Bochner integral in $E_k$. Applying this to $k=0$ yields assertion, we obtain (i).
Moreover, these $C_0$-semigroupssolutions are consistent on the spaces $E_k$-spaces (this followsi.e., for instanceif $f_0 \in E_{k+1} \subseteq E_k$, fromthen the Dyson-Phillips seriesmild solutions in $E_k$ and the fact that each space $E_{k+1}$ embedds continuously into $E_k$are the same).
If $f_0$ decaydecays faster thenthan every polynomial, then we have $f_0 \in E_k$ for each $k$. Hence, for each time $t$ and each $k$, we have $T_0(t)f_0 = T_k(t)f \in E_k$$f(t) \in E_k$. Thus, for each time $t$, the function $T_0(t)f_0$$f(t)$ decays faster in space than any polynomial. $\square$