Here's an answer that leads to explicit formulas whose most complicated term is a sum over a product of binomial coefficients. I haven't yet made any attempt to connect it to the excellent solutions in the answers given earlier.
First some notation. The dihedral group of the square, $D_4,$ has elements $\{e,\rho,\rho^2,\rho^3,h,v,d_1,d_2\},$ where $e$ is the identity, $\rho$ is the $90^\circ$ rotation, $h$ is the reflection about the horizontal axis, $v$ is the reflection about the vertical axis, $d_1$ is the reflection about the forward diagonal, and $d_2$ is the reflection about the backward diagonal.
(Added: Having never learned Pólya enumeration properly, I made my answer overly complicated. I leave my previous answer, but add a more direct route to equation (**) which appears below. In this problem we have $D_4$ acting on the set $X,$ consisting of the $\binom{n^2}{n}$ ways of coloring $n$ small squares in the $n\times n$ grid black, leaving the remaining squares white. The orbit counting theorem leads immediately to (**) since
$$
K(n)=\lvert\mathcal{O}\rvert=\frac{1}{8}\sum_{g\in D_4}C(n,n,\langle g\rangle)
$$
is identical to (**) when one realizes that $\langle\rho\rangle=\langle\rho^3\rangle=C_4.$ The orbit counting theorem is described in many places. For convenience I've appended some background information about it to the end of this post.)
Since the order of $D_4$ is $8,$ all subgroups are of order $1,$ $2,$ $4,$ or $8.$ Either of the two elements of order $4,$ $\rho$ and $\rho_3,$ generates $C_4,$ the rotation group of the square. Each of the five elements of order $2$, $\rho^2,$ $h,$ $v,$ $d_1,$ and $d_2,$ generates a subgroup of order $2.$ There are two additional subgroups of order $4,$ isomorphic to $C_2\times C_2.$ They are $\{e,\rho^2,h,v\}$ and $\{e,\rho^2,d_1,d_2\}.$ Any other group generated by two elements of order $2$ is all of $D_4.$
There are therefore $10$ subgroups with inclusions given by the following Hasse diagram.

Define $C(n,b,G)$ to be the set of configurations of the $n\times n$ grid with $b$ marked squares that are invariant under the group $G.$ We are ultimately interested in $b=n,$ but will start out with slightly greater generality. Define $\overline{C}(n,b,G)$ to be the subset of $C(n,b,G)$ consisting of elements that are not invariant under any group that properly contains $G.$ The number we wish to compute, which we denote $K(n),$ is
$$
\begin{aligned}
(*)\quad K(n)=&\frac{\lvert\overline{C}(n,n,D_4)\rvert}{1}+\frac{\lvert\overline{C}(n,n,\{e,\rho^2,h,v\})\rvert}{2}+\frac{\lvert\overline{C}(n,n,C_4)\rvert}{2}+\frac{\lvert\overline{C}(n,n,\{e,\rho^2,d_1,d_2\})\rvert}{2}\\
&+\frac{\lvert\overline{C}(n,n,\{e,h\})\rvert}{4}+\frac{\lvert\overline{C}(n,n,\{e,v\})\rvert}{4}+\frac{\lvert\overline{C}(n,n,\{e,\rho^2\})\rvert}{4}+\frac{\lvert\overline{C}(n,n,\{e,d_1\})\rvert}{4}\\
&+\frac{\lvert\overline{C}(n,n,\{e,d_2\})\rvert}{4}+\frac{\lvert\overline{C}(n,n,\{e\})\rvert}{8}.
\end{aligned}
$$
(Added: The set $C(n,n,G)$ is the set of colorings whose stabilizer contains $G,$ while $\overline{C}(n,n,G)$ is the set of colorings whose stabilizer equals $G.$ The formula (*) for $K(n),$ the number of orbits, comes from counting each coloring with a weight equal to the reciprocal of the size of its orbit, the weight having been computed using the orbit-stabilizer theorem. (See the end of the post.) For example, if $x$ has stabilizer $D_4,$ then the weight is $1/1=8/8;$ if $x$ has stabilizer $\{e,\rho^2,h,v\},$ then the weight is $1/2=4/8.$ Of course, as noted above, the use of (*) and the inclusion-exclusion argument below is rendered superfluous by the use of the orbit counting theorem. It is interesting that although (*) involves all ten subgroups of $D_4,$ the final expression (**) involves only the seven cyclic subgroups. My naive derivation does not explain this, but the orbit counting theorem makes it manifest that that must be the case.)
From the diagram, we obtain expressions for the $\overline{C}(n,b,G)$ in terms of the $C(n,b,G):$
$$
\begin{aligned}
\lvert\overline{C}(n,b,D_4)\rvert&=\lvert C(n,b,D_4)\rvert\\
\vert\overline{C}(n,b,\{e,\rho^2,h,v\})\rvert&=\lvert C(n,b,\{e,\rho^2,h,v\})\rvert- \lvert C(n,b,D_4)\rvert\\
\vert\overline{C}(n,b,C_4)\rvert&=\lvert C(n,b,C_4)\rvert- \lvert C(n,b,D_4)\rvert\\
\vert\overline{C}(n,b,\{e,\rho^2,d_1,d_2\})\rvert&=\lvert C(n,b,\{e,\rho^2,d_1,d_2\})\rvert- \lvert C(n,b,D_4)\rvert\\
\vert\overline{C}(n,b,\{e,h\})\rvert&=\lvert C(n,b,\{e,h\})\rvert- \lvert C(n,b,\{e,\rho^2,h,v\})\rvert\\
\vert\overline{C}(n,b,\{e,v\})\rvert&=\lvert C(n,b,\{e,v\})\rvert- \lvert C(n,b,\{e,\rho^2,h,v\})\rvert\\
\vert\overline{C}(n,b,\{e,d_1\})\rvert&=\lvert C(n,b,\{e,d_1\})\rvert- \lvert C(n,b,\{e,\rho^2,d_1,d_2\})\rvert\\
\vert\overline{C}(n,b,\{e,d_2\})\rvert&=\lvert C(n,b,\{e,d_2\})\rvert- \lvert C(n,b,\{e,\rho^2,d_1,d_2\})\rvert\\
\vert\overline{C}(n,b,\{e,\rho^2\})\rvert&=\lvert C(n,b,\{e,\rho^2\})\rvert- \lvert \overline{C}(n,b,\{e,\rho^2,h,v\})\rvert- \lvert \overline{C}(n,b,C_4)\rvert\\
&\quad- \lvert \overline{C}(n,b,\{e,\rho^2,d_1,d_2\})\rvert- \lvert \overline{C}(n,b,D_4)\rvert\\
&=\lvert C(n,b,\{e,\rho^2\})\rvert- \lvert C(n,b,\{e,\rho^2,h,v\})\rvert- \lvert C(n,b,C_4)\rvert\\
&\quad- \lvert C(n,b,\{e,\rho^2,d_1,d_2\})\rvert+2 \lvert C(n,b,D_4)\rvert\\
\vert\overline{C}(n,b,\{e\})\rvert&=\lvert C(n,b,\{e\})\rvert- \lvert \overline{C}(n,b,\{e,h\})\rvert- \lvert \overline{C}(n,b,\{e,v\})\rvert- \lvert \overline{C}(n,b,\{e,\rho^2\})\rvert\\
&\quad- \lvert \overline{C}(n,b,\{e,d_1\})\rvert- \lvert \overline{C}(n,b,\{e,d_2\})\rvert- \lvert \overline{C}(n,b,\{e,\rho^2,h,v\})\rvert\\
&\quad- \lvert \overline{C}(n,b,C_4)\rvert- \lvert \overline{C}(n,b,\{e,\rho^2,d_1,d_2\})\rvert- \lvert \overline{C}(n,b,D_4)\rvert\\
&=\lvert C(n,b,\{e\})\rvert- \lvert C(n,b,\{e,h\})\rvert- \lvert C(n,b,\{e,v\})\rvert- \lvert C(n,b,\{e,\rho^2\})\rvert\\
&\quad- \lvert C(n,b,\{e,d_1\})\rvert- \lvert C(n,b,\{e,d_2\})\rvert+2\lvert C(n,b,\{e,\rho^2,h,v\})\rvert\\
&\quad+2\lvert C(n,b,\{e,\rho^2,d_1,d_2\})\rvert
\end{aligned}
$$
These imply that
$$
\begin{aligned}
(**)\quad K(n)=&\frac{\lvert C(n,n,C_4)\rvert}{4}+\frac{\lvert C(n,n,\{e,h\})\rvert}{8}+\frac{\lvert C(n,n,\{e,v\})\rvert}{8}+\frac{\lvert C(n,n,\{e,\rho^2\})\rvert}{8}\\
&+\frac{\lvert C(n,n,\{e,d_1\})\rvert}{8}+\frac{\lvert C(n,n,\{e,d_2\})\rvert}{8}+\frac{\lvert C(n,n,\{e\})\rvert}{8}.
\end{aligned}
$$
It remains to compute the quantities $C(n,b,G)$ for the seven groups that appear in the expression above. (Added: It is usual in Pólya enumeration at this point to use the cycle index and generating function techniques. My more bare-handed approach yields the final expressions in a fairly direct fashion.) We have
$$
C(n,b,\{e\})=\binom{n^2}{b}.
$$
In general, the counting depends on whether $n$ is even or odd. If $n$ is odd, then the action of $C_4$ fixes the center square; every other square has an orbit of size $4.$ Therefore, a configuration that is invariant under the action of $C_4$ is specified by choosing the markings of the center square and of one square from each of the $\frac{n^2-1}{4}$ orbits of size $4.$ (The other three squares in an orbit of size $4$ must be marked the same way as the first square.) Since the total number of marked squares must be $b,$ we have
$$
C(n,b,C_4)=\begin{cases}\binom{(n^2-1)/4}{b/4} & \text{if $b\equiv0\pmod{4}$}\\ \\ \binom{(n^2-1)/4}{(b-1)/4} & \text{if $b\equiv1\pmod{4}$}\\ \\ 0 & \text{otherwise.}\end{cases}
$$
The center square is unmarked when $b\equiv0\pmod{4}$ and marked when $b\equiv1\pmod{4}.$ If $n$ is even then, under the action of $C_4,$ every square has an orbit of size $4.$ Therefore
$$
C(n,b,C_4)=\begin{cases}\binom{n^2/4}{b/4} & \text{if $b\equiv0\pmod{4}$}\\ \\ \binom{n^2/4}{(b-1)/4} & \text{if $b\equiv1\pmod{4}$}\\ \\ 0 & \text{otherwise.}\end{cases}
$$
Specializing to $b=n,$ we get
$$
C(n,n,C_4)=\begin{cases}\binom{n^2/4}{n/4} & \text{if $n\equiv0\pmod{4}$}\\ \\ \binom{(n^2-1)/4}{(n-1)/4} & \text{if $n\equiv1\pmod{4}$}\\ \\ 0 & \text{otherwise.}\end{cases}
$$
The analysis for $G=\{e,\rho^2\}$ is similar. The orbit sizes are now $2$ instead of $4,$ with the exception of the center square in the $n$ odd case, which has orbit size $1.$ The result, for $b=n,$ is
$$
C(n,n,\{e,\rho^2\})=\begin{cases}\binom{n^2/2}{n/2} & \text{if $n$ even}\\ \\ \binom{(n^2-1)/2}{(n-1)/2} & \text{if $n$ odd.}\end{cases}
$$
When $G=\{e,h\}$ and $n$ is odd, the action of $G$ fixes the $n$ squares on the horizontal centerline. All other squares have orbit size $2.$ A configuration invariant under the action of $\{e,h\}$ is specified by marking $j$ of these orbits and $b-2j$ squares on the centerline. When $n$ is even, all squares have orbit size $2,$ so we need only mark $b/2$ of these orbits. The analysis when $G=\{e,v\}$ is similar. In the $b=n$ case, this gives
$$
C(n,n,\{e,h\})=C(n,n,\{e,v\})=\begin{cases}\binom{n^2/2}{n/2} & \text{if $n$ even}\\ \\ \sum_j\binom{(n^2-n)/2}{j}\binom{n}{n-2j} & \text{if $n$ odd.}\end{cases}
$$
When $G=\{e,d_1\}$ or $\{e,d_2\},$ it is a diagonal that is fixed by the action of $G.$ So the analysis is the same as for $G=\{e,h\}$ with $n$ odd. We have
$$
C(n,n,\{e,d_1\})=C(n,n,\{e,d_2\})= \sum_j\binom{(n^2-n)/2}{j}\binom{n}{n-2j}.
$$
Mathematica recognizes this sum, which we denote by $S,$ as a generalized hypergeometric function with argument $-1:$
$$
S=\sum_j\binom{(n^2-n)/2}{j}\binom{n}{n-2j} = \ _3F_2\left[\begin{matrix}1/2-n/2 & -n/2 & (n-n^2)/2\\1/2 & 1 & \end{matrix};-1\right].
$$
Indeed,
$$
\begin{aligned}
S=&\sum_j\binom{(n^2-n)/2}{j}\binom{n}{n-2j} =\sum_j\binom{(n^2-n)/2}{j}\binom{n}{2j}\\
=&\sum_j\frac{\frac{n^2-n}{2}(\frac{n^2-n}{2}-1)\ldots(\frac{n^2-n}{2}-j+1)}{j!}\frac{n(n-1)\ldots(n-2j+1)}{2j(2j-1)\ldots1}\\
=&\sum_j\frac{n-n^2}{2}\left(\frac{n-n^2}{2}+1\right)\ldots\left(\frac{n-n^2}{2}+j-1\right)\frac{\frac{-n}{2}\frac{-n+1}{2}\ldots\frac{-n+2j-1}{2}}{j(j-\frac{1}{2})\ldots\frac{2}{2}\frac{1}{2}}\frac{(-1)^j}{j!}\\
=&\sum_j\frac{n-n^2}{2}\left(\frac{n-n^2}{2}+1\right)\ldots\left(\frac{n-n^2}{2}+j-1\right)\\
&\times\frac{\frac{-n}{2}\left(\frac{-n}{2}+1\right)\ldots\left(\frac{-n}{2}+j-1\right)\frac{1-n}{2}\left(\frac{1-n}{2}+1\right)\ldots\left(\frac{1-n}{2}+j-1\right)}{\frac{1}{2}\frac{3}{2}\ldots\left(\frac{1}{2}+j-1\right)\cdot1\cdot2\ldots\cdot j}\frac{(-1)^j}{j!},
\end{aligned}
$$
which, by definition, is the stated generalized hypergeometric function.
Combining these results, we have
$$
K(n)=\begin{cases}
\frac{1}{8}\binom{n^2}{n} + \frac{1}{4}\binom{n^2/4}{n/4} + \frac{3}{8}\binom{n^2/2}{n/2} +
\frac{1}{4}S & \text{for $n\equiv0\pmod{4}$}\\
\frac{1}{8}\binom{n^2}{n} + \frac{1}{4}\binom{(n^2-1)/4}{(n-1)/4} + \frac{1}{8}\binom{(n^2-1)/2}{(n-1)/2} +
\frac{1}{2}S & \text{for $n\equiv1\pmod{4}$}\\
\frac{1}{8}\binom{n^2}{n} + \frac{3}{8}\binom{n^2/2}{n/2} +
\frac{1}{4}S & \text{for $n\equiv2\pmod{4}$}\\
\frac{1}{8}\binom{n^2}{n} + \frac{1}{8}\binom{(n^2-1)/2}{(n-1)/2} +
\frac{1}{2}S & \text{for $n\equiv3\pmod{4}.$}
\end{cases}
$$
This answer agrees with that of Marko Riedel and the OEIS sequence he cites. (Added: These formulas can be obtained without much difficulty by extracting the relevant coefficient from the final generating function in Marko Riedel's answer.)
Added: (Some background about the orbit-stabilizer theorem and the orbit counting theorem.) Let $G$ be a group acting on a set $X.$. The stabilizer of $x\in X$ is the set of elements of $G$ that fix $x,$ that is, the set of $g\in G$ such that $gx=x.$ The orbit of $x\in X$ is the set of elements of $X$ to which $x$ is sent by some element of $G.$
The orbit-stabilizer theorem states that, for any $x\in X,$
$$
\lvert G\rvert=\lvert\text{orbit}(x)\rvert\cdot\lvert\text{stabilizier}(x)\rvert.
$$
It is a consequence of the fact that $\text{stabilizer}(x)$ is a subgroup of $G$ and the fact that the elements of $\text{orbit}(x)$ are in one-to-one correspondence with the left cosets of $\text{stabilizer}(x).$
The action of $G$ on $X$ partitions $X$ into disjoint orbits. Let $\mathcal{O}$ be the set of orbits. The fixed set of $g\in G$ is the set of elements of $X$ fixed by $g,$ that is, the set of $x\in X$ such that $gx=x.$ The orbit counting theorem, which is often called Burnside's lemma, states that the number of orbits equals the average (over all $g\in G$) of $\lvert\text{fixed}(g)\rvert.$ It is proved by considering the set $P$ of pairs $(g,x)$ satisfying $gx=x.$ We may count the elements of $P$ either by running over $G$ or by running over $X:$
$$
\begin{aligned}
\lvert P\rvert=\sum_{g\in G}\lvert\text{fixed}(g)\rvert&=\sum_{x\in X}\lvert\text{stabilizer}(x)\rvert\\
&=\sum_{x\in X}\frac{\lvert G\rvert}{\lvert\text{orbit}(x)\rvert}\\
&=\lvert G\rvert\sum_{O\in\mathcal{O}}\sum_{x\in O}\frac{1}{\lvert\text{orbit}(x)\rvert}\\
&=\lvert G\rvert\sum_{O\in\mathcal{O}}1\\
&=\lvert G\rvert\cdot\lvert\mathcal{O}\rvert,
\end{aligned}
$$
where the orbit-stabilizer theorem was used to obtain the second line. This implies that
$$
\lvert\mathcal{O}\rvert=\frac{1}{\lvert G\rvert}\sum_{g\in G}\lvert\text{fixed}(g)\rvert,
$$
which is the desired result.