I do not know of specific references, but perhaps you could proceed as follows.
Lemma: Suppose that $p$ is an odd prime and $f: \mathbb{Z}/p\mathbb{Z}\rightarrow \mathbb{F}_2^n$ is injective. Then
\begin{equation*}
\mathbb{P}(f(x+y)=f(x)+f(y))=p^{-\Omega(1)}.
\end{equation*}
Proof: Write $\Gamma:=\{(x,f(x)):x \in \mathbb{Z}/p\mathbb{Z}\} \subset \mathbb{Z}/p\mathbb{Z} \times \mathbb{F}_2^n$ and $\varepsilon$ for the probability we are interested in. Using the language of summation instead of expectation the says
\begin{equation*}
\langle 1_\Gamma\ast 1_\Gamma,1_\Gamma\rangle_{\ell_2} = \varepsilon p^2,
\end{equation*}
where
\begin{equation*}
g\ast h(x):=\sum_{y+z=x}{g(y)h(z)} \text{ and }\langle g,h\rangle_{\ell_2}:=\sum_x{g(x)\overline{h(x)}},
\end{equation*}
for functions $g$ and $h$.
Since $f$ is a function, $|\Gamma|=p$, and so by Cauchy-Schwarz in the usual way we get
\begin{equation*}
\|1_\Gamma\ast 1_\Gamma\|_{\ell_2}^2 \geq \frac{1}{|\Gamma|}\langle 1_\Gamma\ast 1_\Gamma,1_\Gamma\rangle_{\ell_2}^2 \geq \varepsilon^2 \frac{p^4}{|\Gamma|} = \varepsilon^2 |\Gamma|^3.
\end{equation*}
The Balog-Szemerédi-Gowers lemma then tells us that there is a set $\Gamma' \subset \Gamma$ with $|\Gamma'| \geq \varepsilon^{O(1)}|\Gamma|$ and
\begin{equation*}
|\Gamma'+\Gamma'| \leq \varepsilon^{-O(1)}|\Gamma'|.
\end{equation*}
For a set $\Lambda \subset \mathbb{Z}/p\mathbb{Z} \times \mathbb{F}_2^n$ write $2\cdot \Lambda:=\{(2x,2y): (x,y) \in \Lambda\}$. Since $p \neq 2$ the map
\begin{equation*}
\Gamma' \rightarrow 2\cdot \Gamma'; (x,f(x)) \mapsto (2x,2f(x))
\end{equation*}
is a bijection and hence $|2\cdot \Gamma'| = |\Gamma'|=\varepsilon^{O(1)}p$. Now, by Cauchy-Schwarz again we have
\begin{equation}
\tag{1}\label{eq:1}
\|1_{\Gamma'} \ast 1_{2\cdot \Gamma'}\|_{\ell_2}^2 \geq \frac{(|\Gamma'||2\cdot \Gamma'|)^2}{|\Gamma'+2\cdot\Gamma'|} = \varepsilon^{O(1)}p^3,
\end{equation}
where the last inequality follows from the fact that
\begin{equation*}
|\Gamma'+2\cdot \Gamma'| \leq |\Gamma' + \Gamma' + \Gamma'| \leq \varepsilon^{-O(1)}|\Gamma'| = \varepsilon^{-O(1)}p
\end{equation*}
by Plünnecke's inequality.
On the other hand, if $(a,f(a)),(b,f(b)) \in \Gamma'$ and $(2c,2f(c)),(2d,2f(d)) \in 2\cdot\Gamma'$ have
\begin{equation}
\tag{2}\label{eq:2}
(a,f(a))+(2c,2f(c))=(b,f(b)) +(2d,2f(d)),
\end{equation}
then $a+2c=b+2d$ and $f(a)+2f(c)=f(b)+2f(d)$. But the image of $f$ is in $\mathbb{F}_2^n$ and so $2f(c)=0=2f(d)$ whence $f(a)=f(b)$. On the other hand $f$ is injective so $a=b$. It follows from this (and the equation $a+2c=b+2d$) that $2c=2d$, but $p$ is odd and so $c=d$. Thus there are at most $p^2$ quadruples such that \eqref{eq:2} holds. These quadruples are precisely what is counted by the left hand side of \eqref{eq:1} and so
\begin{equation*}
\varepsilon^{O(1)}p^3 \leq p^2,
\end{equation*}
from which we get $\varepsilon = p^{-\Omega(1)}$ as claimed.
