TRANSPOSE and MATMUL interaction

I am confused by the unexpected results from the following:

SUBROUTINE JacobianFwdDiff(m,n,fjac,ldfjac)

IMPLICIT NONE

INTEGER, INTENT(IN) :: m,n,ldfjac
REAL(8), DIMENSION(ldfjac,n), INTENT(OUT) :: fjac
REAL(8), DIMENSION(n,n) :: term1,M3,M4
REAL(8), DIMENSION(n,m) :: M1
REAL(8), DIMENSION(m,n) :: M2

! Code which computes fjac(1:m,1:n)

M1=TRANSPOSE(fjac(1:m,1:n)); M2=fjac(1:m,1:n)
M3=MATMUL(M1, M2); M4=MATMUL(M1(1:n,1:m), M2(1:m,1:n))
term1(1:n,1:n)=MATMUL(TRANSPOSE(fjac(1:m,1:n)), fjac(1:m,1:n))

Specifically, when m=ldfjac=769, n=1, the values for M3 and M4 are
identical, as I expected, but the value for term1 differs from these.
I thought the computation of term1 should be the same as that of M4,
but it is not.

What have I got wrong here?

.