For orientable competitors (integral currents), Federer proved a much more general result in
MR0388226 Reviewed Federer, Herbert A minimizing property of extremal submanifolds. Arch. Rational Mech. Anal. 59 (1975), no. 3, 207–217. (Reviewer: William K. Allard) 49F20
His result concerns more general elliptic integrands and proves the area-minimality of any region $M$ of a small enough area.
For the area-integrand, a sketch of proof for orientable and non-orientable competitors is as follows. The proof is essentially due to Robert Bryant and Gary Lawlor.
Observation 1: for a map $F:W\to V,$ with $W$ a region in $\tilde{M}$ and $V$ a region of $M,$ $F$ is area-non-increasing if and only if $F^\ast dvol_M$ is a calibration. Here $dvol_M$ is any choice of volume form locally if $M$ is non-orientable.
Proof: straightforward calculation.
Observation 2: any multiplicity $1$ integral or mod $p$ ($p\ge 2$ integers) current $T$ that admits an area-non-increasing map $F$ onto it is area-minimizing. Here we require $F$ restricted to the support of $T$ to be identity.
Proof: straightforward calculation.
Remark: mod $2$ currents include unorientable submanifolds and mod $p$ currents for $p\ge 3$ can admit more exotic singularities like triple-junctions, etc.
Thus, to prove the conclusion you want, it suffices to find a suitable map $F$ onto a neighborhood of $p$ in $M,$ so that $F^\ast dvol_M$ is a calibration.
Unfortunately, as pointed out by Robert Bryant, the usual normal bundle projection $\pi$ does not work. Straightforward calculations in Fermi coordinates show that $\pi^\ast dvol_M$ has comass $1+O(y^2),$ with $y$ the distance to $M.$
It is the idea of Robert Bryant and Gary Lawlor to overcome this by tweaking the projection to dominate the $y^2$ term.
One possible route is as follows. Set up a normal coordinate $(x_1,\cdots, x_m)$ centered at $p$ on $M.$ Take orthonormal frames $\nu_1,\cdots,\nu_m$ in the normal bundle of $M$ near $p,$ so that they are parallel with respect to the intrinsic geodesics on $M$ starting at $p$. Set up a Fermi coordinate using this frame, with coordinate label $(x_1,\cdots,x_m,y_1,\cdots,y_n).$ Now consider a coordinate transformation $x_j=(1+C\sum_iy_i^2)u_j,v_j=y_j$, with $C>0$ a constant. It is straightforward to verify that $(u_j,v_j)$ serves as a new coordinate system.
Define $F(u_1,\cdots,u_m,v_1,\cdots,v_n)=(u_1,\cdots,u_m,0,\cdots,0)$. Straightforward calculations back in $(x,y)$ coordinate show that $F^\ast dvol_M$ is a calibration that calibrates $M$ near $p$ for $C$ large enough.
Remark 1: At every step, we need to shrink our coordinate chart a bit. The transformation introduces a $(1+Cy^2)^{-m}$ factor and some other $1+O(y^2)$ factors in the comass, while not introducing any new $1+cy$ factor. A tedious calculation can verify that the $(1+Cy^2)^{-m}$ eventually dominates. See Tubes by Alfred Gray on how to do the calculations.
Remark 2: One might be curious about where is the minimality used. It is hidden in the expansion of comass $\pi^\ast\omega=1+O(y^2).$ If $M$ is not minimal, then the comass will be $1+O(y).$ The above construction breaks down in this case.