With polynomial method you may prove the same bound even if the lines are allowed to coincide. Then it is sharp: consider the vertical and horizontal lines $x=a$ and $y=a$ taken $a+1$ times for $a=1,2,\ldots,n$. Also it works over any field (and with points $(\alpha_i,\beta_j)$ on the place of $(i,j)$, where $\{\alpha_0,\ldots,\alpha_n\}$ and $\{\beta_0,\ldots,\beta_n\}$ are arbitrary subsets of the field.)
For the proof, we adapt the standard Alon -- Furedi trick for the multisets. Namely, we use the multisets analogue of this formula.
If $A\subset \mathbb{R}$ is a multiset, denote by $n(a)$ the multiplicity of the element $a$ in $A$. If $f(x)\in \mathbb{R}[x]$ is a polynomial, we say ``the values of $f$ on $A$'' for the sequence of numbers $f^{(i)}(a)$, where $a\in A$ and $0\leqslant i\leqslant n(a)-1$. It is well known that if $\deg f\leqslant |A|-1$, $f$ is uniquely determined by its values on $A$ (Hermite interpolation). In particular, the coefficient $[x^{|A|-1}] f(x)$ may be written in the form $\Phi_A(f)$, where $\Phi_A$ is a linear combination of the values of $f$ on $A$.
Remark. Whenever $n(a)=1$, the value $f(a)$ participates in $\Phi_A$ non-trivially. This follows from considering the sample polynomial $f(x):=\prod_{b\in A\setminus a} (x-b)^{n(b)}$.
Now if we have two multisets $A,B\subset \mathbb{R}$ and the polynomial $f(x,y)$ such that $\deg f\leqslant |A|+|B|-2$, we may write the formula $$\left[x^{|A|-1}y^{|B|-1}\right]f=\left(\Phi_A\otimes \Phi_B\right) f(x,y).\quad\quad\quad(1)$$
Here we identify the space of polynomials in $(x,y)$ with the tensor product of polynomials in $x$ and polynomials in $y$. So, (1) allows to reconstruct the coefficient of $x^{|A|-1}y^{|B|-1}$ in the polynomial $f$ via the values $((\partial/\partial x)^i (\partial/\partial y)^j f)(a,b)$, $a\in A$, $b\in B$,
$i$, $j$ are less then corresponding multiplicities of $a$ in $A$ and $b$ in $B$. The coefficient of such guy is a product of corresponding coefficients for 1-dimensional reconstruction formulae.
Once stated, (1) is easy to prove: by linearity, we may think that $f(x,y)=x^k y^l$, $k+l\leqslant |A|+|B|-2$. Then RHS of (1) reads as $\Phi_A(x^k) \Phi_B(y^l)$. If both $k\leqslant |A|-1$, $l\leqslant |B|-1$, we may use the 1-dimensional formulae $\Phi_A(x^k)=[x^{|A|-1}]x^k$, $\Phi_B(y^l)=[y^{|B|-1}]y^l$. If, say, $k\geqslant |A|$, the 1-dimensional formula is no longer valid for $x^k$, but we have $l<|B|-1$, and for $y^l$ the formula $\Phi_B(y^l)=[y^{|B|-1}]y^l$ is valid and gives 0. Cause of this 0, it does not really matter what is the value $\Phi_A(x^k)$.
Well, now back to the problem. Assume that we have less than $n(n+3)$ lines. Take the product of their equations, it is a polynomial $f(x,y)$ of degree less than $n(n+3)$. Consider the multisets $A=B=\{0,1,1,2,2,2,\ldots,n,\ldots,n\}$, $|A|=n(n+3)/2+1$. We have $\deg f<|A|+|B|-2$, thus $[x^{|A|-1}y^{|B|-1}]f=0$. But applying (1) we see that there is unique non-zero term in RHS, which corresponds to $f(0,0)$: all other terms have the form $((\partial/\partial x)^i (\partial/\partial y)^j f)(a,b)$, $(a,b)\ne (0,0)$, $i\leqslant a,j\leqslant b$, they do vanish because of the multiplicity condition. The term $f(0,0)$ goes with non-zero coefficient due to Remark. A contradiction.
You may read more about this multisets formula in [G. Károlyi, Z.-L. Nagy, F. V. Petrov, V. Volkov. A new approach to constant term identities and Selberg-type integrals. Adv. Math. 277 (2015), 252-282], earlier the tensor product view appeared in [U. Schauz. Algebraically solvable problems: describing polynomials as equivalent to explicit solutions. Electron. J. Combin. 15 (2008) #R10, 35 p.]
When the lines are all distinct, this seems to be too weak bound, let me prove $\Omega(n^{5/2})$.
If a point is covered by a line, draw a flower at this place. The number of flowers should be at least $cn^3$ for some $c>0$. On the other hand, denote by $r_k$ the number of lines containing at least $k$ points. The vector between consecutive lattice points on such line may be chosen by $O(n^2/k^2)$ ways, the line should pass through a point which may be chosen by $O(n^2)$ ways, and each such line is counted at least $k$ times, totally we get $r_k\leqslant С n^4/k^3$ for certain $C$. Choose large $M$ and consider only $k\geqslant M\sqrt{n}$. We get $\sum_{k\geqslant M\sqrt{n}} r_k\leqslant Cn^4/(M^2n)=(C/M^2)n^3<cn^3/2$ if $M$ is large enough. This means that if label only $\max(s,[M\sqrt{n}])$ flowers corresponding to each line containing exactly $s$ points, the total number $\sum_{k> M\sqrt{n}} r_k$ of not-labelled flowers is at most $cn^3/2$. All other $cn^3/2$ flowers are labelled, and each line contains at most $M\sqrt{n}$ labelled flowers. Thus the number of lines is not less then $\frac{c}{2M} n^{5/2}$.