One way to think about it is that orthogonalize! makes it so that you can just use the tensor on the site you orthogonalize to in order to form the local density matrix. For example, the local density matrix on site j would be:
using ITensors
n = 10
s = siteinds("S=1/2", n)
psi = randomMPS(s)
j = 4
orthogonalize!(psi, j)
# Without calling `orthogonalize!`, this would not be
# the local density matrix! We would have to trace out
# the sites that are not `j`.
rho_j = prime(dag(psi[j]), siteinds(psi, j)) * psi[j]
# Compute an expectation value
(rho_j * op("Sz", s, j))[]
# or:
(prime(dag(psi[j]), siteinds(psi, j)) * op("Sz", s, j) * psi[j])[]
# Same as:
expect(psi, "Sz"; sites=j)
Besides computing local expectation values, orthogonalizing and forming the local density matrix is helpful for a variety of other applications, such as truncating an MPS. For example, internally in the truncate! function, we orthogonalize the MPS to a certain site, take an SVD of the local MPS tensor (which is equivalent to diagonalizing the local density matrix), and then truncate according to the singular values. Without orthogonalizing, the truncation would not be optimal.
I would recommend review articles such as Section 4 of https://arxiv.org/pdf/1008.3477.pdf for more background on MPS orthogonality.