RNA-based thermal regulation mechanism is an important strategy for organisms to cope with temperature changes. Inhabiting on intertidal rocky shore, a key interface of the ocean, atmosphere and terrestrial environments, intertidal species have developed variable thermal adaptation mechanisms, however, adaptions at the RNA level remain largely uninvestigated. For investigating the relationship between mRNA structural stability and species distribution, in the present study, the secondary structures of cytosolic malate dehydrogenase (cMDH) mRNAs of Echinolittorina malaccana, E. radiata and Littorina brevicula were determined using the selective 2'-hydroxyl acylation analyzed by primer extension (SHAPE), and the change in folding free energy of formation (ΔGfold) were calculated. The results showed that ΔGfold increased as the temperature increased. The difference in ΔGfold (ΔΔGfold) between two specific temperatures (0° to 25°C, 0° to 37°C, and 0° to 25°C) among the three species were different, and the ΔΔGfold value of E. malaccana was significantly lower than those of E. radiata and L. brevicula. The number of stems of cMDH mRNA of the snails decreased with temperature increased, and the breakpoint temperature of E. malaccana was the highest among them. The number of loops was also reduced with temperature increased, while the length of the loop structure increased accordingly. Consequently, these structural changes can potentially affect the translational efficiency of mRNA. These results implied that there were interspecific differences in the thermal stability of RNA secondary structures in intertidal snails, and these differences may be related to their distribution.