I know of a quirky but loveable proof that I'll now present, but I'm not sure whether it is in any book. I learned it from Jay Rosen, and variants of it are in his book with Marcus on local times, but I don't think you'll find exactly this there.
We'll apply Kolmogorov's criterion, so we need to show that $E[|L_a(t) -L_{a'}(t)|^n] \leq C |a-a'|^{1+m}$, and then we get Holder continuity of any order less than $m/n$. We'll use the formal identity $\delta_a (x) = \frac{1}{2\pi} \int_{-\infty}^\infty e^{i(x-a)p}dp$; this can all be justified rigorously by choosing an approximate delta and taking a limit.
Now, we calculate n-th moments by multiplying a bunch of integrals (with different dummy variables) together.
\begin{equation*}
\begin{split}
E[|L_a(t) -L_{a'}(t)|^n] & \leq C \int_{\mathbb{R}^n} \int_{[0,t]^n} \prod_{j=1}^n |e^{-iap_j} - e^{-ia'p_j}| E[\prod_{j=1}^n e^{ip_j B_{s_j}}] \prod_{j=1}^n ds_j \prod_{j=1}^n dp_j.
\end{split}
\end{equation*}
Now we'll use the bound $|e^{-iap_j} - e^{-ia'p_j}| \leq C |p_j|^\gamma |a-a'|^\gamma$, valid for any $\gamma \in [0,1]$, however in our case we'll require $\gamma < 1/2$ (the reason for this will be clear later). To analyze the term $E[\prod_{j=1}^n e^{ip_j B_{s_j}}]$, let us assume that $s_n \leq s_{n-1} \leq \ldots \leq s_1$; the contribution from any other configuration will be equal and this can all be swallowed into the constant $C$. We can then write $$\prod_{j=1}^n e^{ip_j B_{s_j}} = e^{ip_1 (B_{s_1}-B_{s_2})}e^{i(p_1+p_2) (B_{s_2}-B_{s_3})}e^{i(p_1+p_2+p_3) (B_{s_3}-B_{s_4})} \ldots e^{i(p_1+p_2+\ldots + p_n) B_{s_n}}$$ and the expectation will factor due to the independence of increments. If we let $u_j = p_1 + \ldots + p_j$ and $v_j = s_j-s_{j+1}$ ($v_n = s_n$), we get
$$E[\prod_{j=1}^n e^{ip_j B_{s_j}}] = \prod_{j=1}^n e^{-u_j^2 v_j}.$$
Putting all this together, we get
\begin{equation*}
\begin{split}
E[|L_a(t) -L_{a'}(t)|^n] & \leq C|a-a'|^{n\gamma} \int_{\mathbb{R}^n} \int_{[0,t]^n} \prod_{j=1}^n |p_j|^\gamma \prod_{j=1}^n e^{-u_j^2 v_j} \prod_{j=1}^n ds_j \prod_{j=1}^n dp_j.
\end{split}
\end{equation*}
In fact, the $ds$ integrals are really over the region $\sum_j v_j \leq t$, but we can bound it by changing the variables taking the $s_j$'s to the $v_j$'s and then integrating over $[0,t]^n$, using the bound $\int_0^t e^{-u_j^2 v_j}dv_j \leq \frac{C}{1+u_j^2}$. We can also write $|p_j| = |u_j - u_{j-1}| \leq |u_j|+|u_{j-1}|$, and do another change of variables, taking the $p_j$'s to $u_j$'s, and we obtain
\begin{equation*}
\begin{split}
E[|L_a(t) -L_{a'}(t)|^n] & \leq C|a-a'|^{n\gamma} \int_{\mathbb{R}^n} \prod_{j=1}^n \frac{(|u_j|+|u_{j-1}|)^\gamma}{1+u_j^2} du_j.
\end{split}
\end{equation*}
Expanding the numerator, we get a sum of many terms of the form
$$
C|a-a'|^{n\gamma} \int_{\mathbb{R}^n} \prod_{j=1}^n \frac{|u_j|^{r_j\gamma}}{1+u_j^2} du_j
$$
where $r_j \in \{0,1,2\}$. Since $\gamma < 1/2$ these integrals are all finite, and the result follows.
