This was detected using a nanoindentation instrument (Hysitron TI980 TriboIndenter, Bruker Corporation) and its association between the displacement and the experimental load was analyzed.

The TPC-1 and KTC-1 cell lines, both derived from humans, were purchased from the Cell Bank of the Chinese Academy of Sciences. Both the cell lines were maintained in RPMI-1640 medium (Invitrogen; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Invitrogen; Thermo Fisher Scientific, Inc.), 100 IU/ml penicillin (Gibco; Thermo Fisher Scientific, Inc.) and 100 µg/ml streptomycin (Gibco; Thermo Fisher Scientific, Inc.), at 37°C in a humidified incubator with 5% CO₂.

The cells were then seeded in a 6-well plate, at a concentration of 1×10⁵ cells/well with 10, 40 and 60 kPa substrate for 12 h, in which the substrate was prepared according to a previous report (13). Subsequently, the cells were scratched vertically with a 100-µl pipette tip the following day. At 48 h, the cells were counted under an inverted phase-contrast microscope in 5 random fields, after washing twice with PBS and placed in serum-free culture medium.

The cells of 1×10⁵/ml were seeded in six-well culture plates with different substrate stiffness. The 10 kPa stiffness group was served as the control. After 24-h cell culture, cells were trypsinized for cell counting using Coulter counter every 12 h, then re-suspended and reseeded in dishes till 72 h. The above experiments were repeated three times for each group.

Computed tomography (CT) data of the neck was collected from 7 patients with PTC for FEA (Department of Maxillofacial and Ear Nose and Throat Oncology, Tianjin Medical University Cancer Institute and Hospital). The 3D image data of the patients with PTC were captured using CT technology (ScanXmate-E090; Comscantecno Co., Ltd.). Following air calibration, X-ray exposure in the neck of each patient was performed with a view field of 25×25 cm and an axial scan, at a thickness of 2.5×2.5 mm interval. The CT data was used to construct the 3D polygonal stereolithography (STL) model of the neck in each case. A part of the STL model, with ~2.5 mm anterior-posterior thickness, was retrieved from the whole model to obtain a segment representing the region of interest, including the loading site. This segmented STL model was converted into computer-aided design software (Catia V5R18; Dassault Systems) to analyze the model in detail and examine minute irregularities. Finally, the 3D volume model was constructed using Ansys finite element software (Ansys version 11.0; ANSYS, Inc.), which was meshed by 10-nodes quadratic tetrahedral element with 3 degrees of freedom.

The thyroid gland was set as viscoelastic and isotropic material, with an initial Young's modulus of 34.85 MPa and Poisson's ratio of 0.49 (14). The trachea was defined as isotropic elastic material with a Young's modulus of 3.33 MPa and Poisson's ratio of 0.49 (15). The density of the thyroid was set as 1,150 kg/m³, and the trachea was considered as cartilage only with a density of 1,400 kg/m³, which was obtained using the density measurement function of Mimics version 8.1 software (Materialise; http://www.materialise.com/).

The present clinical study included 998 patients who were diagnosed with PTMC by a pathologist following surgery at the Tianjin Medical University Cancer Hospital (TJMUCH) between 1st June 2016 and 1st December 2016; 709 patients had a single tumor and 289 patients had multifocal disease. The mean age was 48.2±11.3 years and there were 261 males and 737 females. The present study was conducted according to the principles outlined in the Helsinki Declaration and was approved by the Research Ethics Committee at TJMUCH (no. 2018090). Written informed consent was provided by all the patients. The clinicopathological data, including sex, age and the presence of thyroiditis, were also collected. All the patients were subjected to total thyroidectomy or unilateral thyroidectomy according to the National Comprehensive Cancer Network guidelines (version 2016) (16,17); cases in which the tumor was close to the tracheal region were marked as ‘interior’, and those with tumors 2–4 cm away from the trachea were marked as ‘lateral’.

A total of 33 samples, each containing PTMC located in interior regions and adjacent normal tissue, were collected for RNA-seq analysis (3 samples) and RT-qPCR (30 samples). A total of 10 samples, including multifocal PTMC located in the interior and lateral regions, were collected for RT-qPCR for comparing differences in gene expression between PTMC located in different regions of the trachea. All tissue samples were obtained by thyroidectomy and stored at −80°C until further use. RNA sequencing (RNA-seq) and subsequent analysis was performed by the BGI-Shenzhen Company (http://www.genomics.cn/en/).

Gene expression levels of the transcripts were measured using the reads per kilobase of transcript, per million mapped reads method. Subsequently, the edgeR package (edgeR 3.14.0) (18) tool was utilized to identify the DEGs between the PTMC and adjacent normal tissue from 3 samples of PTMC. During the differential analysis, the negative binomial model was used to calculate the significance of the differentially expressed mRNAs, followed by the adjustment of P-values using the Benjamini Hochberg method (19). The cut-off values of the DEG selection was a false discovery rate adjusted P<0.05 and |log₂ fold change |≥1. These results were determined based on the comparison with The Cancer Genome Atlas (TCGA). Venn diagram analysis was performed using online software at the following URL: http://bioinformatics.psb.ugent.be/webtools/Venn/.

To further investigate the functions and pathways of the DEGs, Gene Ontology (GO; http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (www.kegg.jp) pathway enrichment analyses were performed using ‘Term Finder’ (https://yeastgenome.org/goTermFinder), with a threshold corrected P≤0.05 for the identification of significant GO terms and pathways using default settings. Furthermore, Search Tool for the Retrieval of Interacting Genes/Proteins (http://string-db.org/) was used to further investigate the associations in the DEGs, at the protein level (20). The criterion for the construction of the PPI network was based on the confidence score ≥0.90. The Cfinder software (version 2.0.6, http://www.cfinder.org/) was used to extract the functional modules of the PPI network with default parameters (21).

30 samples tissues, including PTMC and adjacent normal tissue, were selected for RT-PCR test. Three biomechanical genes, which were overexpressed from RNA-seq analysis were analyzed using RT-qPCR for validation and the TransStart Top Green qPCR SuperMix (Beijing Transgen Biotech Co., Ltd.). The PCR was performed using the ABI7500 Real-Time qPCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.) and the following thermocycling conditions: Initial denaturation at 95°C for 20 sec, followed by 40 cycles of 95°C for 15 sec, annealing at 60°C for 30 sec and extension at 70°C for 90 sec. The specificity of RT-qPCR was examined using the dissociation curve and the relative expression of the selected DEGs was normalized with the 18S rRNA gene. The cycle threshold (Cq) 2^−∆∆Cq method was used to calculate the relative expression level (22). All the gene specific primers used are listed in Table SI.

All the data are presented as the mean ± SEM. Gene expression analysis was based on sex, age or thyroiditis. For the data from the patients with PTC, the categorical variables were analyzed using a χ² or Fisher's exact test. Continuous variables were analyzed using a paired and unpaired Student's t-tests as appropriate. Comparisons of >2 groups were performed using ANOVA followed by Bonferroni's post hoc test. P<0.05 was considered to indicate a statistically significant difference. ImageJ software (version 1.42; National Institutes of Health) was used for wound healing analysis. Origin software (version 9.0; OriginLab Corporation) was used for all other data analysis.

Thyroid images were captured and used to construct thyroid model (Fig. S1). When the load reaches the maximum value (hmax), the displacement also reaches the maximum value (Pmax), that is the maximum indentation depth. After unloading, the displacement finally returns to a fixed value (S). At this time, the depth is called the residual indentation depth (hr), that is the permanent plastic deformation left by the indenter on the sample (Fig. S2A). As the load reaches the maximum value of 0.743 mN, the displacement reached the maximum indentation depth of 984.63 nm; after unloading, the displacement finally returned to the residual indentation depth of 881.36 nm (Fig. S2C). The results showed that PTMC tissue stiffness ranged from 20–70 kPa due to calcification or fibrosis (Fig. S2D).

To evaluate the migratory and proliferative abilities of the PTC cells on different substrates, these characteristics were analyzed using different experiments. As shown in Fig. 1, the scratch areas in the control and 60 kPa groups were significantly narrower compared with that in the 10 and 40 kPa groups for the TPC-1 (Fig. 1A) and KTC-1 (Fig. 1B) cell lines, suggesting a stronger migratory ability of these cells on the stiffer substrate compared with that on the softer substrate (P<0.001; Fig. 1C and D). Furthermore, the density of the TPC-1 cells on the 60 kPa substrate following 72 h of culture was 7.19±0.2×10⁵/ml. However, the cell densities were 2.88±0.06×10⁵/ml and 1.2±0.03×10⁵/ml on the 40 and 10 kPa substrates, respectively (Fig. 1E). It was notably different among the different stiffness groups compared with 10 kPa groups, (P<0.001; Fig. 1F). For the KTC-1 cell line, similar results were obtained, with densities of 7.36±0.29×10⁵/ml, 4.59±0.11×10⁵/ml and 1.75±0.05×10⁵/ml on the 60, 40 and 10 kPa substrates, respectively (Fig. 1G). A statistically significant difference was observed among these groups (P<0.001; Fig. 1H), which indicated a stronger proliferation ability of the cells on the stiffer substrate compared with that on the softer substrate.

Von Mises stress, an equivalent stress, was selected to reveal the distribution of elastic stress in the biological tissue. As shown in Fig. 2C, the von Mises stress in the thyroid was mainly produced by stiffness and the deformation of the trachea, which was significantly higher in the ‘interior’ compared with that in the ‘lateral’ regions (10.18±3.35 vs. 7.71±2.89 kPa, respectively; P<0.05).

As shown in Table I, the percentages of PTMC distribution in the ‘interior’ and ‘lateral’ sections were 68.58 and 31.42% in males, and 66.22 and 33.78% in females, respectively. There were no statistically significant differences in the parameters of sex, age or thyroiditis (P>0.05). The number of PTMC cases in the ‘interior’ and ‘lateral’ sections was 916 and 458, respectively, suggesting a possible association between PTMC distribution and Von Mises stress in the thyroid.

To evaluate biomechanical gene expression in different tissues, including PTMC and adjacent normal tissue, a total of 18,154 genes were detected from gene expression profiling. Flow diagram analysis of all the expressed genes indicated that 15,743 (86.7%) were expressed in all 3 samples for RNA-seq (Fig. 3B). The heatmap revealed the relative gene expression level in different samples (Fig. 3B). Based on the criteria, 1,504 genes were identified as DEGs between the PTC and adjacent normal tissues, of which 794 were upregulated and 710 were downregulated (Fig. S4A). In addition, the results of RT-qPCR indicated that the overexpressed genes in PTC were biomechanical sensors, including Piezo2, TRPV4 and CDH3 (Fig. 4A), ECM stiffness and EMT related genes, including CLDN1, CLDN16, Runx2, Twist and G3BP3 (Fig. 4B) compared with that in the lateral normal tissues. However, some of the G-protein-related genes, including ADORA1 and GABBR2, exhibited a higher expression in PTC compared with that in adjacent normal tissues; however, there was no notable difference in the expression of RHOA among the groups (Fig. 4C). The expression level of these genes did not exhibit any significant difference by sex, age or thyroiditis, as shown in Fig. 5.

With respect to the Piezo2, CHD3, Runx2 and Twist1 genes, the expression levels in PTMC, located in the interior section were notably higher compared with that in the lateral regions (Fig. S3).

GO functional enrichment analysis in PTC, in comparison with TCGA, revealed the top 10 GO enrichment terms for tumors for ‘biological processes’ (BP) (Fig. 4D), ‘cellular component’ (CC) (Fig. 4E) and ‘molecular function’ (MF) (Fig. 4F).

The upregulated genes were significantly enriched in ‘cell adhesion’ (GO: 0007155), ‘immune system process’ (GO: 0002376), ‘ECM organization’ (GO: 0030198), ‘osteoclast formation’, ‘anatomical structure development’ (GO: 0048856), ‘single-organism cellular process’ (GO: 0044763), ‘response to stress’ (GO: 0006950), ‘cell proliferation’ (GO: 0008283) (Fig. S4).

The most notable BP terms of the downregulated genes were ‘anatomical structure development’ (GO: 0048856), ‘single-organism developmental process’ (GO: 0044767), ‘cell differentiation’ (GO: 0030154), ‘cellular developmental process’ (GO: 0048869), ‘extracellular matrix organization’ (GO: 0030198), ‘cell morphogenesis’ (GO: 0000902), ‘cell adhesion’ (GO: 0007155) (Fig. S4C). KEGG pathway annotation revealed that signal transduction and immune system-related genes mainly participated in PTMC (Fig. S4D).