Here is an idea of a partial answer; cf. the answer by Christian Remling. Assume that $(1)$ holds for all $\theta$ in (say) a right neighborhood of $0$ (including $0$), where $|f(x)|=O(1/(1+|x|^k))$ for each $k>0$ (note that the latter condition is not satisfied by the function $f=I_{(0,\infty)}-c$ in the example in the question). Then one can differentiate the integral in $(1)$ in $\theta$ to the right of $0$ to get \begin{equation} \int_{\mathbb R}\frac{\partial^k g_\theta(x)}{\partial\theta^k}\Big|_{\theta=0}\,f(x)\,dx=0, \tag{2} \end{equation}\begin{equation} \left. \int_{\mathbb R}\frac{\partial^k g_\theta(x)}{\partial\theta^k} \right|_{\theta=0}\,f(x)\,dx=0, \tag{2} \end{equation} for all $k=0,1,\dots$. Expanding $g_\theta(x)/g_0(x)=e^{(\theta x)-(\theta x)^2/2}$ into the Maclaurin series, we have \begin{equation*} \frac{\partial^k g_\theta(x)}{\partial\theta^k}\Big|_{\theta=0}=a_k x^k\,e^{-1/2}, \end{equation*}\begin{equation*} \left. \frac{\partial^k g_\theta(x)}{\partial\theta^k} \right|_{\theta=0}=a_k x^k\,e^{-1/2}, \end{equation*} where \begin{equation*} a_k:=\left(-\frac{1}{2}\right)^k k! \sum _{n=\left\lceil k/2\right\rceil }^k \frac1{n!}\,\binom{n}{k-n} (-2)^n. \end{equation*} The sequence $(a_k)$ is the sequence A001464 at http://oeis.org/A001464.
It appears that $a_k\ne0$ for all $k=0,1,\dots$ except $k=2$. If that is indeed so, then, by $(2)$, we would have \begin{equation} \int_{\mathbb R}x^kf(x)\,dx=0 \tag{3} \end{equation} for all $k=0,1,\dots$ except $k=2$. Letting $y:=x^4$, we see that $x^2=\sqrt y$ can be approximated arbitrarily closely by polynomials in $y=x^4$. So, one would have $(3)$ for all $k=0,1,\dots$, including $k=2$, if $|f(x)|$ decreases fast enough as $|x|\to\infty$. This would imply that $f=0$ almost everywhere.