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 bounded.
(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 with $p$ is a bounded linear operator on $E_k$ for each $k$. Hence, by standard perturbation theory for $C_0$-semigroups, the differential operator on the right hand side of the PDE generates a $C_0$-semigroup on $E_k$. Applying this to $k=0$ yields assertion (i).
Moreover, these $C_0$-semigroups are consistent on the spaces $E_k$ (this follows, for instance, from the Dyson-Phillips series and the fact that each space $E_{k+1}$ embedds continuously into $E_k$).
If $f_0$ decay faster then 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$. Thus, for each time $t$, the function $T_0(t)f_0$ decays faster in space than any polynomial. $\square$