I assume by square lattice you mean a 3D cubic lattice because there's no translational symmetry along the $z$-axis for a 2D square lattice.
Suppose the masses are located at $(n_x, n_y, n_z)$ where $n_{x,y,z}\in\mathbb Z$. Let's also define the unit of mass and length so that $m=l=1$.
Consider the total force acted on the mass point at (0, 0, 0) just due to the 1st octant $(x>0,y>0,z>0)$:
$$\begin{aligned}
\mathbf F_{+++} &= -G \sum_{n_x=1}^\infty \sum_{n_y=1}^\infty \sum_{n_z=1}^\infty
\frac{n_x \hat{\mathbf x} + n_y \hat{\mathbf y} + n_z \hat{\mathbf z}}{(n_x^2+n_y^2+n_z^2)^{3/2}} \\
&= -G \left( \hat{\mathbf x} \sum_{n_x=1}^\infty \sum_{n_y=1}^\infty \sum_{n_z=1}^\infty \frac{n_x}{(n_x^2+n_y^2+n_z^2)^{3/2}} + \dotsb \right),
\end{aligned}$$
however, the sum actually diverges, since,
$$\begin{aligned}\sum_{n_x=1}^\infty \sum_{n_y=1}^\infty \sum_{n_z=1}^\infty \frac{n_x}{(n_x^2+n_y^2+n_z^2)^{3/2}} &\ge \int_1^\infty \int_1^\infty \int_1^\infty \frac{n_x}{(n_x^2+n_y^2+n_z^2)^{3/2}} dn_x dn_y dn_z \\
&= \int_1^\infty \int_1^\infty \frac1{\sqrt{1+n_y^2+n_z^2}} dn_y dn_z \\ &= \infty, \end{aligned} $$
so while symmetry suggests that the force at center is 0, mathematically it is not well defined.
Of course, if we assume the net force can be well-defined as 0 (e.g. the gravity actually decays faster than $1/r^3$!), then Points 1 and 3 are correct. When we remove a particle from the lattice, the contribution $-\frac{GMm\hat{\mathbf r}}{r^\alpha}$ will be subtracted from it, so it is as if there is a particle of negative mass $-m$ put to that point. This is because force is additive and gravity is proportional to mass.
(Yeah this is stating the obvious.)
