I think you would want the second command to be this:

3dAllineate \
    -master ${anat}  \
    -input ${maskdir}/language_comprehension_fromNeuroSynth.nii \
    -prefix ${wdir}/WernickeinNative_NEW \
    -1Dmatrix_apply ${wdir}/transformation_matrix.aff12.1D
    -final NN

.... where note the use of the option "-master .." instead of "-base .."---this makes the output grid of the "-prefix .." dset match ${anat}, rather than just be resampled at input/source dset resolution. Also note that I added "-final NN". Since you said earlier that the input/source dset was a mask, this would have it remain integer-valued; here is part of the 3dAllineate help for more detail:

-final iii  = Defines the interpolation mode used to create the
               output dataset.  [Default == 'cubic']
            ** N.B.: If you are applying a transformation to an
                       integer-valued dataset (such as an atlas),
                       then you should use '-final NN' to avoid
                       interpolation of the integer labels.

Again, if you are happy with the 3dAllineate-estimated alignment for your purposes, that is fine. Just note for other cases where more precise alignment is needed, then using nonlinear warping (e.g., 3dQwarp or @SSwarper, for example) would be recommended.

--pt