We present the first simulations of the tidal disruption of stars with realistic structures and compositions by massive black holes (BHs). We build stars in the stellar evolution code MESA and simulate their disruption in the 3D adaptive-mesh hydrodynamics code FLASH, using an extended Helmholtz equation of state and tracking 49 elements. We study the disruption of a 1$M_odot$ star and 3$M_odot$ star at zero-age main sequence (ZAMS), middle-age, and terminal-age main sequence (TAMS). The maximum BH mass for tidal disruption increases by a factor of $sim$2 from stellar radius changes due to MS evolution; this is equivalent to varying BH spin from 0 to 0.75. The shape of the mass fallback rate curves is different from the results for polytropes of Guillochon & Ramirez-Ruiz (2013). The peak timescale $t_{rm peak}$ increases with stellar age, while the peak fallback rate $dot M_{rm peak}$ decreases with age, and these effects diminish with increasing impact parameter $beta$. For a $beta=1$ disruption of a 1$M_odot$ star by a $10^6 M_odot$ BH, from ZAMS to TAMS, $t_{rm peak}$ increases from 30 to 54 days, while $dot M_{rm peak}$ decreases from 0.66 to 0.14 $M_odot$/yr. Compositional anomalies in nitrogen, helium, and carbon can occur before the peak timescale for disruptions of MS stars, which is in contrast to predictions from the frozen-in model. More massive stars can show stronger anomalies at earlier times, meaning that compositional constraints can be key in determining the mass of the disrupted star. The abundance anomalies predicted by these simulations provide a natural explanation for the spectral features and varying line strengths observed in tidal disruption events.