Edit on 10/13/2016
The number of similarity classes $|I|$ of $\mathrm{GL}_3(\mathbb{F}_p)$ below easily follows from this answer. In fact, without having to consider $8$ different types, we have
$$
|I|=\sum_{\lambda \in P_3} p^{\ell(\lambda)}
$$
where $P_3$ is the set of partition of $3$, and $\ell(\lambda)$ is the number of parts of $\lambda$.
Then
$$
|I|=p^3+p^2+p.
$$
An Upper bound for $p\geq 5$
Let $A\in M_3(\mathbb{Z}/p\mathbb{Z})$. Denote by $\mathbb{F}_p=\mathbb{Z}/p\mathbb{Z} $ and $\mathbb{F}_p[x]$ the polynomial ring over $\mathbb{F}_p$. We write $M^A=\mathbb{F}_p^3$ the finitely generated module over $\mathbb{F}_p[x]$ where $x \cdot v = A \cdot v$ for any $v\in\mathbb{F}_p^3$. It is well-known that the dimension of the space of commuting matrices $C_A = \{ B \in M_3(\mathbb{Z}/p\mathbb{Z}) | AB=BA\}$ depends only on the similarity class of $A$. The similarity class of $A$ falls into one of the following eight types:
$aaa$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)}$ where $a\in\mathbb{F}_p$.
$aa^2$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-a)^2}$ where $a\in\mathbb{F}_p$.
$a^3$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)^3}$ where $a\in\mathbb{F}_p$.
$C$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(C(x))}$ where $C(x)\in\mathbb{F}_p[x]$ is a monic cubic irreducible polynomial over $\mathbb{F}_p$.
$abb$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}$ where $a, b\in\mathbb{F}_p$ and $a\neq b$.
$ab^2$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)^2}$ where $a, b\in\mathbb{F}_p$ and $a\neq b$.
$aB$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(B(x))}$ where $a \in\mathbb{F}_p$ and $B(x)$ is a monic quadratic irreducible polynomial over $\mathbb{F}_p$.
$abc$-type: $M^A \simeq \frac{\mathbb{F}_p[x]}{(x-a)}\oplus \frac{\mathbb{F}_p[x]}{(x-b)}\oplus \frac{\mathbb{F}_p[x]}{(x-c)}$ where $a,b,c\in\mathbb{F}_p$ are distinct.
Let $\{I_1,I_2,\ldots, I_8\}$ be the set of the above eight types similarity classes, and let $I=\cup_{j=1}^8 I_j$ be the set of all distinct similarity classes. Define $G=GL_3(\mathbb{F}_p)$, and $G_A=\{H\in G| AH=HA\}$. By the orbit-stabilizer formula,
$$
|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|=\sum_{A\in I} |G/G_A| |C_A|\leq \sum_{A\in I} |G| \frac{|C_A|}{|C_A|-3p^{\dim C_A -1}}.
$$
For the last inequality, let $A\in I_j$ and $\dim C_A = m$. Let $\{I=B_1, \ldots, B_m\}$ be a basis for $C_A$. Consider
$$
\det (c_1 B_1 + \cdots + c_m B_m) = 0. \ \ \ (*)
$$
For any given $c_2,\ldots, c_m$, the equation $(*)$ is a cubic equation in $c_1$. Therefore, the number of solution to $(*)$ is at most $3p^{m-1}$.
We remark that with a lot of effort, we will be able to find explicit formula for $|G_A|$. This is certainly more difficult than finding $|C_A|$.
Note that
$$
1\leq\frac{|C_A|}{|G_A|}\leq\frac{|C_A|}{|C_A|-3p^{\dim C_A -1}}=\frac1{1-\frac3p}=\frac p{p-3}.
$$
The number of elements in eacy $I_j$ can be easily found:
$$|I_1|= |I_2|= |I_3|=p, \ |I_4|= \frac{p^3-p}3, \ |I_5| =|I_6|=p(p-1), \ |I_7|=\frac{p(p^2-p)}2,$$
$$|I_8|=\frac{p(p-1)(p-2)}6.$$
Therefore,
$$
|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|\leq \frac p{p-3}|G|\left(3p+\frac{p^3-p}3 +2p(p-1)+\frac{p(p^2-p)}2+\frac{p(p-1)(p-2)}6\right)
$$
$$
=\frac p{p-3} |G| (p^3+p^2+p).
$$
Combining these, we obtain that
$$
\frac{ |G|(p^3+p^2+p)}{(p^9)^2}\leq \frac{|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|}{(p^9)^2}\leq \frac p{p-3}\frac{|G|(p^3+p^3+p)}{(p^9)^2}.
$$
Since $|G|=(p^3-1)(p^3-p)(p^3-p^2)$, we have as $p\rightarrow\infty$,
$$
\frac{|\{(A,B)\in M_3(\mathbb{F}_p)^2 | AB=BA\}|}{(p^9)^2}\sim \frac1{p^6}.
$$