Yes.

Consider the adjacency matrices
$$ A = \left[\begin{array}{rrrrrrrrrrr}
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\\\
0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
\end{array}\right] $$
and
$$ B = \left[ \begin{array}{rrrrrrrrrrr}
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\\\
0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0
\end{array}\right]. $$
These are both the adjacency matrices of trees, and both have characteristic polynomial
$$\lambda^{11}-10\lambda^9+34\lambda^7-47\lambda^5+25\lambda^3-4\lambda.$$
Each tree
has exactly two vertices of degree 3, separated by a path of length 1 in the case of $A$ but length 2 in the case of $B$. In particular, the trees are not isomorphic.

Now consider the [EDIT: improved, much nicer] matrix
$$ C = \left[\begin{array}{rrrrrrrrrrr}
1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 \\\\
0 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\\\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 \\\\
0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 \\\\
\end{array}\right] $$
with determinant $-1$.

Since $C^{-1}AC = B$, the two trees (on 11 vertices) are non-isomorphic but
have adjacency matrices that are conjugate over $\mathbb Z$.

Now to explain the where the example comes from. The pair of graphs
was constructed by a method, attributed to Schwenk, that I found in Doob's
chapter of *Topics in algebraic graph theory* (edited by Beineke and Wilson).
The first 9 rows and columns of $A$, in common with $B$, come from a particular tree on 9 vertices that has a pair of attachment points such that extending the tree in the same way from either point gives isomorphic spectra.
Adding a single pendant vertex cannot work for this problem, as I found using Brouwer and van Eijl's trick, mentioned
by Chris Godsil, of comparing the Smith normal forms of (very) small polynomials in $A$ and $B$, in this case $A+2I$ and $B+2I$. When a path of length two is added at either of the two special vertices, however, there doesn't seem to be any obstruction of this type.

I then set about trying to conjugate both $A$ and $B$, separately, to the companion matrix of their mutual characteristic polynomial, by looking for a random small integer vector $x$ for which the matrix $X_A = [ x\ Ax\ A^2x\ \ldots\ A^{10}x]$ has determinant $\pm 1$, and similarly $y$ giving $Y_B$. (The fact that I succeeded fairly easily may have something to do with the fact that $A+I$ is invertible over $\mathbb Z$.) The matrix
$X_AY_B^{-1}$ then acts like the $C$ above.

[EDIT: The actual matrix $C$ I found at random and first posted was not nearly so pretty, with a Frobenius norm nearly ten times the current example. But taking powers 0 to 10 of $A$ times $C$ gave a $\mathbb Q$-basis for the full space of conjugators, whose Smith normal form (as 11 vectors in $\mathbb R^{121}$) was all 1's down the diagonal, so in fact it was a $\mathbb Z$-basis. Performing an LLL reduction on this lattice basis then gave a list of smaller-norm matrices, the third of which is the more illuminating $C$ given above, of determinant $-1$. The other determinants from the reduced basis were all $0$ and $\pm 8$.]

Taking rational $x$ and not restricting the determinant of $X_A$ gives a space of possible rational matrices $C$ of dimension 11, which are generically invertible; varying $y$ gives the same space [EDIT: as does multiplying on the left by powers (or in the more general case commutants) of $A$]. Since the spectrum of $A$ has no repeated roots, this is also the dimension of the commutant of $A$, and every matrix conjugating $A$ to $B$ lies in this space. Starting with a rational basis, it is not hard to find an exact basis for the integer sublattice, and taking the determinant of a general point in the integer lattice gives an integer polynomial in 11 variables which takes the value $1$ or $-1$ if and only if the matrices $A$ and $B$ are conjugate over $Z$. If there are repeated roots, you have to work a little harder; in general the full space has dimension the sum of the squares of the multiplicities, and is generated by multiplying on the left by a basis for the commutator space of $A$. A basis for the commutant can be produced (for a diagonalizable matrix) by first conjugating $A$ to a direct sum of companion matrices for the irreducible factors of the characteristic polynomial, and then one at a time, for each $k$-by-$k$ block corresponding to a $k$-times repeated factor of degree $m$, replacing each of the $k^2$ blocks with powers $0$ to $m-1$ of the companion matrix for that factor, with $0$ everywhere elsewhere.