I believe that I have a soft analysis argument that nonconstructively addresses my question in all dimensions. The idea of the argument is to use $H^{1}$-$BMO$ duality together with the fact that the space of essentially bounded functions is a proper subspace of $BMO$ to show that $H^{1}(\mathbb{R}^{n})$ is a proper subspace of $L^{1}(\mathbb{R}^{n})$ functions with vanishing integral.
In the interest of definiteness, we will take the $H^{1}$-norm $\|\cdot\|_{H^{1}}$ to be the atomic norm:
$$\|f\|_{H^{1}}:=\inf\left\{\sum_{k}|\lambda_{k}| : f=\sum_{k=1}^{\infty}\lambda_{k}a_{k}\right\},$$
where the infimum is taken over all decompositions $f=\sum_{k}\lambda_{k}a_{k}$ as a sum of $(1,\infty)$-atoms $a_{k}$.
Denote the subspace of $L^{1}(\mathbb{R}^{n})$ consisting of functions with mean zero by $L_{0}^{1}(\mathbb{R}^{n})$. By endowing this subspace with the $L^{1}$-norm, we obtain a Banach space. Using the basic result, $(L^{1})^{*}\cong L^{\infty}$, we obtain the analogous statement for $(L_{0}^{1})^{*}$.
Lemma 1. Let $\sim$ denote the equivalence relation on measurable functions defined by $f\sim g$ if and only if $f-g$ is almost everywhere (a.e.) equal to a constant. Then $(L_{0}^{1})^{*}\cong L^{\infty}/\sim$. More precisely, for every every continuous linear functional $L$ on $L_{0}^{1}$, there exists a unique equivalence class $[g]_{\sim}\in L^{\infty}/\sim$, such that
$$L(f)=\int_{\mathbb{R}^{n}}f(x)g(x)dx,\quad\forall f\in L_{0}^{1}(\mathbb{R}^{n})$$
and $\|L\|_{*}=\|g\|_{L^{\infty}}$.
Proof. It is evident that $L^{\infty}/\sim$ embeds in $(L_{0}^{1})^{*}$. Let $L\in (L_{0}^{1})^{*}$. Fix a nonnegative $C_{c}^{\infty}$ function $\eta$ with $\int\eta=1$, and define a linear functional $L_{\eta}$ on $L^{1}$ by
$$L_{\eta}(f):=L(f-(\int f)\eta),\quad f\in L^{1}(\mathbb{R}^{n})$$
I claim that $L_{\eta}\in (L^{1})^{*}$. Indeed,
$$\|L_{\eta}(f)\|=\|L(f-(\int f)\eta)\|\leq\|f-(\int f)\eta\|_{L^{1}}\leq\|f\|_{L^{1}}+\|f\|_{L^{1}}\|\eta\|_{L^{1}}=2\|f\|_{L^{1}}$$
Thus, by $L^{p}$-$L^{p'}$ duality, there exists a $g\in L^{\infty}$ such that
$$L_{\eta}(f)=\int_{\mathbb{R}^{n}}f(x)g(x)dx,\quad\forall f\in L^{1}(\mathbb{R}^{n})$$
If we instead choose a different $\gamma\in C_{c}^{\infty}$, nonnegative, and with $\int\gamma=1$, then $\eta-\gamma\in L_{0}^{1}$ and we have that there exists an $h\in L^{\infty}$ such that
\begin{align*}
\int_{\mathbb{R}^{n}}f(x)h(x)dx&=L_{\gamma}(f)=L\left((f-(\int f)\eta+(\int f)(\eta-\gamma)\right)\\
&=L_{\eta}(f)+\underbrace{L(\eta-\gamma)}_{:=c}\int f\\
&=\int_{\mathbb{R}^{n}}f(x)[g(x)+c]dx
\end{align*}
From Lebesgue differentiation theorem, we conclude that $f-h=c$ a.e. $\Box$
Lemma 2. Let $f\in L_{loc}^{2}(\mathbb{R}^{n})$ be supported in a cube $Q$, and suppose that
$$\int_{Q}f(x)g(x)dx=0$$
for all compactly supported $g\in L_{0}^{2}(\mathbb{R}^{n})$ (square integrable zero mean functions). Then there is a constant $c\in\mathbb{C}$ such that $f=c$ a.e. in $Q$.
Proof. Let $\eta\in C_{c}^{\infty}$ be supported in $Q$, nonnegative, and satisfy $\int\eta=1$. Set $\eta_{\epsilon}:=\epsilon^{-n}\eta(\cdot/\epsilon)$. For a.e. $x\in Q^{o}$, we have that
$$f(x)=\lim_{\epsilon\rightarrow 0}\int_{Q}f(y)\eta_{\epsilon}(x-y)dy$$
For such an $x$ and $\epsilon>0$ sufficiently small, we have that the function $\varphi_{\epsilon}(y):=\eta_{\epsilon}(x-y)-|Q|^{-1}\chi_{Q}(y)$ is compactly supported and in $L_{0}^{2}$. Whence,
$$\int_{Q}f(y)\varphi_{\epsilon}(y)dy=0\Longrightarrow\int_{Q}f(y)\eta_{\epsilon}(x-y)dy=|Q|^{-1}\int_{Q}f(y),\quad\forall 0<\epsilon\ll 1$$
Letting $\epsilon\rightarrow 0$, we obtain that
$$f(x)=|Q|^{-1}\int_{Q}f(y)dy,$$
and we see that this equality holds for a.e. $x\in Q$. $\Box$
Assume for the sake of a contradiction that $L_{0}^{1}(\mathbb{R}^{n})=H^{1}(\mathbb{R}^{n})$. Recall that we have the embedding
$$(H^{1}(\mathbb{R}^{n}),\|\cdot\|_{H^{1}})\hookrightarrow (L_{0}^{1}(\mathbb{R}^{n}),\|\cdot\|_{L^{1}}),\quad \left\|f\right\|_{L^{1}}\leq\left\|f\right\|_{H^{1}}$$
By hypothesis, this embedding is in fact bijective. Whence by the Bounded Inverse Theorem, the inverse map is also bounded, showing that $\|\cdot\|_{H^{1}}$ and $\|\cdot\|_{L^{1}}$ are equivalent norms. Now recall that $BMO(\mathbb{R}^{n})$ is the dual of $H^{1}(\mathbb{R}^{n})$ in the sense that for any $L\in (H^{1})^{*}$, there exists a BMO representative function $g$ (whose equivalence class is unique), such that
$$L(f)=\int_{\mathbb{R}^{n}}f(x)g(x)dx$$
for all compactly supported functions in $L_{0}^{2}(\mathbb{R}^{n})$ and $\left\|L\right\|_{*}\sim\|g\|_{BMO}$. Conversely, any BMO function extends to a bounded linear functional on $H^{1}(\mathbb{R}^{n})$ by being initially defined by integration against this subspace of $H^{1}$. Since $L$ also defines a continuous linear functional on $L_{0}^{1}$, by the first lemma, there exists an $h\in L^{\infty}$ such that
$$L(f)=\int_{\mathbb{R}^{n}}f(x)h(x)dx,\quad\forall f\in L_{0}^{1}(\mathbb{R}^{n})$$
Since BMO functions are locally square integrable as a consequence of the John-Nirenberg inequality and $L^{\infty}$ functions are trivially locally square integrable, we can apply Lemma 2 to $g-h$ to conclude that $g-h=c$ is a.e. equal to a constant in a cube $Q$. By taking a sequence of increasing cubes, we obtain that $g-h=c\in\mathbb{C}$ a.e. in $\mathbb{R}^{n}$. But if we take $g(x):=\log|x|\in BMO$, we obtain a contradiction, as $g$ is clearly not essentially bounded.
