Some other things that are true:
A square matrix of full rank has linearly independent columns.
A square matrix that can be row-reduced to the identity matrix has full rank.
The matrix for taking the standard basis to the standard basis is the identity matrix.
Row-reduction operations are reversible, and can be represented by invertible matrices. For example, in the $3 \times 3$ case, the row-operation that adds twice the third row to the second is left-multiplication by the matrix:
$\begin{bmatrix}1&0&0\\0&1&2\\0&0&1\end{bmatrix}$
which has the inverse which subtracts twice the third row form the second:
$\begin{bmatrix}1&0&0\\0&1&-2\\0&0&1\end{bmatrix}$
The net effect of taking to row-reduced echelon form is to left-multiply a matrix $A$ by an invertible matrix $P$.
Now if $PA = I$, multiplying by $P^{-1}$ on the left gives us:
$P^{-1}(PA) = P^{-1}I$
$(P^{-1}P)A = P^{-1}$
$IA = P^{-1}$
$A = P^{-1}$
Hence $A^{-1} = (P^{-1})^{-1} = P$
So, if $A$ can be row-reduced to the identity, then $A$ is invertible, which means in particular that its range (as the function which takes a vector $v$ to $Av$) is all of $\Bbb R^n$ (that is, $A$ is surjective, or onto, as a function). This means that its columns (which span its range) must be linearly independent, since we have a spanning set of $n$ vectors and the dimension of $\Bbb R^n$ is $n$.
On the other hand, suppose the columns of an $n \times n$ matrix $A$ are linearly independent. Then these form a basis for $\Bbb R^n$, since $A$ has $n$ columns.
If we label these columns as $u_1,\dots,u_n$, linear independence means that:
$c_1u_1 +\cdots + c_nu_n = 0 \implies c_1 = \dots = c_n = 0$.
Since any vector $Av$ is in the column space spanned by the $u_i$, this means the ONLY solution to $Av = 0$ is $v = 0$.
Recall that the "general solution" of $Av = b$, is:
$v = x_0 + x$, where $x_0$ is one PARTICULAR solution, and $x$ is ANY solution of $Av = 0$.
For our particular $A$, we know there is only one solution of $Av = 0$, namely $v = 0$. This means there is only one particular solution $x_0$ (for if there were two, say $x_0,x_1$, with $x_0 \neq x_1$, we would have:
$Ax_0 = b$ and $Ax_1 = b$ so that $Ax_0 - Ax_1 = b - b = 0$, and thus $A(x_0 - x_1) = 0$, meaning $x_0 - x_1 = 0$, by the uniqueness of this solution. But that means $x_0 = x_1$, contradicting that we have two solutions).
Now since the columns of $A$ span $\Bbb R^n$ (being a LI set of cardinality $n$), we always have at least ONE solution to $Av = b$, and now we see we have EXACTLY one.
Now let's look at what this means for row-reduction: in particular for the leading 1's of our rows. These always move "to the right" in successive rows.
Claim 1: each leading 1 is exactly "one right" of the row above. For if not, as we go down, we would run out of room "to go right" before we reached the last row, leading to at least one zero row at the bottom. Such a zero row would lead to:
$PAx = \begin{bmatrix}1&\ast&\dots&\ast\\0&\ast&\dots&\ast\\ \vdots&\vdots&\ddots&\vdots\\0&0&\dots&0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\ \vdots\\x_n\end{bmatrix} = \begin{bmatrix}y_1\\y_2\\ \vdots\\0\end{bmatrix}$
In particular, there is NO $x \in \Bbb R^n$ with $PAx = e_n$.
Why is this a problem? Well, since $P$ is invertible, there is no $x \in \Bbb R^n$, with $Ax = P^{-1}e_n$ (for if there were, then $PAx = P(P^{-1}e_n) = Ie_n = e_n$, contradiction).
But EVERY vector in $\Bbb R^n$ is $Av$ for some $v$ (because the column space spans $\Bbb R^n$ since the columns are a basis), including whatever $P^{-1}e_n$ is. So such a zero row cannot happen.
Claim 1 means our row-reduced echelon form is:
$PA = \begin{bmatrix}1&\ast&\dots&\ast\\0&1&\dots&\ast\\ \vdots&\vdots&\ddots&\vdots\\0&0&\dots&1\end{bmatrix}$
Since we are in row-reduced echelon form, all the $\ast$ above each one are zeros:
$PA = I$.
In other words, $A$ is an $n \times n$ matrix of basis vectors if and only if the row-reduced echelon form of $A$ is the identity matrix, as you indeed conjectured.
As an aside, you remark that "row-reduction doesn't change the space". It doesn't change the SOLUTION space, it can change the COLUMN space. The first is the elements that $A$ acts on, the second is the RESULTS. What NEVER changes is the DIMENSION of these spaces, and if $A$ is ONTO (if it indeed has columns that form a basis), then the column space doesn't change since "all of $\Bbb R^n$" is still all of $\Bbb R^n$ no matter how you are slicing and dicing it.
If $A$ isn't onto (if it row-reduces to something other than the identity), then to get the column space of $A$ you have to reference the "pivot columns" of the row-reduced matrix in the original $A$, not the columns of the row-reduced matrix.
Aside number 2: linear transformations are ALGEBRA, matrices are ARITHMETIC. We can "mirror" the algebraic structure of linear transformations in the arithmetic of matrices (which use actual numbers we can add, and multiply, etc.), but to do so, we have to "pick" a basis, and this choice can be done in "more than one way". Mathematicians prefer, in general, proofs that don't require such a choice of basis, because these show the underlying structure inherent in just the "vectorness of vectors", whereas numerical or algorithmic arguments depend on our basis coordinates, which are somewhat arbitrarily chosen (Imagine modelling vector forces in physical space. The choice of $x,y,z$ axes and the origin's location aren't part of the physical problem, they are something we IMPOSE upon the situation to make life easier for us. Sometimes, halfway through this process, we realize a different choice makes for easier calculations-the vectors (the underlying algebra) remain the same, but the coordinates (the numbers we calculate with) can be totally different).
